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

Question

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

data(mtcars)
library(multcomp)
lm1 <- lm(mpg ~ cyl + hp, data = mtcars)
summary(glht(lm1, linfct = 'cyl + hp = 0'))

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

但是，假设我想在第三个变量上对标准错误进行聚类：

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

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

Error in modelparm.default(model, ...) : no ‘coef’ method for ‘model’ found!

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

谢谢！

Answer 1

glht(ct1, linfct = 'cyl + hp = 0') 不会起作用，因为 ct1 不是 glht 对象，不能通过 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)

如果我们在 LinearCombTest 中不指定 .vcov，它与 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升级为Get p-value for group mean difference without refitting linear model with a new reference level，这里我们可以测试任意组合，组合系数为alpha:

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

而不只是总和

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

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

How to conduct linear hypothesis test on regression coefficients with a clustered covariance matrix?

r

regression

linear-regression

covariance

lm