루프를 두 배로 만드는 방법은 무엇입니까?

Question

다음 코드는 효과적이지만, 코드를 실행하기 전에 매번 샘플 크기 n = 25, 50, ... 및 분산 추정기를 변경해야합니다. 루프를 사용하여이 문제를 해결하고 싶습니다.

다음은 코드를 간단히 설명합니다. 코드 내에서 주어진 샘플 크기 n에 대한 1000 개의 회귀 모델이 생성됩니다. 그런 다음, 1000 개 중 각 회귀 모델은 OLS에 의해 추정됩니다. 그 후 1000 개의 샘플 중 x3의 다른 베타 값을 기반으로 t 통계를 계산합니다. null hypothesisation은 H0 : beta03 = beta3, 즉 x3의 계산 된 베타 값은 1로 정의 된 '실제'값과 같습니다. 마지막 단계에서는 null 가설이 거부되는 빈도 (유의 수준 = 0.05)를 확인합니다. 최종 목표는 각 표본 크기 및 분산 추정치에 대한 무효 가설의 procentual rejection rate을 뱉어내는 코드를 만드는 것입니다. 여러분 중 누구라도 저를 도울 수 있다면 기뻐할 것입니다.

#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