[R-sig-ME] Apparent false convergence in lmer with some BIBD data
Phillip Chapman
pchapman at stat.colostate.edu
Thu Oct 29 20:05:22 CET 2009
Greetings ME Board members,
Is it possible to manipulate the convergence parameters in lmer?
The control options on the lmer help page are very limited. I would like
to alter the convergence criteria. Also, it would be nice to be able to
input my own starting values. Can someone tell me whether this is
possible and
point me to an appropriate reference.
I am using the "Pillow" BIBD data from page 1066 of Ott and
Longnecker, 5th edition. My lmer program seems to have a false
convergence with
the random block estimate at zero. At the risk of committing heresy, I
have posted
SAS analysis of the same data. The SAS solution has a lower
REMLDeviance than the lmer solution.
I also ran the SAS code with the lmer solution to confirm that the
problem is with the
optimization, not with the data or the REML likelihood. lmer and SAS
give the same
REMLDeviance at the lmer solution. I have plotted the likelihood
surface (vertical=block
variance, horizontal=error variance) and see that there is a ridge
running from SW to NE on
plot. The lmer solution is not on that ridge.
Any help or suggestions would be appreciated.
Thanks,
Phil Chapman
Here is the data:
blk pillow firmness
1 A 59
1 B 26
1 C 38
2 D 85
2 E 92
2 F 69
3 G 74
3 H 52
3 I 27
4 A 62
4 D 70
4 G 68
5 B 27
5 E 98
5 H 59
6 C 31
6 F 60
6 I 35
7 A 63
7 E 85
7 I 30
8 B 22
8 F 73
8 G 75
9 C 45
9 D 74
9 H 51
10 A 52
10 F 76
10 H 43
11 B 18
11 D 79
11 I 41
12 C 41
12 E 84
12 G 81
Here is my program and some of its output.
data <- read.table(file="c:/temp/Pillow data.txt",header=T)
attach(data)
blk <- factor(blk)
install.packages("lme4")
library(lme4)
h <- lmer(firmness ~ pillow + (1|blk), REML=TRUE, verbose=TRUE)
print(h, digits=10)
Linear mixed model fit by REML
Formula: firmness ~ pillow + (1 | blk)
AIC BIC logLik deviance REMLdev
207.3991367 224.8178450 -92.69956833 220.2067619 185.3991367
Random effects:
Groups Name Variance Std.Dev.
blk (Intercept) 2.6584e-10 1.6305e-05
Residual 3.5398e+01 5.9496e+00
Number of obs: 36, groups: blk, 12
Here is a SAS program with the same model:
proc mixed data=pillow;
class blk pillow;
model firmness=pillow;
random blk;
The Mixed Procedure
Iteration History
Iteration Evaluations -2 Res Log Like Criterion
0 1 185.39913672
1 2 185.35185521 0.00000016
2 1 185.35184420 0.00000000
Convergence criteria met.
Covariance Parameter
Estimates
Cov Parm Estimate
blk 2.0411
Residual 33.3815
Fit Statistics
-2 Res Log Likelihood 185.4
AIC (smaller is better) 189.4
AICC (smaller is better) 189.9
BIC (smaller is better) 190.3
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
pillow 8 16 57.59 <.0001
Here is SAS code with the parameters fixed at the lmer variance estimates:
proc mixed data=pillow;
class blk pillow;
model firmness=pillow;
random blk;
parms (2.6584e-10) (3.5398e+01) /noiter;
run;
The Mixed Procedure
Parameter Search
CovP1 CovP2 Variance Res Log Like -2 Res Log Like
2.66E-10 35.3980 35.3981 -92.6996 185.3991
Covariance Parameter
Estimates
Cov Parm Estimate
blk 2.66E-10
Residual 35.3981
Fit Statistics
-2 Res Log Likelihood 185.4
AIC (smaller is better) 185.4
AICC (smaller is better) 185.4
BIC (smaller is better) 185.4
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
pillow 8 16 56.60 <.0001
More information about the R-sig-mixed-models
mailing list