Как создать двойную петлю?

Следующий код работает довольно хорошо, НО: я должен менять размер выборки n = 25, 50,... и оценку дисперсии каждый раз перед запуском кода. Я хотел бы решить эту проблему с помощью цикла.

Далее я кратко опишу код. В коде создается 1000 регрессионных моделей для заданного размера выборки n. Затем каждая регрессионная модель из 1000 оценивается с помощью МНК. После этого я вычисляю статистику t на основе различных значений бета x3 из 1000 выборок. Нулевая гипотеза гласит: 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

r loops economics
person targa    schedule 18.11.2016    source источник
comment
Вам не нужна петля. Создайте список n<-list(25, 50, 100, 250, 500, 1000). Затем используйте lapply или, если у вас есть два списка, mapply. Это выведет список результатов для каждого размера выборки n.   -  person Morgan Ball    schedule 18.11.2016

Ответы (1)


arrow_upward
0
arrow_downward

Напишите функцию, которая запускает симуляцию

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
person ExperimenteR    schedule 18.11.2016
comment
Большое спасибо! К сожалению, было бы лучше, если бы я получил результаты решений различных оценщиков дисперсии. Таким образом, результатом будет матрица. В лучшем случае было бы, если бы на выходе для каждой переменной (x1, x2, x3) было бы три матрицы. Вы тоже умеете это делать? - person targa; 18.11.2016