I've been playing around with both `nlme::lme`

and `lme4::lmer`

.

I fitted a simple random intercepts model using `lme()`

and `lmer()`

. As you can see below, I got completely different results from `lmer()`

and `lme()`

. Even signs of coefficients are different! Am I doing something wrong? I also fitted an empty model with the two packages. In this case, the results were practically the same (results not shown). Would you educate me to understand this issue? Unless I made a mistake, I think there is something wrong with the `lme4`

package.

` multi <- structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95, 4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12, 2.07, 1.05, 1, 1.07, 1.16, 1.11), y = c(3.2, 3.21, 3.256, 3.25, 3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628, 3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954, 4, 3.994, 4), pid = 1:25, gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L)), class = "data.frame", row.names = c(NA, -25L)) #lme > lme(y~x, random=~1|gid,data=multi,method="REML") Linear mixed-effects model fit by REML Data: multi Log-restricted-likelihood: 41.76745 Fixed: y ~ x (Intercept) x 4.1846756 -0.1928357 #lmer lmer(y~x+(1|(gid)), data=multi, REML=T) Linear mixed model fit by REML ['lmerMod'] Formula: y ~ x + (1 | (gid)) Data: multi REML criterion at convergence: -78.4862 Random effects: Groups Name Std.Dev. (gid) (Intercept) 0.70325 Residual 0.02031 Number of obs: 25, groups: (gid), 5 Fixed Effects: (Intercept) x 2.8152 0.2638 `

**Contents**hide

#### Best Answer

As it was noted in this answer, and also mentioned in one of the comments, the problem seems to be a local maximum. To see this more clearly I have written below a simple code to calculate the negative log-likelihood of this model and do the optimization using `optim()`

. Starting with different initial values leads to the two different solutions:

`# data multi <- structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95, 4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12, 2.07, 1.05, 1, 1.07, 1.16, 1.11), y = c(3.2, 3.21, 3.256, 3.25, 3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628, 3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954, 4, 3.994, 4), pid = 1:25, gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L)), class = "data.frame", row.names = c(NA, -25L)) # function to calculate the negative log-likelihood of the random intercepts model library("mvtnorm") logLik <- function (thetas, y, X, id) { ncX <- ncol(X) betas <- thetas[seq_len(ncX)] sigma_b <- exp(thetas[ncX + 1]) sigma <- exp(thetas[ncX + 2]) eta <- c(X %*% betas) unq_id <- unique(id) n <- length(unq_id) lL <- numeric(n) for (i in seq_len(n)) { id_i <- id == unq_id[i] n_i <- sum(id_i) V_i <- matrix(sigma_b^2, n_i, n_i) diag(V_i) <- diag(V_i) + sigma^2 lL[i] <- dmvnorm(y[id_i], mean = eta[id_i], sigma = V_i, log = TRUE) } - sum(lL, na.rm = TRUE) } # optimization using as initial values 0 for the fixed effects, # and 1 for the variance components opt <- optim(rep(0, 4), logLik, method = "BFGS", y = multi$y, X = cbind(1, multi$x), id = multi$gid) opt$par[1:2] # fixed effects #> [1] 2.855872 0.250341 exp(opt$par[3]) # sd random intercepts #> [1] 0.6029724 exp(opt$par[4]) # sd error terms #> [1] 0.01997889 # optimization using as initial values 4 & -0.2 for the fixed effects, # and 0.0003 and 0.034 for the variance components opt2 <- optim(c(4, -0.2, -8, -3.4), logLik, method = "BFGS", y = multi$y, X = cbind(1, multi$x), id = multi$gid) opt2$par[1:2] # fixed effects #> [1] 4.1846965 -0.1928397 exp(opt2$par[3]) # sd random intercepts #> [1] 0.000270746 exp(opt2$par[4]) # sd error terms #> [1] 0.03239167 `

### Similar Posts:

- Solved – R lme4 interpret MLM confidence intervals
- Solved – Can a variable be included in a mixed model as a fixed effect and as a random effect at the same time
- Solved – Proportion of variance in dependent variable accounted for by predictors in a mixed effects model
- Solved – Proportion of variance in dependent variable accounted for by predictors in a mixed effects model
- Solved – REML or ML to compare two mixed effects models with differing fixed effects, but with the same random effect