2

I have fitted a mixed effects model considering both functions widely used in R, namely: the lme function from the nlme package and the lmer function from the lme4 package. To readjust the model from lme to lme4, following the same reparametrization, I used the following information from this topic, being that is only possible to do this in lme4 in a hackable way.: Heterocesdastic model of mixed effects via lmer function

I apologize for hosting the data in a link, however, I couldn't find an internal R database that has variables that might match my problem.

Data: https://drive.google.com/file/d/1jKFhs4MGaVxh-OPErvLDfMNmQBouywoM/view?usp=sharing

The fitted models were:

library(nlme)
library(lme4)
ModLME = lme(Var1~I(Var2)+I(Var2^2),
                         random = ~1|Var3,
                         weights = varIdent(form=~1|Var4),
                         Dataone, method="REML")

ModLMER = lmer(Var1~I(Var2)+I(Var2^2)+(1|Var3)+(0+dummy(Var4,"1")|Var5),
                Dataone, REML = TRUE,
                control=lmerControl(check.nobs.vs.nlev="ignore",
                check.nobs.vs.nRE="ignore"))

Which are equivalent, see:

all.equal(REMLcrit(ModLMER), c(-2*logLik(ModLME))) 
[1] TRUE
all.equal(fixef(ModLME), fixef(ModLMER), tolerance=1e-7)
[1] TRUE

> summary(ModLME)

Linear mixed-effects model fit by REML
  Data: Dataone 
        AIC       BIC   logLik
  -209.1431 -193.6948 110.5715

Random effects:
 Formula: ~1 | Var3
        (Intercept)   Residual
StdDev:  0.05789852 0.03636468

Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | Var4 
 Parameter estimates:
       0        1 
1.000000 5.641709 
Fixed effects:  Var1 ~ I(Var2) + I(Var2^2) 
                 Value  Std.Error DF  t-value p-value
(Intercept)  0.9538547 0.01699642 97 56.12093       0
I(Var2)     -0.5009804 0.09336479 97 -5.36584       0
I(Var2^2)   -0.4280151 0.10038257 97 -4.26384       0

summary(ModLMER)
Linear mixed model fit by REML. t-tests use Satterthwaites method [lmerModLmerTest]
Formula: Var1 ~ I(Var2) + I(Var2^2) + (1 | Var3) + (0 + dummy(Var4, "1") |  
    Var5)
   Data: Dataone
Control: lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE = "ignore")

REML criterion at convergence: -221.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.1151 -0.5891  0.0374  0.5229  2.1880 

Random effects:
 Groups   Name             Variance  Std.Dev. 
 Var3     (Intercept)      6.466e-12 2.543e-06
 Var5     dummy(Var4, "1") 4.077e-02 2.019e-01
 Residual                  4.675e-03 6.837e-02
Number of obs: 100, groups:  Var3, 100; Var5, 100

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)  0.95385    0.01700 95.02863  56.121  < 2e-16 ***
I(Var2)     -0.50098    0.09336 92.94048  -5.366 5.88e-07 ***
I(Var2^2)   -0.42802    0.10038 91.64017  -4.264 4.88e-05 ***

However, when observing the residuals of these models, note that they are not similar. See that in the model adjusted by lmer, mysteriously appears a residue with the shape of a few points close to a straight line. So, how could you solve such a problem so that they are identical? I believe the problem is in the lme4 model.

aa=plot(ModLME, main="LME")
bb=plot(ModLMER, main="LMER")
gridExtra::grid.arrange(aa,bb,ncol=2)

enter image description here

user55546
  • 37
  • 1
  • 15

1 Answers1

3

I can tell you what's going on and what should in principle fix it, but at the moment the fix doesn't work ...

The residuals being plotted take all of the random effects into account, which in the case of the lmer fit includes the individual-level random effects (the (0+dummy(Var4,"1")|Var5) term), which leads to weird residuals for the Var4==1 group. To illustrate this:

plot(ModLMER, col = Dataone$Var4+1)

enter image description here

i.e., you can see that the weird residuals are exactly the ones in red == those for which Var4==1.

In theory we should be able to get the same residuals via:

res <- Dataone$Var1 - predict(ModLMER, re.form = ~(1|Var3))

i.e., ignore the group-specific observation-level random effect term. However, it looks like there is a bug at the moment ("contrasts can be applied only to factors with 2 or more levels").

An extremely hacky solution is to construct the random-effect predictions without the observation-level term yourself:

## fixed-effect predictions
p0 <- predict(ModLMER, re.form = NA)
## construct RE prediction, Var3 term only:
Z <- getME(ModLMER, "Z")
b <- drop(getME(ModLMER, "b"))
## zero out observation-level components
b[101:200] <- 0
## add RE predictions to fixed predictions
p1 <- drop(p0 + Z %*% b)
## plot fitted vs residual
plot(p1, Dataone$Var1 - p1)

For what it's worth, this also works:

library(glmmTMB)
ModGLMMTMB <- glmmTMB(Var1~I(Var2)+I(Var2^2)+(1|Var3),
                      dispformula = ~factor(Var4),
                      REML = TRUE,
                      data = Dataone)
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Hello Ben! Would there be a solution with which I could still use the ``lmer`` function? The first option given by the subtraction seems interesting to me, however, in fact it is giving a bug involving the error already mentioned by you. The second option I am facing error, see: ``Error in .Call("getParameterOrder", data, parameters, new.env(), PACKAGE = DLL) :`` – user55546 Jul 09 '22 at 23:59
  • I'm pretty sure the `glmmTMB` solution will work with the development version (`remotes::install_github("glmmTMB/glmmTMB/glmmTMB")`). I will try to fix the `lme4` bug sometime soon (https://github.com/lme4/lme4/issues/691) but can't promise a timeline. What do you need to do with the residuals ... ? – Ben Bolker Jul 10 '22 at 00:03
  • Hello Ben! I'm doing a simulation study with the model fitted in ``lme4``. The problem is that as it is carrying out these "atypical" residuals, in my simulation the error is carried to it. That is, I have a well-tuned model in ``lme`` but not in ``lmer`` due to the problem highlighted by you. – user55546 Jul 10 '22 at 00:13
  • Although you presented the "extremely hacky" solution, the residuals of the model doing it manually do not match that of the model fitted by the ``lme`` function. I need to include the heterogeneity of variances, since it is through this that the additional effects are captured and then I get a satisfactory model. – user55546 Jul 10 '22 at 00:41
  • I understand that you want to fit a model with heteroscedasticity in the residuals. One point is that `lme`'s default plot is based on *Pearson* residuals, so we would also have to divide by the residual SE (or sqrt(resid var + extra var) in the Var4==1 group). The other problem is that this is a very poorly constrained (unidentifiable?) - try `intervals(ModLME)` - which will easily lead to different estimates for the residual variance and hence for the Pearson residuals – Ben Bolker Jul 10 '22 at 16:55
  • Right. I understand, but, taking into account the model fitted by ``lmer``, there is no other way (that is bug free) that I can plot the residuals and they will be similar to the residuals of the model fitted in ``lme``? – user55546 Jul 10 '22 at 19:51
  • It's going to be very difficult because the Pearson residuals are scaled according to the residual SD. If the residual SD estimates are highly unstable, then it's going to be difficult to to get stable Pearson residuals. Fixing the bug in `lme4` won't fix any of this - it will just get you to the same point that the hack (i.e., adding the random-effect component for `Var3` to the fixed-effect predictions manually) will. (If you want to discuss this further it might be worth posting to `r-sig-mixed-models@r-project.org`, as Stack Overflow is explicitly not a discussion forum ...) – Ben Bolker Jul 10 '22 at 20:33