模拟数千个回归并获得 p 值
Simulating thousands of regressions and obtaining p-values
我希望在 R 中进行一些基本模拟,以检查 p 值的性质。我的目标是查看大样本量是否趋向于小 p 值。我的想法是生成 1,000,000 个数据点的随机向量,将它们相互回归,然后绘制 p 值的分布并寻找偏差。
到目前为止我是这么想的:
x1 = runif(1000000, 0, 1000)
x2 = runif(1000000, 0, 1000)
model1 = lm(x2~x1)
使用从另一个线程获取的代码:
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
lmp(model1)
0.3874139
关于我如何为 1000 个或更多模型执行此操作的任何建议?谢谢!
参见 ?replicate
...但是您正在计算的 p 值假设高斯误差不是均匀误差(并不是说这在 n=10^6 时很重要)
具体来说,是这样的:
nrep <- 1000
ndat <- 1000000
results <- replicate(nrep, {
x1=runif(ndat, 0, 1000);
x2=runif(ndat, 0, 1000);
model1=lm(x1 ~ x2);
lmp(model1)
})
应该可以,但是
运行.
需要很长时间
我建议缩小 nrep 和 ndat 来尝试一下。
我希望在 R 中进行一些基本模拟,以检查 p 值的性质。我的目标是查看大样本量是否趋向于小 p 值。我的想法是生成 1,000,000 个数据点的随机向量,将它们相互回归,然后绘制 p 值的分布并寻找偏差。
到目前为止我是这么想的:
x1 = runif(1000000, 0, 1000)
x2 = runif(1000000, 0, 1000)
model1 = lm(x2~x1)
使用从另一个线程获取的代码:
lmp <- function (modelobject) {
if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
f <- summary(modelobject)$fstatistic
p <- pf(f[1],f[2],f[3],lower.tail=F)
attributes(p) <- NULL
return(p)
}
lmp(model1)
0.3874139
关于我如何为 1000 个或更多模型执行此操作的任何建议?谢谢!
参见 ?replicate
...但是您正在计算的 p 值假设高斯误差不是均匀误差(并不是说这在 n=10^6 时很重要)
具体来说,是这样的:
nrep <- 1000
ndat <- 1000000
results <- replicate(nrep, {
x1=runif(ndat, 0, 1000);
x2=runif(ndat, 0, 1000);
model1=lm(x1 ~ x2);
lmp(model1)
})
应该可以,但是 运行.
需要很长时间我建议缩小 nrep 和 ndat 来尝试一下。