我想复制
NeweyWest
包中 sandwich
函数的 p 值。如果有任何解释文档可以告诉我吗?
# Question: R: Replicating P-value of `NeweyWest` of `sandwich` package.
library(sandwich)
library(lmtest)
# Sample dataset with 10 predictor variables
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
l <- lm(eq, d)
# Create Newey-West covariance matrix estimates
nw_vcov <- NeweyWest(l, lag=1, prewhite=FALSE, adjust=TRUE)
# compute the square root of the diagonal elements in vcov
nw_se <- sqrt(diag(nw_vcov))
# Apply estimates to linear model
l_nw_se <- coeftest(l, vcov.=nw_vcov)
# Extract the results
s <- as.data.frame(l_nw_se[, ])
# Question: how do I replicate the p-value?
# Could you let me know if there is any documentation for the explanation?
P 值通常使用学生 t 分布
pt()
以及 t 统计数据和自由度来计算。对于通常的双尾测试,也适用于lmtest::coeftest
,我们可以这样做:
> (p_values <- 2*pt(q=-abs(s$`t value`), df=l$df.residual))
[1] 9.286402e-38 9.345499e-23 1.359084e-05 9.134816e-02 2.707676e-40
[6] 2.393224e-15 7.685513e-12 5.457149e-13 1.957496e-15 6.150043e-04
[11] 5.182621e-21
> stopifnot(identical(p_values, s$`Pr(>|t|)`))