6
$\begingroup$

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 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 
$\endgroup$
  • 3
    $\begingroup$ Check out this answer and try repeating with optimizer="nloptwrap" in lmer. If there are multiple optima in the REML criterion than you might not find the true optimum otherwise. Please then edit your post to specify whether or not that fixed your problem. $\endgroup$ – EdM Jan 4 at 18:39
  • $\begingroup$ Data in graph is not good, because we cannot convert it into dataset. How about the variances of estimated intercept and slope? From the data, it seems that the relation between X and Y is very weak, so it is possible that tiny difference in algorithm can result in big difference on point estimates. $\endgroup$ – user158565 Jan 5 at 3:05
  • $\begingroup$ @user158565 Thank you. I edited my post and included the data. The correlation bewteen X and Y is indeed very strong (Pearson r=-0.99). $\endgroup$ – user233068 Jan 5 at 5:13
  • $\begingroup$ @Edm Thanks. I read it and tried a couple things but the results from lme() and lmer() are still completly different. Some of my colleagues used this example to argue that it is important to use a random intercept model (based on the results of lmer). However, I think it is a result due to malfunctioning of lmer or misusage of lmer. In this case, lme showed a correct (?) result. Would you tell me what's excatly going on and whether there is a way to produce similar results between lmer and lme? Thank you. $\endgroup$ – user233068 Jan 5 at 5:41
  • $\begingroup$ It seems lme has problem. I fit the model without random effect and got the exact the same results as lme. stdDev of random intercept is 7.800953e-07, which is wrong based on the estimate of gid effect as fixed. I will compare the results with SAS PROC MIXED. $\endgroup$ – user158565 Jan 5 at 5:50
13
$\begingroup$

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
$\endgroup$
8
$\begingroup$

I agree with @DimitrisRizopoulos's answer, and have a few more points to make.

  • I will start by saying that I am unhappy that lmer doesn't find the best answer - even though I suspect this situation is probably limited to small, unusual (see below) data sets. One of the reasons that lme may do better is that it fits on the log-standard-deviation scale, which may make the minimum near zero "broader".
  • You can get lmer to replicate the lme results by setting an explicit, lower starting value for the scaled standard deviation (start=...); based on the explorations below, start=8 or any lower value should work OK. For what it's worth, this will lead to an estimated random effects variance of 0 (and a "singular fit" message, and an answer that's equivalent to leaving out the random effects component entirely and using lm() ...)
  • In this particular case using the "nloptwrap" optimizer doesn't help; in fact all of the optimizers that lmer can use, starting from the default starting values ($\theta$ (scaled standard deviation) = 1.0), find the higher local minimum away from zero.
  • here is code equivalent to the approach lmer uses to find the starting value by default, when only intercept-valued random effects are present (see here):
v0 <- with(multi,var(ave(y,gid)))  ## variance among group values  
v.e <- var(multi$y)-v0             ## residual var ~ total var - group variance
sqrt(v0/v.e)                       ## convert to scaled standard deviation

This leads to a starting value of $\theta=10.8$.

  • We can see systematically how different starting values give different results:
m0 <- lmer(y~x+(1|(gid)), data=multi, REML=TRUE)
tvec2 <- seq(0,20,length=51)
ff <- function(t0) getME(update(m0,start=t0),"theta")
v <- sapply(tvec2,ff)
plot(tvec2,v)
abline(v=10.8,col="red")

enter image description here

  • We can also explicitly visualize the (negative log-)likelihood surface:
## helper function to capture fitting trajectory
cfun  <- function(...) {
    cc <- capture.output(x <- do.call(lmer,c(list(...),list(verbose=100))))
    gfun <- function(x,s) {
        as.numeric(gsub(s,"",grep(s,x,value=TRUE)))
    }
    it <- gfun(cc,"iteration: +")
    xval <- gfun(cc,"\tx = ")
    fval <- gfun(cc,"\tf\\(x\\) = +")
    attr(x,"optvals") <- data.frame(it,xval,fval)
    return(x)
}

c0 <- cfun(y~x+(1|(gid)), data=multi, REML=TRUE)
c1 <- cfun(y~x+(1|(gid)), data=multi, REML=FALSE)

f <- as.function(m0)
tvec <- seq(0,100,length=101)
dvec <- sapply(tvec,f)
m3 <- update(m0,REML=FALSE)
f2 <- as.function(m3)
dvec2 <- sapply(tvec,f2)
par(las=1,bty="l")
matplot(tvec,cbind(dvec,dvec2),type="l",
        ylab="deviance/REMLcrit",
        xlab="scaled standard dev")
with(attr(c0,"optvals"),text(xval,fval,it))
with(attr(c1,"optvals"),text(xval,fval,it,col=2))
legend("bottomright",c("REML","ML"),
       col=1:2,lty=1:2)

enter image description here

The numbers show the sequence of values tried. We can see that it is only the slightly different shape of the ML curve that tips the optimizer toward the boundary fit rather than the interior fit.

  • Are these data artificial? The left plot below shows the data by group; the right plot shows the values with their group means subtracted. There is almost no variation among the 5 values within each group ...

enter image description here

  • If we simulate data with the same properties (starting from the estimated coefficients), but where the variation is actually Gaussian, we don't get the same kind of multimodal surface at all:
multi_sim <- transform(multi,y=simulate(m0,seed=101)[[1]])
f3 <- as.function(update(m0,data=multi_sim))
dvec3 <- sapply(tvec,f3)
plot(tvec,dvec3,type="l")

enter image description here

$\endgroup$
  • 2
    $\begingroup$ Thank you for your answer. Yes, the data are artificial extracted from a book. The author of the book provided this example to maintain that it is important to use a random intercept model by comparing the results from lm to lmer. Now that I read your comments, lmer with its default option failed to find the best solution. $\endgroup$ – user233068 Jan 8 at 3:26
  • $\begingroup$ maybe contact the author of the book and let them know ... ? (I'm curious what the book is myself, up to you if you feel like saying) $\endgroup$ – Ben Bolker Jan 8 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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