具备多路聚类方差协方差矩阵的影响包。

Question

我想使用“effects”包绘制带有误差线的交互作用效应。但是，我想基于多路聚类协方差矩阵计算误差，就像使用“multiwayvcov”包所能计算的那样。“effects”可以使用函数来计算协方差矩阵（vcov.=），但似乎该函数不接受任何其他参数。

示例：

library(effects)
library(ggplot2)
library(multiwayvcov)
m1 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
tmp <- as.data.frame(effect("hp:wt", model, vcov=vcov, se=TRUE, xlevels=list(wt=c(2.2,3.2,4.2))))
ggplot(data=tmp, aes(x=hp, y=fit, colour=as.factor(wt))) +
   geom_line() +
   geom_ribbon(aes(ymin=lower, ymax=upper, fill=as.factor(wt)), alpha = 0.5) +
   labs(colour="wt")

上面的代码使用非聚类标准误差生成了一个类似的交互图。我想要做同样的事情，但是使用聚类标准误差。可以尝试如下代码（假设"gear"和"cyl"是mtcars数据集中的聚类id）：

m1 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
cluster.ids <- data.frame(i = mtcars$gear, j = mtcars$cyl)
tmp <- as.data.frame(effect("hp:wt", m1, vcov=cluster.vcov(m1, cluster = cluster.ids), xlevels=list(wt=c(2.2,3.2,4.2))))

任何帮助都非常感激。

- Richard Benton

---

在创建 tmp 时，看起来应该用 m1 替换 model。 - eipi10

那个示例代码中有一个错误。我已经修复了，但仍然存在同样的问题。 - Richard Benton

当你运行第二块代码时，实际上会发生什么？ 当我运行它时，我得到了“Effect.lm(predictors, mod, vcov. = vcov., ...)中的错误：无法找到函数"vcov."”。 - eipi10

我得到了相同的错误。看起来发生的情况是effect()函数中的vcov参数不会接受你提供的用于估计协方差矩阵的函数的其他参数。在文档中，vcov将接受car包中的hccm函数来计算异方差-校正协方差矩阵：library(car) m2 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars) tmp2 <- as.data.frame(effect("hp:wt", m2, vcov=hccm, xlevels=list(wt=c(2.2,3.2,4.2))))但是如果您尝试为hccm提供其他参数，则会收到相同的错误。 - Richard Benton

Answer 1

我构思了一个小技巧，可能对将来的其他人有用。看起来effect函数中的vcov参数只接受带有一个参数的函数，而且必须是effect的第一部分lm（或其他）对象。如果编写一个外部函数，在内部调整协方差矩阵函数的其他参数，它就会起作用。

以上是例子： 标准协方差矩阵

m1 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
tmp1 <- as.data.frame(effect("hp:wt", m1, vcov=vcov, se=TRUE,  xlevels=list(wt=c(2.2,3.2,4.2))))

ggplot(data=tmp, aes(x=hp, y=fit, colour=as.factor(wt))) +
geom_line() +
geom_ribbon(aes(ymin=lower, ymax=upper, fill=as.factor(wt)), alpha = 0.5) +
labs(colour="wt")

library(car)
m2 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
hccmfunc <- function(x) {  
return(hccm(x, type="hc0"))}
tmp2 <- as.data.frame(effect("hp:wt", m2, vcov=hccmfunc,   xlevels=list(wt=c(2.2,3.2,4.2))))

最后，我们来看多路聚类协方差矩阵（假设在mtcars数据中gear和cyl是聚类ID）：

m3 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
cluster.ids <- data.frame(i = mtcars$gear, j = mtcars$cyl)
mwfunc <- function(x) {return(cluster.vcov(x, cluster= cluster.ids))}
tmp <- as.data.frame(effect("hp:wt", m3, vcov=mwfunc, xlevels=list(wt=c(2.2,3.2,4.2))))

最后一个例子会抛出NaNs错误，但这是由于该例子的特殊情况。

- Richard Benton

Answer 2

根据2021年6月的更新和影响包版本4.2-0： 基础代码略有更改，因此需要按照以下要求进行源代码更改：

trace(effects:::Effect.lm, edit = T) # this is not permanent but has to be redone every time RStudio is loaded
# change line 162 in the windows that opens from
V <- vcov.(mod, complete = FALSE) 
# to
change to V <- vcov.(mod)
# click save

在这个更改之后，您可以使用Richard Benton建议的方法，并进行轻微调整：现在效果公式中使用vcov. [带点]而不是仅使用vcov。

mwfunc <- function(x) {return(cluster.vcov(x, cluster= cluster.ids))}
tmp <- as.data.frame(effect("hp:wt", m3, vcov.=mwfunc, xlevels=list(wt=c(2.2,3.2,4.2))))

- JNWHH