3

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 .

1 Answer 1

3

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

3
  • Hi Ben, thank you very much for your feedback. However, when I take the model presented by you and try to plot the half-normal plot using the hnp function of the hnp package for diagnostic analysis, the following message appears: Erro: number of levels of each grouping factor must be < number of observations (problems: obs). What is the reason for the error?
    – user55546
    Apr 22, 2021 at 13:34
  • 1
    Well, I don't know what the hnp function is doing. Presumably it's trying to re-run the model, and failing to pass the lmerControl specification ... ?
    – Ben Bolker
    Apr 22, 2021 at 13:41
  • Probably yes, in this situation is there any way that I can make this graph taking advantage of the help you provided regarding the adjustment in lmer? I really liked your answer, however, I need to perform the half-normal plot of this model, but for this, the function only works if a mixed effects model is adjusted using lmer, due to this I would like to use your answer .
    – user55546
    Apr 22, 2021 at 13:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.