title: Arrange monte carlo p-value into a matrix for different sample size and variance estimators
body: <p>The following code works out quite well (based on <a href="https://stackoverflow.com/q/40674367/4891738">my previous question</a>). But I have to change the variance estimator (<code>ols</code>, <code>hc0</code>, <code>hc1</code>, <code>hc2</code>, <code>hc3</code>) every time before I run the code. I would like to solve this problem with a loop.</p>
  
  <p>Hereafter, I briefly describe the code. Within the code, 1000 regression models for each sample size (<code>n = 25, 50, 100, 250, 500, 1000</code>) are created. Then, each regression model out of the 1000 is estimated by OLS. After that, I calculate t-statistics based on the different beta values of <code>x3</code> out of the 1000 samples. The null hypothesis reads: <code>H0: beta03 = beta3</code>, that is the calculated beta value of <code>x3</code> equals the 'real' value which I defined as 1. In the last step, I check how often the null hypothesis is rejected (significance level = 0.05). My final goal is to create a code which spits out the procentual rejection rate of the null hypothesis for each sample size and variance estimator. Thus, the result should be a matrix whereas right now I get a vector as a result. I would be pleased if anyone of you could help me with that. Here you can see my code:</p>
  
  <pre><code>library(car)
  sample_size = c("n=25"=25, "n=50"=50, "n=100"=100, "n=250"=250, "n=500"=500, "n=1000"=1000)
  
  B &lt;- 1000
  beta0 &lt;- 1
  beta1 &lt;- 1
  beta2 &lt;- 1
  beta3 &lt;- 1
  alpha &lt;- 0.05
  
  simulation &lt;- function(n, beta3h0){
  t.test.values &lt;- rep(NA, B)
  #simulation of size
  for(rep in 1:B){
  #data generation 
  d1 &lt;- runif(n, 0, 1)
  d2 &lt;- rnorm(n, 0, 1)
  d3 &lt;- rchisq(n, 1, ncp=0)
  x1 &lt;- (1 + d1)
  x2 &lt;- (3*d1 + 0.6*d2)
  x3 &lt;- (2*d1 + 0.6*d3)
  # homoskedastic error term: exi &lt;- rchisq(n, 4, ncp = 0)
  exi &lt;- sqrt(x3 + 1.6)*rchisq(n, 4, ncp = 0)
  y &lt;- beta0 + beta1*x1 + beta2*x2 + beta3*x3 + exi
  mydata &lt;- data.frame(y, x1, x2, x3)
  #ols estimation
  lmobj &lt;- lm(y ~ x1 + x2 + x3, mydata)
  #extraction
  betaestim &lt;- coef(lmobj)[4]
  betavar   &lt;- vcov(lmobj)[4,4]
  #robust variance estimators: hc0, hc1, hc2, hc3
  betavar0 &lt;- hccm(lmobj, type="hc0")[4,4]
  betavar1 &lt;- hccm(lmobj, type="hc1")[4,4]
  betavar2 &lt;- hccm(lmobj, type="hc2")[4,4]
  betavar3 &lt;- hccm(lmobj, type="hc3")[4,4]
  #t statistic
  t.test.values[rep] &lt;- (betaestim - beta3h0)/sqrt(betavar)
  }
  mean(abs(t.test.values) &gt; qt(p=c(1-alpha/2), df=n-4))
  }
  
  sapply(sample_size, simulation, beta3h0 = 1)
  </code></pre>