将 monte carlo p 值排列到不同样本大小和方差估计量的矩阵中

Question

以下代码运行良好（基于）。但是每次运行代码之前，我都必须更改方差估计器 (ols、hc0、hc1、hc2、hc3) .我想用一个循环来解决这个问题。

下面简单介绍一下代码。在代码中，为每个样本大小 (n = 25, 50, 100, 250, 500, 1000) 创建了 1000 个回归模型。然后，1000 个回归模型中的每个回归模型都由 OLS 进行估计。之后，我根据 1000 个样本中 x3 的不同 beta 值计算 t 统计量。原假设为：H0: beta03 = beta3，即 x3 的计算 beta 值等于我定义为 1 的 'real' 值。在最后一步中，我检查了原假设出现的频率拒绝（显着性水平 = 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：

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

现在，simulation 将 return 一个长度为 5 的向量。对于每个样本大小，monte carlo t 统计 p 值的估计是 returned对于所有 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

将 monte carlo p 值排列到不同样本大小和方差估计量的矩阵中

Arrange monte carlo p-value into a matrix for different sample size and variance estimators

r

loops

regression

linear-regression

p-value