在R中替换回归模型中的标准误差

Question

我正在寻找一种方法，可以直接用我的标准误差替换回归模型中的标准误差，以便在另一个R包中使用鲁棒模型。该R包没有自己的鲁棒选项，只能提供特定类型的模型而不是coeftest格式。

假设我有一个线性模型：

model <- lm(data=dat, Y ~ X1 + X2 + X3)

我希望您能为我提供强健的标准误差结果：

robust <- coeftest(model, vcov=sandwich)

接下来，我需要在一个特定的包中使用这个模型，但该包无法提供coeftest，并且没有自己的鲁棒标准误选项。我想在将其馈入该包之前，替换原始模型的标准误差（以及p值、t统计量等），以便进行考虑。

要访问原始模型中的标准误差，我使用：

summary(model)$coefficients[,2]

为了从coeftest中提取标准误差，我使用以下方法：

coeftest.se <- robust[, 2]

然而，下面的方法在尝试替换模型的标准误差时返回错误，因为它将“summary”本身视为一个命令。

summary(model)$coefficients[,2] <- coeftest.se


Error in summary(M3)$coefficients[, 2] <- seM3 : could not find function "summary<-"

具体内容

我正在尝试使用 Mediation R 包运行中介分析。该包将使用“mediate”函数进行一种聚类标准误差，但我希望使用双向聚类标准误差。

为了获得双向聚类标准误差，我正在使用 Mahmood Arai 的 mclx 函数（代码可以在这里找到（第4页））。

我的想法是向该包的 mediate 函数提供已经报告正确标准误差的模型。根据文档，中介分析包接受以下模型类别：lm、polr、bayespolr、glm、bayesglm、gam、rq、survreg 和 merMod。

- firebird17139

---

因此，summary(model)是模型的一个函数，而不是模型的一部分。标准误差不是模型的一部分，摘要函数会即时计算它们（请参见summary.lm以了解详细信息）。 - Gregor Thomas

因此，最简单的解决方案可能是修改您的“特定包”以使用强健的标准误差，或者它可以被馈送到一个rlm对象中。无论哪种方式，在不知道您要在哪里使用它之前，都无法提供太多帮助。 - Gregor Thomas

这真的取决于您将对模型进行的函数。lm对象是一个“活动”的对象。它知道用来拟合它的数据和原始调用。修改其内部属性并不保证会产生您想要的行为。您确实应该提供一个更完整的可复现示例，以展示您希望在模型下游做什么。 - MrFlick

我已经包含了一些具体的细节。我主要是从理论角度考虑这个问题，但很快我会提供一个可重现的例子。 - firebird17139

Answer 1

我的建议是：

对我有用的方法是：

    model <- lm(data=dat, Y ~ X1 + X2 + X3)

    library(sandwich)

    SE_robust <- sqrt(diag(vcovHC(model, type="HC2")))

    model2 <- summary(model)

    model2$coefficients[,2] <- SE_robust

- My Seppo

Answer 2

@MySeppo的回答很好，但我认为在此提供一个函数，使所有内容（SE、t-统计量、p-值）更加健壮以增强通用性（我也使用HC1以与Stata保持一致）是值得的：

model <- lm(mpg~cyl,data=mtcars)

robustsummary <- function(model) { 
  library(sandwich)
  library(lmtest)
  coeftest <- coeftest(model, vcov = vcovHC(model, type="HC1"))
  summ <- summary(model)
  summ$coefficients[,2] <- coeftest[,2]
  summ$coefficients[,3] <- coeftest[,3]
  summ$coefficients[,4] <- coeftest[,4]
  summ
}
summary(model)
robustsummary(model)

或者如果您不想使用coeftest

robustsummary <- function(model) { 
  library(sandwich)
  SE_robust <- sqrt(diag(vcovHC(model, type="HC1")))
  summ <- summary(model)
  summ$coefficients[,2] <- SE_robust
  summ$coefficients[,3] <- summ$coefficients[,1]/SE_robust # get t-stats
  summ$coefficients[,4] <- 2*(1-pt(abs(summ$coefficients[,3]),summ$df[2])) # apply t distribution with the right degrees of freedom to get p-values
  summ
}

- Richard DiSalvo

---

嗨@Richard DiSalvo，你的获取强大摘要的解决方案非常完美！谢谢。不过，有没有可能在“模型”内部调整数值？这样当我运行summary(model)时，我能得到一个强大的摘要？我问这个问题是因为我想进一步输入模型来预测数值和计算逆米尔斯比。你能给我指点一下吗？非常感谢。 - undefined

@fsure 我不认为这是可能的，因为我认为summary()或coeftest()使用模型来计算标准误差。鲁棒标准误差只是使用了与summary()默认使用的不同公式。然而，我不认为鲁棒标准误差会影响模型的预测值。但它会影响条件均值和预测区间的置信区间。对于逆米尔斯比率，我不太确定。我建议您查看您将模型输入的函数，并查看它们是否有一种方法来使用鲁棒（或更重要的是，聚类）标准误差公式。 - undefined