ANOVA with lm_robust objects?

#1

I’m wondering if there is any way to get lm_robust objects to work with ANOVA. Alternatively, how would I go about trying to achieve the following:

First I run the following:

fit1 <- lm_robust(outcome ~ treatment / moderator - 1, 
                 cluster = teamID,
                 se = "stata",
                 data = data)

fit2 <- lm_robust(outcome ~ treatment - 1, 
                 cluster = teamID,
                 se = "stata",
                 data = data)

Now I want to see if there is an overall effect of the moderator by comparing fit1 with fit2 in an ANOVA.

anova(rescue_events, rescue_events2)

However, I get the error that lm_robust objects don’t work with anova. Is there any way to make them compatible?

#2

This isn’t available out-of-the box as far as I know, but wouldn’t be too hard to implement anova generics for lm_robust. Since it’s 99% compatible with lm, we can try running lm_robust objects through lm’s anova method.


require(estimatr)
m1 <- lm_robust(extra~group, sleep)
m2 <- lm_robust(extra~group+ID, sleep)

debugonce(stats:::anova.lmlist)
stats:::anova.lmlist(m1,m2)

When you run that, it will crash because estimatr doesn’t provide a deviance method. If we hijack also the lm deviance method, it runs, although the output is garbage:

deviance.lm_robust <- stats:::deviance.lm
stats:::anova.lmlist(m1,m2)
Analysis of Variance Table

Model 1: extra ~ group
Model 2: extra ~ group + ID
  Res.Df RSS Df Sum of Sq F Pr(>F)
1     18   0                      
2      9   0  9         0

This is returning zeros, because lm_robust does not provide a residuals method either. We can fill that in:

residuals.lm_robust <- function(object, ...) {
  model.response(model.frame(m1)) - fitted(object)
}

And then we have something matching the output of lm.

> stats:::anova.lmlist(m1,m2)
Analysis of Variance Table

Model 1: extra ~ group
Model 2: extra ~ group + ID
  Res.Df    RSS Df Sum of Sq       F    Pr(>F)   
1     18 64.886                                  
2      9  6.808  9    58.078 8.53085 0.0019014 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#3

Also note that using base R’s anova method will decompose RSS to create it’s test statistics; this will basically ignore the robust standard errors specified on the model.

#4

Thanks so much Neal!

Two quick questions:

I am running the following:

residuals.lm_robust <- function(object, ...) {
  model.response(model.frame(m1)) - fitted(object)
}

debugonce(stats:::anova.lmlist)
stats:::anova.lmlist(fit1,fit2)

And getting the error:
Error in stats:::anova.lmlist(fit1, fit2) :
(list) object cannot be coerced to type ‘double’

Any thoughts?

Also, what do you mean in your second message? Do you mean the cluster-robust standard errors I specified in lm_robust won’t be reflected in the resulting ANOVA output?

#5

ANOVA tests relies strongly on the assumption of variance homogeneity, and the test statistics are basically likelihood ratio tests. That approach doesn’t work if you also want to allow for heteroskedasticity of un-specified form, where (as far as I know) you’d have to use a Wald test or score test. For one-way ANOVA, the onewaytests package provides a bunch of different variations on a Wald test. For two-way ANOVA, the Wald_test() function from clubSandwich package can be used to test constraints on Model 1.

#6

Typo in the residuals function, change “m1” to “object”. Also make sure you did

deviance.lm_robust <- stats:::deviance.lm

Correct, the above literally does the exact same thing lm would by literally reusing lm’s code.

@James_Pustejovsky’s answer is solid, I would also add that you can also get (some of) those tests out of car::linearHypothesis as a third option. There’s some functionality for that built in to estimatr, which I haven’t really used personally, but should integrate more directly with DeclareDesign if you are using that for simulation. see ?estimatr::lh_robust.

#7

Building off @nfultz, would this work for you?


library(estimatr)
library(lmtest)
m1 <- lm_robust(extra~group, sleep)
m2 <- lm_robust(extra~group+ID, sleep)
waldtest(m1, m2)
#8

Thanks so much, @Alex_Coppock. That’s really helpful