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

[Phillip Chapman]

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;

Thanks,

Phil Chapman