[R-sig-ME] Convergence problem with lmer with non-centered WWheat data

[Douglas Bates]

On Thu, Dec 3, 2009 at 4:18 PM, Phillip Chapman
<pchapman at stat.colostate.edu> wrote:
> Dear All,
>
> I am trying to run a random coefficients regression using the winter wheat
> data from the "SAS System for Mixed Models" book. lmer appears to be
> converging to an erroneous solution for the second model (diagonal cov
> matrix for random effects).   Is there a better way to specify my lmer model
> (h1 in the text below)?  Note: this convergence problem goes away when I
> center (i.e. subtract the mean of) the continuous predictor, Moisture.  I
> would usually center predictors in random coefficients models anyway.
>
> The data frame is WWheat from the "SASmixed" package.  The response is
> Yield; the continuous predictor is Moisture; and the factor with 10 levels
> is Variety.
>
> I have no problem with the first model in the book; it fits random
> intercept/slope pairs for each variety.  The following three analyses
> produce almost same results, out to a reasonable number of significant
> digits:
>
> 1.   Using lme4:       g1 <- lmer(Yield ~ Moisture + (Moisture | Variety),
> data=WWheat)
>
> 2.   Using nlme:        g2 <- lme(Yield ~ Moisture, random = ~ Moisture |
> Variety, data=WWheat)
>
> 3.  Using SAS:
>  proc mixed data=wheat scoring=8;
>  class variety;
>  model yield = moist/solution;
>  random int moist/sub=variety type=un solution G Gcorr;
>  run;
>
> The second model restricts the covariance matrix of the intercept/slope
> pairs to be diagonal.  For this model I get the same result using lme and
> SAS, but lmer converges to a different solution, which appears to be
> erroneous:
>
> 1.  Using lme4:       h1 <- lmer(Yield ~ Moisture + (1 | Variety) +
> (Moisture -1 | Variety), data=WWheat)
> 2. Using nlme (I could only figure out how to do this after creating a
> grouped data object):
>                     WWheat2 <- groupedData(Yield ~ Moisture | Variety,
> data=WWheat)
>                      h2 <- lme(Yield ~ Moisture, random =  pdDiag(~
> Moisture), data=WWheat2)
>
> 3. Using SAS:
>  proc mixed data=wheat scoring=8;
>  class variety;
>  model yield = moist/solution;
>  random int moist/sub=variety type=un(1) solution G Gcorr;
>  run;
>
> I took the variance estimates from lmer and put them into SAS Proc Mixed to
> see if the REML deviances match: they do. Both SAS and lmer give the same
> deviance of 205.7, which is much greater than the optimum value of 187.1
> produced by SAS and lme.  The fixed effects estimates from SAS and lmer also
> match when these variance estimates are used.
> proc mixed data=wheat scoring=8;
>  class variety;
>  model yield = moist/solution;
>  random int moist/sub=variety type=un(1) solution G Gcorr;
>  parms (1.6854e+01)(0) ( 2.6133e-16) (7.5295e-01)/noiter; *the solution from
> lmer;
>  run;

It does look like a convergence problem because of the 2.6133e-16 for
the variance of the random effects for slope.  Could you tell us what
version of the lme4 package you are using?  It is easiest just to send
the output of sessionInfo()

This problem will occur because of the need to use a constrained
nonlinear optimizer to determine the ML or REML estimates.  The
evaluation of the deviance is done very carefully so as to give the
optimizer a good shot at the getting the optimum but it doesn't always
work.

In the development branch of the lme4 package (the lme4a branch) there
is the capability of switching to a different optimizer.  I have found
the bobyqa optimizer in the package by Kate Mullen and John Nash to be
more robust and faster in most cases than is nlminb, which is the
default. However, even nlminb converges to the correct optimum in that
branch.