targa - 9 months ago 38

R Question

I would like to know how to create a double loop.

In my code, I do a multiple regression on 1000 samples (each sample size: 25)

Then, I create t test values for each sample out of the 1000 with the nullhypothesis: beta3 value from the sample = 'real' beta3 value. I know the 'real' beta3 value from a Monte Carlo Simulation (beta 3 = value of third coefficient of the regression).

However, the code works so far.

Now, I want to do the same procedure for sample sizes of 50, 100, 250, 500, and 1000 (each sample size 1000 times).

How can I realise this goal with a loop. I would be pleased if you could help me! Here you can see my code:

`n <- 25`

B <- 1000

beta3 <- 1.01901 #'real' beta3 value

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

for(rep in 1:B){

##data generation

d1 <- runif(25, 0, 1)

d2 <- rnorm(25, 0, 1)

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

x1 <- (1 + d1)

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

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

exi <- rchisq(25, 5, ncp = 0)

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

## estimation

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

## extraction

betaestim <- coefficients(lmobj)[2:4]

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

## t-test

t.test.values[rep] <- (betaestim[3] - beta3)/sqrt((betavar)[9])

}

Answer Source

We can use a `data.frame`

to store the results. Also note that you didn't include values for `beta0`

, `beta1`

, or `beta2`

, so I just used placeholder values.

```
n <- c(50,100,250,500,1000) #how big are our sample sizes?
B <- 1000
beta3 <- 1.01901 #'real' beta3 value
#other beta values (note that these were not included in your question)
beta1 <- 2
beta2 <- 4
beta0 <- 6
iter <- 1
#initialize our results data.frame
result_df <- data.frame(sample_size = numeric(length(n) * B),
t.test.values = numeric(length(n) * B)
)
for(size in n){
for(rep in 1:B){
##data generation
d1 <- runif(size, 0, 1)
d2 <- rnorm(size, 0, 1)
d3 <- rchisq(size, 1, ncp=0)
x1 <- (1 + d1)
x2 <- (3 * d1 + 0.6 * d2)
x3 <- (2 * d1 + 0.6 * d3)
exi <- rchisq(size, 5, ncp = 0)
y <- beta0 + beta1*x1 + beta2*x2 + beta3*x3 + exi
## estimation
lmobj <- lm(y ~ x1 + x2 + x3)
## extraction
betaestim <- coefficients(lmobj)[2:4]
betavar <- vcov(lmobj)[2:4, 2:4]
## store our values
result_df[iter, 1] <- size
result_df[iter, 2] <- (betaestim[3] - beta3)/sqrt((betavar)[9])
iter = iter + 1 #iterate
}
}
```

As long as you use something to track the iterations (I've used `iter`

here), a double `for`

loop isn't too bad. Just make sure to initialize a `data.frame`

to the correct size. It might be helpful to look at the `replicate`

function and the `*apply`

class of functions if you're planning to do more simulations.