ループを2倍にするにはどうすればいいですか？

Question

次のコードはうまくいきますが、コードを実行する前に毎回サンプルサイズn = 25、50、...と分散推定を変更する必要があります。 私はループでこの問題を解決したいと思います。

以下、コードについて簡単に説明します。コード内で、与えられたサンプルサイズnの1000個の回帰モデルが作成されます。次に、1000のうちの各回帰モデルがOLSによって推定される。その後、1000サンプルのうちx3の異なるベータ値に基づいてt統計量を計算します。この仮説は、H0：beta03 = beta3となり、x3の計算されたベータ値は1と定義された「真の」値に等しい。最後のステップでは、虚偽仮説がどのくらいの頻度で棄却されるかを調べる（有意水準= 0.05）。 私の最終的な目標は、各サンプルサイズと分散推定値についてのヌル仮説の賛成拒絶率を吐き出すコードを作成することです。あなたの誰かが私にそれを手伝ってもらえれば嬉しいです。

#sample size n = 25, 50, 100, 250, 500, 1000 
n <- 50 
B <- 1000 

#'real' beta values 
beta0 <- 1 
beta1 <- 1 
beta2 <- 1 
beta3 <- 1 

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


} 

alpha <- 0.05 
test.decision <- abs(t.test.values) < qt(p=c(1-alpha/2), df=n-4) 
length(test.decision[test.decision==FALSE])/B 

Answer 1

は、シミュレーションに

library(car) 
sample_size = c("n=25"=25, "n=50"=50, "n=100"=100, "n=250"=250, "n=500"=500, "n=1000"=1000) 
B <- 100 
beta0 <- 1 
beta1 <- 1 
beta2 <- 1 
beta3 <- 1 
alpha <- 0.05 

sim <- 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) 
    exi <- 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, sim, beta3h0 = 0.7) 
# n=25 n=50 n=100 n=250 n=500 n=1000 
# 0.92 0.88 0.92 0.79 0.44 0.24