CODEWITHSUNDEEP
rstatisticsregressionlme4mixed-models
user55546
user55546

Reputation: 62

Equivalence of a mixed model fitted by lme and lmer

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

Upvotes: 1

Views: 669

Answers (1)

Ben Bolker
Ben Bolker

Reputation: 226247

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)

Upvotes: 2

Related Questions