targa targa - 18 days ago 5
R Question

How to create a double loop?

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

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.

Comments