具有多路聚类方差 – 协方差矩阵的特效包

2023年8月14日 282次阅读

我想使用“效果”包绘制与错误功能区的交互效果.但是,我想基于多路聚类协方差矩阵计算误差,可以使用“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")

上面为使用非聚集se的交互产生了类似的图.我想要的东西是相同的,但与群集se.类似的东西(想象“齿轮”和“圆柱”是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))))

任何帮助深表感谢

最佳答案 我做了一些可能对其他人有用的工作.看起来效果函数中的vcov参数只接受带有一个参数的函数,并且必须是效果第一部分的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")

hccm协方差矩阵(来自汽车包装)带有一个选项：

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数据中的簇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错误,但这是因为示例的细节.

点赞