targa - 1 year ago 73

R Question

The following code works out quite well, BUT: I have to change the variance estimator (

`ols`

`hc0`

`hc1`

`hc2`

`hc3`

Hereafter, I briefly describe the code. Within the code, 1000 regression models for each sample size (

`n = 25, 50, 100, 250, 500, 1000`

`x3`

`H0: beta03 = beta3`

`x3`

`library(car)`

sample_size = c("n=25"=25, "n=50"=50, "n=100"=100, "n=250"=250, "n=500"=500, "n=1000"=1000)

B <- 1000

beta0 <- 1

beta1 <- 1

beta2 <- 1

beta3 <- 1

alpha <- 0.05

simulation <- function(n, beta3h0){

t.test.values <- rep(NA, B)

#simulation of size

for(rep in 1:B){

#data generation

d1 <- runif(n, 0, 1)

d2 <- rnorm(n, 0, 1)

d3 <- rchisq(n, 1, ncp=0)

x1 <- (1 + d1)

x2 <- (3*d1 + 0.6*d2)

x3 <- (2*d1 + 0.6*d3)

# homoskedastic error term: exi <- rchisq(n, 4, ncp = 0)

exi <- sqrt(x3 + 1.6)*rchisq(n, 4, ncp = 0)

y <- beta0 + beta1*x1 + beta2*x2 + beta3*x3 + exi

mydata <- data.frame(y, x1, x2, x3)

#ols estimation

lmobj <- lm(y ~ x1 + x2 + x3, mydata)

#extraction

betaestim <- coef(lmobj)[4]

betavar <- vcov(lmobj)[4,4]

#robust variance estimators: hc0, hc1, hc2, hc3

betavar0 <- hccm(lmobj, type="hc0")[4,4]

betavar1 <- hccm(lmobj, type="hc1")[4,4]

betavar2 <- hccm(lmobj, type="hc2")[4,4]

betavar3 <- hccm(lmobj, type="hc3")[4,4]

#t statistic

t.test.values[rep] <- (betaestim - beta3h0)/sqrt(betavar)

}

mean(abs(t.test.values) > qt(p=c(1-alpha/2), df=n-4))

}

sapply(sample_size, simulation, beta3h0 = 1)

`exi <- sqrt(x3 + 1.6)*rchisq(n, 4, ncp = 0)`

`exi <- rchisq(n, 4, ncp = 0)`

Recommended for you: Get network issues from **WhatsUp Gold**. **Not end users.**

Answer Source

You don't need a double nested loop. Just make sure you get a matrix inside your loop. Update your current `simulation`

with the following:

```
## set up a matrix
## replacing `t.test.values <- rep(NA, B)`
t.test.values <- matrix(nrow = 5, ncol = B) ## 5 estimators
## update / fill a column
## replacing `t.test.values[rep] <- (betaestim - beta3h0)/sqrt(betavar)`
t.test.values[, rep] <- abs(betaestim - beta3h0) / sqrt(c(betavar, betavar0, betavar1, betavar2, betavar3))
## row means
## replacing `mean(abs(t.test.values) > qt(p=c(1-alpha/2), df=n-4))`
rowMeans(t.test.values > qt(1-alpha/2, n-4))
```

`simulation`

would return a vector of length 5. For each sample size, the "rejection rate" of t-statistic is returned for all 5 variance estimators. Then, when you call `sapply`

, you get a matrix result:

```
sapply(sample_size, simulation, beta3h0 = 1)
# n=25 n=50 n=100 n=250 n=500 n=1000
#[1,] 0.132 0.237 0.382 0.696 0.917 0.996
#[2,] 0.198 0.241 0.315 0.574 0.873 0.994
#[3,] 0.157 0.220 0.299 0.569 0.871 0.994
#[4,] 0.119 0.173 0.248 0.545 0.859 0.994
#[5,] 0.065 0.122 0.197 0.510 0.848 0.993
```

Recommended from our users: **Dynamic Network Monitoring from WhatsUp Gold from IPSwitch**. ** Free Download**