If you want to test whether the change in *R*^{2} is statistically significant for nested linear models with heteroscedasticity-consistent (HC) standard errors (e.g., hierarchical regression), then you can use `vcovHC()`

from the `sandwich`

package and `waldtest()`

from the `lmtest`

package.

Once you have installed and loaded the `sandwich`

and `lmtest`

packages, you can use the following command to compare two linear models with HC standard errors and return an ANOVA:

waldtest(model1, model2, vcov = function(x) {vcovHC(x, type = "HC3")})

Your options for the estimation type (`type`

) are: `type = c("HC3", "const", "HC", "HC0", "HC1", "HC2", "HC4", "HC4m", "HC5").`

`"HC3"`

is the default and `"const"`

is used for the unadjusted standard errors directly from the `lm`

objects.

If you have any issues with using `waldtest()`

, particularly if you’re running it in a function, then you will need to fit the linear models using the `do.call()`

function:

# Fits model1 model1 = do.call("lm", list(formula = y ~ x, data = yourData)) # Fits model2 by adding 'z' term to model1 model2 = do.call("update", list(object = model2, formula = "~ . + z")) # Alternatively, you can enter the full formula for model2 model2 = do.call("lm", list(formula = y ~ x + z, data = yourData))

The above code should look familiar because I used it in the PROCESS Model 1 function described in the previous post.