如何使用聚类协方差矩阵对回归系数进行线性假设检验？

Question

我有兴趣计算 R 中线性回归后系数线性组合的估计值和标准误差。例如，假设我有回归和测试：

\n\n

data(mtcars)\nlibrary(multcomp)\nlm1 <- lm(mpg ~ cyl + hp, data = mtcars)\nsummary(glht(lm1, linfct = 'cyl + hp = 0'))\n

\n\n

cyl这将估计和上的系数之和的值hp，并根据 生成的协方差矩阵提供标准误差lm。

\n\n

但是，假设我想将我的标准错误聚集在第三个变量上：

\n\n

data(mtcars)\nlibrary(multcomp)\nlibrary(lmtest)\nlibrary(multiwayvcov)\nlm1 <- lm(mpg ~ cyl + hp, data = mtcars)\nvcv <- cluster.vcov(lm1, cluster = mtcars$am)\nct1 <- coeftest(lm1,vcov. = vcv)\n

\n\n

ct1包含我的聚类的 SE am。但是，如果我尝试使用ct1中的对象glht，则会收到错误消息

\n\n

\n

modelparm.default(model, ...) 中的错误：\n 未找到 \xe2\x80\x98model\xe2\x80\x99 的 \xe2\x80\x98coef\xe2\x80\x99 方法！

\n

\n\n

关于如何使用聚类方差协方差矩阵进行线性假设有什么建议吗？

\n\n

谢谢！

\n

Answer 1

glht(ct1, linfct = 'cyl + hp = 0')不起作用，因为ct1不是一个glht对象，不能被强制为这样的 via as.glht。我不知道是否有一个包或现有的函数可以做到这一点，但是我们自己解决这个问题并不困难。下面的小函数可以做到这一点：

LinearCombTest <- function (lmObject, vars, .vcov = NULL) {
  ## if `.vcov` missing, use the one returned by `lm`
  if (is.null(.vcov)) .vcov <- vcov(lmObject)
  ## estimated coefficients
  beta <- coef(lmObject)
  ## sum of `vars`
  sumvars <- sum(beta[vars])
  ## get standard errors for sum of `vars`
  se <- sum(.vcov[vars, vars]) ^ 0.5
  ## perform t-test on `sumvars`
  tscore <- sumvars / se
  pvalue <- 2 * pt(abs(tscore), lmObject$df.residual, lower.tail = FALSE)
  ## return a matrix
  matrix(c(sumvars, se, tscore, pvalue), nrow = 1L,
         dimnames = list(paste0(paste0(vars, collapse = " + "), " = 0"),
                         c("Estimate", "Std. Error", "t value", "Pr(>|t|)")))
  }

我们来测试一下：

data(mtcars)
lm1 <- lm(mpg ~ cyl + hp, data = mtcars)
library(multiwayvcov)
vcv <- cluster.vcov(lm1, cluster = mtcars$am)

如果我们.vcov在 中未指定LinearCombTest，则与 相同multcomp::glht：

LinearCombTest(lm1, c("cyl","hp"))

#              Estimate Std. Error   t value     Pr(>|t|)
#cyl + hp = 0 -2.283815  0.5634632 -4.053175 0.0003462092

library(multcomp)
summary(glht(lm1, linfct = 'cyl + hp = 0'))

#Linear Hypotheses:
#              Estimate Std. Error t value Pr(>|t|)    
#cyl + hp == 0  -2.2838     0.5635  -4.053 0.000346 ***

如果我们提供协方差，它就会满足您的要求：

LinearCombTest(lm1, c("cyl","hp"), vcv)

#              Estimate Std. Error  t value    Pr(>|t|)
#cyl + hp = 0 -2.283815  0.7594086 -3.00736 0.005399071

评论

LinearCombTest在获取组均值差异的 p 值处升级，无需使用新的参考水平重新拟合线性模型，我们可以在其中测试具有组合系数的任何组合alpha：

alpha[1] * vars[1] + alpha[2] * vars[2] + ... + alpha[k] * vars[k]

而不仅仅是总和

vars[1] + vars[2] + ... + vars[k]