有没有办法提取 coeftest() 函数的 $R^2$?

Is there a way to extract $R^2$ of the coeftest() function?

我有一个关于函数 coeftest() 的问题。我想弄清楚从哪里可以得到这个函数的 R 平方的任何结果。我正在拟合一个标准的多元线性回归,如下所示:

Wetterstation.lm <- lm(temp~t+temp_auto+dum.jan+dum.feb+dum.mar+dum.apr+dum.may+dum.jun+dum.aug+dum.sep+dum.oct+dum.nov+dum.dec+
                      dum.jan*t+dum.feb*t+dum.mar*t+dum.apr*t+dum.may*t+dum.jun*t+dum.aug*t+dum.sep*t+dum.oct*t+dum.nov*t+dum.dec*t)

预先我分别定义了每个变量,我的结果如下:


> summary(Wetterstation.lm)

Call:
lm(formula = temp ~ t + temp_auto + dum.jan + dum.feb + dum.mar + 
    dum.apr + dum.may + dum.jun + dum.aug + dum.sep + dum.oct + 
    dum.nov + dum.dec + dum.jan * t + dum.feb * t + dum.mar * 
    t + dum.apr * t + dum.may * t + dum.jun * t + dum.aug * t + 
    dum.sep * t + dum.oct * t + dum.nov * t + dum.dec * t)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.9564  -1.3214   0.0731   1.3621   9.9312 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.236e+00  9.597e-02  33.714  < 2e-16 ***
t            1.206e-05  3.744e-06   3.221  0.00128 ** 
temp_auto    8.333e-01  2.929e-03 284.503  < 2e-16 ***
dum.jan     -3.550e+00  1.252e-01 -28.360  < 2e-16 ***
dum.feb     -3.191e+00  1.258e-01 -25.367  < 2e-16 ***
dum.mar     -2.374e+00  1.181e-01 -20.105  < 2e-16 ***
dum.apr     -1.582e+00  1.142e-01 -13.851  < 2e-16 ***
dum.may     -7.528e-01  1.106e-01  -6.809 9.99e-12 ***
dum.jun     -3.283e-01  1.106e-01  -2.968  0.00300 ** 
dum.aug     -2.144e-01  1.094e-01  -1.960  0.05005 .  
dum.sep     -8.009e-01  1.103e-01  -7.260 3.96e-13 ***
dum.oct     -1.752e+00  1.123e-01 -15.596  < 2e-16 ***
dum.nov     -2.622e+00  1.181e-01 -22.198  < 2e-16 ***
dum.dec     -3.287e+00  1.226e-01 -26.808  < 2e-16 ***
t:dum.jan    2.626e-06  5.277e-06   0.498  0.61877    
t:dum.feb    2.479e-06  5.404e-06   0.459  0.64642    
t:dum.mar    1.671e-06  5.277e-06   0.317  0.75145    
t:dum.apr    1.357e-06  5.320e-06   0.255  0.79872    
t:dum.may   -3.173e-06  5.276e-06  -0.601  0.54756    
t:dum.jun    2.481e-06  5.320e-06   0.466  0.64098    
t:dum.aug    5.998e-07  5.298e-06   0.113  0.90985    
t:dum.sep   -5.997e-06  5.321e-06  -1.127  0.25968    
t:dum.oct   -5.860e-06  5.277e-06  -1.110  0.26683    
t:dum.nov   -4.215e-06  5.320e-06  -0.792  0.42820    
t:dum.dec   -2.526e-06  5.277e-06  -0.479  0.63217    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.12 on 35744 degrees of freedom
Multiple R-squared:  0.9348,    Adjusted R-squared:  0.9348 
F-statistic: 2.136e+04 on 24 and 35744 DF,  p-value: < 2.2e-16

现在我尝试使用函数 coeftest() 和 vcovHAC 调整异方差和自相关,如下所示:

library("lmtest")
library("sandwich")
Wetterstation.lm.HAC <- coeftest(Wetterstation.lm, vcov = vcovHAC)

这些结果是:

> Wetterstation.lm.HAC

t test of coefficients:

               Estimate  Std. Error  t value  Pr(>|t|)    
(Intercept)  3.2356e+00  7.8816e-02  41.0529 < 2.2e-16 ***
t            1.2059e-05  2.7864e-06   4.3280 1.509e-05 ***
temp_auto    8.3334e-01  2.9798e-03 279.6659 < 2.2e-16 ***
dum.jan     -3.5505e+00  1.1843e-01 -29.9789 < 2.2e-16 ***
dum.feb     -3.1909e+00  1.2296e-01 -25.9507 < 2.2e-16 ***
dum.mar     -2.3741e+00  1.0890e-01 -21.8002 < 2.2e-16 ***
dum.apr     -1.5821e+00  9.5952e-02 -16.4881 < 2.2e-16 ***
dum.may     -7.5282e-01  8.8987e-02  -8.4599 < 2.2e-16 ***
dum.jun     -3.2826e-01  8.2271e-02  -3.9899 6.622e-05 ***
dum.aug     -2.1440e-01  7.7966e-02  -2.7499  0.005964 ** 
dum.sep     -8.0094e-01  8.4456e-02  -9.4835 < 2.2e-16 ***
dum.oct     -1.7519e+00  9.2919e-02 -18.8538 < 2.2e-16 ***
dum.nov     -2.6224e+00  1.0028e-01 -26.1504 < 2.2e-16 ***
dum.dec     -3.2873e+00  1.1393e-01 -28.8546 < 2.2e-16 ***
t:dum.jan    2.6256e-06  5.2429e-06   0.5008  0.616517    
t:dum.feb    2.4790e-06  5.5284e-06   0.4484  0.653850    
t:dum.mar    1.6713e-06  4.8632e-06   0.3437  0.731107    
t:dum.apr    1.3567e-06  4.5670e-06   0.2971  0.766423    
t:dum.may   -3.1734e-06  4.2970e-06  -0.7385  0.460209    
t:dum.jun    2.4809e-06  4.1490e-06   0.5979  0.549880    
t:dum.aug    5.9983e-07  4.0379e-06   0.1485  0.881910    
t:dum.sep   -5.9975e-06  4.1675e-06  -1.4391  0.150125    
t:dum.oct   -5.8595e-06  4.4635e-06  -1.3128  0.189265    
t:dum.nov   -4.2151e-06  4.6555e-06  -0.9054  0.365263    
t:dum.dec   -2.5257e-06  4.9871e-06  -0.5065  0.612539    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

但是因为我想在 table 中添加 R 平方来总结我的结果,所以我不知道如何得到它。现在我想知道是否有人可以帮助解决这个问题并告诉我可以从哪里获得信息。也许我只是愚蠢,但由于我对 R 很陌生,我很乐意得到任何帮助。

非常感谢您。

简单的回答:不,没有。而且也没有理由这样做。 coeftest() 函数正在使用给定模型的值。对于 stats4::coefcoeftest 函数采用模型的系数。

如果函数打算这样做,则可以提取 r^2 值。另外imtestcoeftest()只是returns一个贴,所以不能提取值。

总结一下:lmtest::coeftest() 使用的是 lm 的值,因此 r^2 没有改变


关于标准误差的背景:lm 使用稍微不同的方法来计算标准误差。在source code你可以提取:

  R <- chol2inv(Qr$qr[p1, p1, drop = FALSE])
  se <- sqrt(diag(R))

所以这给我们带来了以下内容:lm 使用 Cholesky 组合(我也认为 R 默认使用它)。

coeftest()同时使用sqrt()残差的方差-协方差矩阵(见here)。所以自相关。 vcov 提取给定模型的方差-协方差矩阵(如lm

se <- vcov.
se <- sqrt(diag(se))

我个人认为 lm 的用户出于速度原因正在使用 chalesky 组合,因为您不必反转矩阵。但是 lmtest 包的作者并没有意识到这一点。但这只是猜测。

因为 t 值是在两个包中使用估计值和标准误差计算的,如下所示:

tval <- as.vector(est)/se

并且 p-value 基于 tval:

2*pt(abs(tval), rdf, lower.tail = FALSE)

所有差异都是基于不同的估计误差。请注意,估计是相同的,因为 coeftest() 只是继承了它们。