Search code examples
rstatisticslme4mixed-modelsnlme

Heterocesdastic model of mixed effects via lmer function


I am adjusting a mixed effects model which, due to the observed heteroscedasticity, it was necessary to include an effect to accommodate it. Therefore, using the lme function of the nlme package, this was easy to be solved, see the code below:

library(nlme)
library(lme4)
Model1 <- lme(log(Var1)~log(Var2)+log(Var3)+
                     (Var4)+(Var5),
                    random = ~1|Var6, Data1, method="REML",
                   weights = varIdent(form=~1|Var7))
#Var6: It is a factor with several levels.
#Var7: It is a Dummy variable.

However, I need to readjust the model described above using the lme4 package, that is, using the lmer function. It is known and many are the materials that inform some limitations existing in the lme4, such as, for example, modeling heteroscedasticity. What motivated me to readjust this model is the fact that I have an interest in using a specific package that in cases of mixed models it only accepts if they are adjusted through the lmer function. How could I resolve this situation? Below is a good part of the model adjusted using the lmer function, however, this model is not considering the effect to model the observed heteroscedasticity.

Model2 <- lmer(log(Var1)~log(Var2)+log(Var3)+
                          (Var4)+(Var5)+(1|Var6),
                    Data1, REML=T)

Regarding the choice of the random effect (Var6) and the inclusion of the effect to consider the heterogeneity by levels of the variable (Var7), these were carefully analyzed, however, I will not put here the whole procedure so as not to be an extensive post and to be more objective .


Solution

  • This is hackable. You need to add an observation-level random effect that is only applied to the group with the larger residual variance (you need to know this in advance!), via (0+dummy(Var7,"1")|obs); this has the effect of multiplying each observation-level random effect value by 1 if the observation is in group "1" of Var7, 0 otherwise. You also need to use lmerControl() to override a few checks that lmer does to try to make sure you are not adding redundant random effects.

    Data1$obs <- factor(seq(nrow(Data1)))
    Model2 <- lmer(log(Var1)~log(Var2)+log(Var3)+
                       (Var4)+(Var5) + (1|Var6) +
                       (0+dummy(Var7,"1")|obs),
                   Data1, REML=TRUE,
                   control=lmerControl(check.nobs.vs.nlev="ignore",
                                       check.nobs.vs.nRE="ignore"))
    
    all.equal(REMLcrit(Model2), c(-2*logLik(Model1))) ## TRUE
    all.equal(fixef(Model1), fixef(Model2), tolerance=1e-7)
    

    If you want to use this model with hnp you need to work around the fact that hnp doesn't pass the lmerControl option properly.

    library(hnp)
    d <- function(obj) resid(obj, type="pearson")
    s <- function(n, obj) simulate(obj)[[1]]
    f <- function(y.) refit(Model2, y.)
    
    hnp(Model2, newclass=TRUE, diagfun=d, simfun=s, fitfun=f)
    

    You might also be interested in the DHARMa package, which does similar simulation-based diagnostics.

    half-normal plot