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

Answer 1

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.

Answer 2

Update as per June 2021 and effects package version 4.2-0: The underlying code has slightly changed so that a source code change is required as per the following:

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

After this change you can use the approach suggested by Richard Benton with a slight adjustment: it is now vcov. [with a dot] instead of only vcov in the effects formula.

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