如何创建双循环? [英] How to create doubled a loop?
问题描述
以下代码效果很好,但是:在运行代码之前,我必须每次更改样本大小n = 25、50,...和方差估计量. 我想用循环来解决这个问题.
The following code works out quite well, BUT: I have to change the sample size n = 25, 50, ... and the variance estimator everytime before I run the code. I would like to solve this problem with a loop.
此后,我简要介绍一下代码.在代码中,为给定样本大小n创建了1000个回归模型.然后,通过OLS估算1000个样本中的每个回归模型.之后,我根据1000个样本中x3的不同beta值计算t统计量.零假设假设为:H0:beta03 = beta3,即x3的计算得出的beta值等于我定义为1的真实"值.在最后一步,我检查零假设假设被拒绝的频率(显着性水平= 0.05). 我的最终目标是创建一个代码,为每个样本量和方差估计量吐出零假设的百分比拒绝率.如果你们中的任何一个可以帮助我,我将很高兴.在这里您可以看到我的代码:
Hereafter, I briefly describe the code. Within the code, 1000 regression models for a given sample size n 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 x3 out of the 1000 samples. The nullhypothessis reads: H0: beta03 = beta3, that is the calculated beta value of x3 equals the 'real' value which I defined as 1. In the last step, I check how often the nullhypothesis is rejected (significance level = 0.05). My final goal is to create a code which spits out the procentual rejection rate of the nullhypothesis for each sample size and variance estimator. I would be pleased if anyone of you could help me with that. Here you can see my code:
#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
推荐答案
编写一个运行模拟的函数
Write a function that runs a simulation
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))
}
并使用lapply
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
这篇关于如何创建双循环?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!