lme() and lmer() giving conflicting results

Question

I have been working with some data that has some problems with repeated measurements. In doing so I noticed very different behavior between lme() and lmer() using my test data and want to know why.

The fake data set I created has height and weight measurements for 10 subjects, taken twice each. I set up the data so that between subjects there would be a positive relationship between height and weight, but a negative relationship between the repeated measures within each individual.

set.seed(21)
Height=1:10; Height=Height+runif(10,min=0,max=3) #First height measurement
Weight=1:10; Weight=Weight+runif(10,min=0,max=3) #First weight measurement

Height2=Height+runif(10,min=0,max=1) #second height measurement
Weight2=Weight-runif(10,min=0,max=1) #second weight measurement

Height=c(Height,Height2) #combine height and wight measurements
Weight=c(Weight,Weight2)

DF=data.frame(Height,Weight) #generate data frame
DF$ID=as.factor(rep(1:10,2)) #add subject ID
DF$Number=as.factor(c(rep(1,10),rep(2,10))) #differentiate between first and second measurement

Here is a plot of the data, with lines connecting the two measurements from each individual.

So I ran two models, one with lme() from the nlme package and one with lmer() from lme4. In both cases I ran a regression of weight against height with a random effect of ID to control for the repeated measurements of each individual.

library(nlme)
Mlme=lme(Height~Weight,random=~1|ID,data=DF)
library(lme4)
Mlmer=lmer(Height~Weight+(1|ID),data=DF)

These two models often (though not always depending on the seed) generated completely different results. I have seen where they generate slightly different variance estimates, calculate different degrees of freedom, etc., but here the coefficients are in opposite directions.

coef(Mlme)
#   (Intercept)    Weight
#1   1.57102183 0.7477639
#2  -0.08765784 0.7477639
#3   3.33128509 0.7477639
#4   1.09639883 0.7477639
#5   4.08969282 0.7477639
#6   4.48649982 0.7477639
#7   1.37824171 0.7477639
#8   2.54690995 0.7477639
#9   4.43051687 0.7477639
#10  4.04812243 0.7477639

coef(Mlmer)
#   (Intercept)    Weight
#1     4.689264 -0.516824
#2     5.427231 -0.516824
#3     6.943274 -0.516824
#4     7.832617 -0.516824
#5    10.656164 -0.516824
#6    12.256954 -0.516824
#7    11.963619 -0.516824
#8    13.304242 -0.516824
#9    17.637284 -0.516824
#10   18.883624 -0.516824

To illustrate visually, model with lme()

And model with lmer()

Why are these models diverging so much?

Answer 1

tl;dr if you change the optimizer to "nloptwrap" I think it will avoid these issues (probably).

Congratulations, you've found one of the simplest examples of multiple optima in a statistical estimation problem! The parameter that lme4 uses internally (thus convenient for illustration) is the scaled standard deviation of the random effects, i.e. the among-group std dev divided by the residual std dev.

Extract these values for the original lme and lmer fits:

(sd1 <- sqrt(getVarCov(Mlme)[[1]])/sigma(Mlme))
## 2.332469
(sd2 <- getME(Mlmer,"theta")) ## 14.48926

Refit with another optimizer (this will probably be the default in the next release of lme4):

Mlmer2 <- update(Mlmer,
  control=lmerControl(optimizer="nloptwrap"))
sd3 <- getME(Mlmer2,"theta")   ## 2.33247

Matches lme ... let's see what's going on. The deviance function (-2*log likelihood), or in this case the analogous REML-criterion function, for LMMs with a single random effect takes only one argument, because the fixed-effect parameters are profiled out; they can be computed automatically for a given value of the RE standard deviation.

ff <- as.function(Mlmer)
tvec <- seq(0,20,length=101)
Lvec <- sapply(tvec,ff)
png("CV38425.png")
par(bty="l",las=1)
plot(tvec,Lvec,type="l",
     ylab="REML criterion",
     xlab="scaled random effects standard deviation")
abline(v=1,lty=2)
points(sd1,ff(sd1),pch=16,col=1)
points(sd2,ff(sd2),pch=16,col=2)
points(sd3,ff(sd3),pch=1,col=4)
dev.off()

I continued to obsess further over this and ran the fits for random seeds from 1 to 1000, fitting lme, lmer, and lmer+nloptwrap for each case. Here are the numbers out of 1000 where a given method gets answers that are at least 0.001 deviance units worse than another ...

          lme.dev lmer.dev lmer2.dev
lme.dev         0       64        61
lmer.dev      369        0       326
lmer2.dev      43        3         0

In other words, (1) there is no method that always works best; (2) lmer with the default optimizer is worst (fails about 1/3 of the time); (3) lmer with "nloptwrap" is best (worse than lme 4% of the time, rarely worse than lmer).

To be a little bit reassuring, I think that this situation is likely to be worst for small, misspecified cases (i.e. residual error here is uniform rather than Normal). It would be interesting to explore this more systematically though ...