effects package with multiway clustered variance-covariance matrix

Question

I'd like to plot an interaction effect with an error ribbon using the "effects" package. However, I want to calculate errors based on a multi-way clustered covariance matrix as can be calculated with the "multiwayvcov" package. "effects" can take a function to calculate a covariance matrix (vcov.=) but it looks like the function won't accept any additional arguments.

Example:

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")

The above produces a similar plot for the interactions using non-clustered se. I want something that does the same but with clustered se. Something like (imagine that "gear" and "cyl" are cluster ids in the mtcars data):

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))))

Any help is much appreciated

It looks like model should be replaced with m1 when you create tmp. — eipi10
That was an error in the example code. I fixed it but still have the same problem. — Richard Benton
What actually happens when you run the second block of code? When I run it I get Error in Effect.lm(predictors, mod, vcov. = vcov., ...) : could not find function "vcov.". — eipi10
I get the same error. What seems to be happening is that vcov argument in the effect() function won't accept additional arguments from the function you feed it to estimate the covariance matrix. In the documentation, vcov will accept the hccm function from the car package to calculate heteroscedasticity-corrected covariance matrices: 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)))) But if you try to feed hccm additional arguments you get the same error. — Richard Benton

Richard Benton Richard Benton · Accepted Answer · 2015-12-16T16:54:19

I cooked up a little work around that might be useful for others in the future. It looks like the vcov argument in the effect function will only accept functions with one argument and that has to be the lm (or whatever) object from the first part of effect. If you write an external function that adjusts the other arguments for the covariance matrix function internally it will work.

examples from above: Standard Covariance Matrix

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 covariance matrix (from car package) with an option:

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))))

And finally, multiway clustered covariance matrix (again imagining that gear and cyl are cluster ids in the mtcars data):

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))))

The last example throws you a NaNs error but that is because of the particulars of the example.

effects package with multiway clustered variance-covariance matrix

2 Answers