# different linear mixed model coefficients in lme4 vs nlme

i'm trying to compare the coefficients for the same linear mixed model in lme4 vs nlme, see this example using the penguins dataset.

I can't work out why they are different? why is the intercept the same across the 3 groups when using nlme?

library(tidyverse)
library(palmerpenguins)
library(lme4)
library(nlme)

db <- penguins %>%
filter(!is.na(flipper_length_mm), !is.na(bill_length_mm), !is.na(body_mass_g))

lme4_fit <- lme4::lmer(
body_mass_g ~ flipper_length_mm + bill_length_mm + (1+flipper_length_mm|species),
REML = TRUE,
data = db
)

nlme_fit <- nlme::lme(
body_mass_g ~ flipper_length_mm + bill_length_mm,
random = ~ 1+flipper_length_mm|species,
method = "REML",
data = db
)

coef(lme4_fit)
coef(nlme_fit)


• For me, the lme4 model does not converge. Having only three groups as a random intercept is generally going to cause problems in a model like this
– sjp
Commented Aug 8, 2022 at 17:24
• i have a feeling it's more me not extracting the coefficients correctly. When i extract the preds from both models they are very close. Commented Aug 8, 2022 at 17:41

Okay, so after looking at the models a little bit more closely, the issue is that the model fit using nlme has a singular fit. The estimate of the SD for the random intercepts is tiny, so that's why almost the same number comes back for each group in the data when you call coef(nlme_fit).

See this post and the comments: nlme estimates near zero variance for the random effects

In this specific example, if we call summary() instead of coef() and look at all parts of the model reported there. That's where I saw that the Intercept in the random effects was 322.54 in the model fit with lme4, whereas in the one fit with nlme it was 4.042139e-04. I've commented the lines where this information is located in each model summary.

Let us know if you have any other questions.

summary(lme4_fit)

Linear mixed model fit by REML ['lmerMod']
Formula: body_mass_g ~ flipper_length_mm + bill_length_mm + (1 + flipper_length_mm |      species)
Data: db

REML criterion at convergence: 4946.3

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.46124 -0.66304 -0.09222  0.61054  3.12864

Random effects:
Groups   Name              Variance  Std.Dev. Corr
species  (Intercept)       104030.92 322.54            # ESTIMATE LME4 HERE
flipper_length_mm     14.59   3.82   -0.98
Residual                   115095.78 339.26
Number of obs: 342, groups:  species, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)       -3943.198    577.935  -6.823
flipper_length_mm    26.749      3.835   6.975
bill_length_mm       60.460      7.097   8.520

Correlation of Fixed Effects:
(Intr) flpp__
flppr_lngt_ -0.850
bll_lngth_m -0.018 -0.401
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge: degenerate  Hessian with 1 negative eigenvalues

summary(nlme_fit)

Linear mixed-effects model fit by REML
Data: db
AIC    BIC    logLik
4960.718 4987.5 -2473.359

Random effects:
Formula: ~1 + flipper_length_mm | species
Structure: General positive-definite, Log-Cholesky parametrization
StdDev       Corr
(Intercept)       4.042139e-04 (Intr)  # ESTIMATE NLME HERE
flipper_length_mm 2.308314e+00 0.941
Residual          3.394369e+02

Fixed effects:  body_mass_g ~ flipper_length_mm + bill_length_mm
Value Std.Error  DF   t-value p-value
(Intercept)       -4025.568  550.7142 337 -7.309723       0
flipper_length_mm    27.117    3.4150 337  7.940672       0
bill_length_mm       60.773    7.0958 337  8.564632       0
Correlation:
(Intr) flpp__
flipper_length_mm -0.794
bill_length_mm    -0.022 -0.448

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-2.43161474 -0.67327147 -0.08989753  0.62206892  3.11662877

Number of Observations: 342
Number of Groups: 3