さまざまなサンプルサイズと分散推定値の行列にモンテカルロP値を配置します

Question

次のコードは非常にうまくいきます（my previous questionに基づいています）。しかし、コードを実行する前に、分散推定量（ols、hc0、hc1、hc2、hc3）を変更する必要があります。私はループでこの問題を解決したいと思います。

以下、コードについて簡単に説明します。コード内で、各標本サイズ（n = 25, 50, 100, 250, 500, 1000）の1000個の回帰モデルが作成されます。次に、1000のうちの各回帰モデルがOLSによって推定される。その後、1000サンプルのうちの異なるベータ値に基づいてt統計量を計算します。帰無仮説はH0: beta03 = beta3と計算され、の計算ベータ値は1と定義された「実数」の値と等しくなります。最後のステップでは、帰無仮説が棄却される頻度（有意水準= 0.05）を確認します。私の最終的な目標は、各標本サイズと分散推定量について帰無仮説の棄却率を吐き出すコードを作成することです。したがって、結果は行列になりますが、今は結果としてベクトルが得られます。あなたの誰かが私にそれを手伝ってもらえれば嬉しいです。ここで私のコードを見ることができます：

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) 

Answer 1

ダブルネストループは必要ありません。ループの中に行列があることを確認してください。次を使用して、現在のsimulationを更新：今、simulationは、各サンプルのサイズについて長さ5のベクトルを返します

## 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)) 

、t統計量のp値のモンテカルロ推定値は、すべての5つの分散推定量のために返されます。それでsapplyに電話すると、マトリックス結果が得られます：

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