2

This question asks the same question, but hasn't been answered. My question relates to how to specify the model with the lm() function and is therefore a programming (not statistical) question.

I have a mixed design (2 repeated and 1 independent predictors). Participants were first primed into group A or B (this is the independent predictor) and then they rated how much they liked 4 different statements (these are the two repeated predictors). There are many great online resources how to model this data. However, my data is heterscedastic. So I like to use heteroscedastic-consistent covariance matrices. This paper explains it well. The sandwich and lmtest packages are great. Here is a good explanation how to do it for a indpendent design in R with lm(y ~ x).

It seems that I have use lm, else it wont work?

Here is the code for a regression model assuming that all variances are equal (which they are not as Levene's test comes back significant). 

fit3 <- nlme:::lme(DV ~ repeatedIV1*repeatedIV2*independentIV1, random =  ~1|participants, df) ##works fine

Here is the code for an indepedent model correcting for heteroscedasticity, which works.

fit3 <- lm(DV ~ independentIV1)
library(sandwich)
vcovHC(fit3, type = 'HC4', sandwich = F)
library(lmtest)
coef(fit3, vcov. = vcovHC, type = 'HC4')

So my question really is, how to specify my model with lm? Alternative approaches in R how to fit my model accounting for heteroscedasticity are welcome too!

Thanks a lot!!!

Simone
  • 497
  • 5
  • 19
  • 2
    `lm` can't fit mixed models. You'll need to use the `nlme` or `lme4` package or another package that can fit mixed models. [**This tutorial on mixed models in R**](https://idaejin.github.io/bcam-courses/neiker-2016/material/mixed-models/#heteroskedasticity) might be helpful. The link will take you right to the heteroscedasticity section. – eipi10 Jan 27 '18 at 20:59
  • 2
    Also check out the CRAN task views on [Robust Statistical Methods](https://cran.r-project.org/web/views/Robust.html) and [Econometrics](https://cran.r-project.org/web/views/Econometrics.html). – eipi10 Jan 27 '18 at 21:05
  • @eipi10 The weighted least squares regression (tutorial on mixed models in R) requires me to correctly identify the pattern of heteroscedasticity in my residuals. I am not sure I can do that with my current knowledge and as usual there is a deadline limiting the amount of time I can invest in reading up on this. That’s the beauty of the HC approach – it doesn’t require to specify assumptions. – Simone Jan 29 '18 at 09:56
  • @eipi10 The first link gives a couple of alternative options on how to use HC covariance matrices (robustbase, robus, MASS package). Will have a look at them. Thanks!! The second link mainly talks about HC and HAC based on the sandwich and lmtest packages, which only work with lm() output, but not with nlme. Achim has e-mailed back saying that merDeriv works with the package lme4. Have asked him to post the information here. Again thanks a lot for sharing the links. Really useful! – Simone Jan 29 '18 at 09:59

1 Answers1

2

My impression is that your problems come from mixing various approaches for various aspects (repeated measurements/correlation vs. heteroscedasticity) that cannot be mixed so easily. Instead of using random effects you might also consider fixed effects, or instead of only adjusting the inference for heteroscedasticity you might consider a Gaussian model and model both mean and variance, etc. For me, it's hard to say what is the best route forward here. Hence, I only comment on some aspects regarding the sandwich package:

The sandwich package is not limited to lm/glm only but it is in principle object-oriented, see vignette("sandwich-OOP", package = "sandwich") (also published as doi:10.18637/jss.v016.i09.

There are suitable methods for a wide variety of packages/models but not for nlme or lme4. The reason is that it's not so obvious for which mixed-effects models the usual sandwich trick actually works. (Disclaimer: But I'm no expert in mixed-effects modeling.)

However, for lme4 there is a relatively new package called merDeriv (https://CRAN.R-project.org/package=merDeriv) that supplies estfun and bread methods so that sandwich covariances can be computed for lmer output etc. There is also a working paper associated with that package: https://arxiv.org/abs/1612.04911

Achim Zeileis
  • 15,710
  • 1
  • 39
  • 49