最近几天我一直在阅读很多关于如何从引导程序中获取回归模型的 p 值（而不是通过排列）。对于模型的每个系数，零假设是系数等于 0，H1 是它不等于 0（双边检验）。
最引人注目的类似主题是stackexchange上的以下两个问题，但答案让我很困惑：
根据我的阅读，我注意到了几种方法，但我不知道哪种方法是有效的：
我的引导函数应该返回为每个样本计算的测试统计量还是估计值？
我应该计算高于 0 还是高于基本模型的点估计的检验统计量/估计的比例？
我应该将结果乘以 2，因为测试是双边的还是使用绝对值？
假设我们要使用 bootstrap 检验回归系数 = 0 的原假设，并假设我们决定 0.05 为显着性水平。现在，我们可以使用 bootstrap 生成每个系数的采样分布。很容易检查 0 是否在 95% 的置信区间内，因此我们可以很容易地决定是否可以拒绝空值。
要获得 p 值，我们需要检查采样分布中 0 的分位数是多少。（我使用的是基于分位数的方法，还有其他方法可以在此处找到Fox on Regression）获得分位数 Q 后，p 值为 2*Q 或 2*(1-Q)，具体取决于Q > 0.5 还是小于 0.5。
作为该方法的说明，请考虑这个
library(faraway)
mdl <- lm(divorce ~ ., data = divusa)
bootTest <- sapply(1:1e4,function(x){
rows <- sample(nrow(divusa),nrow(divusa),replace = T)
mdl <- lm(divorce ~ ., data = divusa[rows,])
return(mdl$coefficients)
})
这是模型
summary(mdl)
Call:
lm(formula = divorce ~ ., data = divusa)
Residuals:
Min 1Q Median 3Q Max
-2.9087 -0.9212 -0.0935 0.7447 3.4689
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 380.14761 99.20371 3.832 0.000274 ***
year -0.20312 0.05333 -3.809 0.000297 ***
unemployed -0.04933 0.05378 -0.917 0.362171
femlab 0.80793 0.11487 7.033 1.09e-09 ***
marriage 0.14977 0.02382 6.287 2.42e-08 ***
birth -0.11695 0.01470 -7.957 2.19e-11 ***
military -0.04276 0.01372 -3.117 0.002652 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.513 on 70 degrees of freedom
Multiple R-squared: 0.9344, Adjusted R-squared: 0.9288
F-statistic: 166.2 on 6 and 70 DF, p-value: < 2.2e-16
注意 p 值
使用 bootstrap 生成的回归系数子集
bootTest[,1:5]
[,1] [,2] [,3] [,4] [,5]
(Intercept) 335.70970574 372.525260160 569.85830341 338.70069977 344.69261238
year -0.18107568 -0.201798080 -0.30579380 -0.18125215 -0.18328105
unemployed -0.02916575 0.006828023 0.01197723 -0.05610887 -0.11463230
femlab 0.79078784 0.842924808 1.02607863 0.77527548 0.76472406
marriage 0.17372382 0.199571033 0.18782967 0.15289119 0.15693996
birth -0.11613752 -0.118507758 -0.11998122 -0.11666450 -0.13344442
military -0.04051730 -0.056277118 -0.04062756 -0.05167556 -0.07251748
使用引导程序生成 p 值
pvals <- sapply(1:nrow(bootTest),function(x) {
distribution <- ecdf(bootTest[x,])
qt0 <- distribution(0)
if(qt0 < 0.5){
return(2*qt0)
} else {
return(2*(1-qt0))
}
})
比较来自 t-test 和 bootstrap 的 p 值
summary(mdl)$coefficients[,4]
(Intercept) year unemployed femlab marriage birth military
2.744830e-04 2.966776e-04 3.621708e-01 1.085196e-09 2.419284e-08 2.191964e-11 2.652003e-03
> pvals
[1] 0.0008 0.0008 0.2196 0.0000 0.0000 0.0000 0.0188
系数 < 1e-8 的高度显着 p 值在 bootstrap 和 1e4 次迭代中均为 0。此外，p 值的排名也具有可比性。