[R-sig-ME] Apparent false convergence in lmer with some BIBD data
Douglas Bates
bates at stat.wisc.edu
Thu Oct 29 21:40:29 CET 2009
On Thu, Oct 29, 2009 at 2:05 PM, Phillip Chapman
<pchapman at stat.colostate.edu> wrote:
> Greetings ME Board members,
There's an ME Board? I didn't know.
> 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.
The news that readers of this list have come to fear from me are, "I'm
in the midst of reformulating the lmer function" and, well, I am. In
the formulation being developed in the lme4a package under the R-forge
site, lme4.r-forge.r-project.org there is an option for the lmer
function to return an object that has a function to evaluate the
deviance for various choices of the relative variance parameters. In
this model there is only one relative variance parameter which is
expressed as the ratio of sigma_b (the standard deviation of the
random effects) to sigma (the residual standard deviation).
What is happening here is that the likelihood surface is unusually
flat so a wide range of values provide a likelihood close to that at
the mle. Consider the enclosed plot of the profiled deviance as a
function of theta, this relative standard deviation parameter.
Because this is the profiled deviance we would compare the change in
the profiled deviance (the difference between the value shown here and
the minimum) to a chisquared distribution with 11 degrees of freedom
(from 9 fixed-effects parameters and two variance components). A 95%
joint confidence region would include all the values shown here and
many more. The difference between the value at theta = 0 and at the
optimal theta is negligible.
You are correct that the nlminb optimizer used in the released value
of lmer produces a false optimum for this problem but, as I said, the
difference in the profiled deviance is negligible.
> 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
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
R version 2.11.0 Under development (unstable) (2009-10-28 r50245)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
Natural language support but running in an English locale
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> sessionInfo()
R version 2.11.0 Under development (unstable) (2009-10-28 r50245)
x86_64-unknown-linux-gnu
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=C LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
> library(lme4a)
Loading required package: stats4
Loading required package: Matrix
Loading required package: lattice
Loading required package: splines
> pillow <- within(read.table("./pillow.txt", header = TRUE),
+ blk <- factor(blk))
> h <- lmer(firmness ~ pillow + (1|blk), pillow, doFit = FALSE)
> (opt <- optimize(h at setPars, c(0,1)))
$minimum
[1] 0.2472657
$objective
[1] 185.3518
> h at setPars(0)
[1] 185.3991
> h at setPars(0.5)
[1] 185.6404
> thvals <- seq(0, 1, 0.01)
> xyplot(sapply(thvals, h at setPars) ~ thvals, type = "l",
+ xlab = expression(theta), ylab = "Profiled deviance")
>
> proc.time()
user system elapsed
12.660 0.150 13.283
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Rplots.pdf
Type: application/pdf
Size: 5147 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20091029/d50277c9/attachment.pdf>
More information about the R-sig-mixed-models
mailing list