如何复现sandwich包NeweyWest函数的P值？求相关文档说明

一、复现P值的步骤 #

你用coeftest(l, vcov.=nw_vcov)得到的P值，是基于t检验的双侧P值，计算逻辑可拆解为以下几步：

1. 计算t统计量：用线性模型的回归系数，除以Newey-West方法得到的标准误

   t_stat <- coef(l) / nw_se
   

2. 确定t分布的自由度：等于模型的残差自由度，即样本量n减去模型参数总数（含截距项）。你的示例中n=100，参数共11个（截距+10个自变量），因此自由度为100 - 11 = 89
3. 计算双侧P值：利用t分布累积分布函数pt()，公式为2 * (1 - pt(abs(t统计量), 自由度))

完整复现代码 #

将上述步骤整合到你的代码中，可验证复现结果与coeftest输出完全一致：

library(sandwich)
library(lmtest)

# 你的示例数据集代码
set.seed(123) 
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
x4 <- rnorm(n)
x5 <- rnorm(n)
x6 <- rnorm(n)
x7 <- rnorm(n)
x8 <- rnorm(n)
x9 <- rnorm(n)
x10 <- rnorm(n)

y <- 2 + 1.5*x1 - 0.5*x2 + 0.3*x3 + 2*x4 + x5 - x6 + 0.8*x7 - 1.2*x8 + 0.5*x9 + 
  1.1*x10 + rnorm(n)

d <- data.frame(x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, y)
eq <- y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10

# 拟合模型并计算Newey-West协方差矩阵
l <- lm(eq, d) 
nw_vcov <- NeweyWest(l, lag=1, prewhite=FALSE, adjust=TRUE)
nw_se <- sqrt(diag(nw_vcov))

# 复现P值
df_resid <- df.residual(l)  # 直接调用模型内置的残差自由度，更通用
t_stat <- coef(l) / nw_se
p_values <- 2 * (1 - pt(abs(t_stat), df = df_resid))

# 和coeftest的结果对比
l_nw_se <- coeftest(l, vcov.=nw_vcov) 
s <- as.data.frame(l_nw_se[, ])

# 验证一致性（结果完全相同）
cbind(复现的P值 = p_values, coeftest输出的P值 = s$Pr...t..)

二、相关说明文档 #

1. R包内置帮助文档
   • 输入?NeweyWest可查看sandwich包中Newey-West协方差矩阵的计算细节，包括滞后阶数、adjust=TRUE调整项的处理逻辑
   • 输入?coeftest可查看lmtest包中系数检验的方法，明确了P值是基于t分布（线性模型默认）计算的双侧检验结果
2. sandwich包官方手册
   运行vignette("sandwich")可查看包的详细使用手册，其中包含异方差自相关一致（HAC）协方差估计（含Newey-West）的理论与实现说明
3. 原始文献
   Newey-West方法的核心逻辑来自Newey, W. K., & West, K. D. (1987). A Simple, Positive Semi-Definite, Heteroskedasticity and Autocorrelation Consistent Covariance Matrix. Econometrica, 55(3), 703-708.

内容的提问来源于stack exchange，提问作者J.K.