[R-sig-ME] nlme optim control option
Stephan Moratti
moratti at med.ucm.es
Tue Apr 29 16:48:34 CEST 2008
Thank you very much for your fast response. I will first describe the
data and then the models applied. Further, I will attach an .RData
object with all the data, so you can directly use it if you like.
The data are amplitude values obtained from magnetoencephalografic data
(like EEG). There are 8 conditions and 12 subjects (I know, the number
of conditions is quite high with respect to number of subjects). The
basic data matrix looks like this (called "mat" in the RData object):
V1 V2 V3 V4 V5 V6 V7 V8
1 2.851582 3.314363 3.279070 2.394195 2.914685 1.876838 2.180715 3.203826
2 4.521648 6.171458 4.124039 5.010868 4.923190 5.221255 4.745925 5.477719
3 3.957449 4.743048 3.015077 2.890330 3.428396 3.078638 3.015924 2.921296
4 5.139640 5.272191 4.067972 3.151300 5.739077 3.279299 3.524670 6.110507
5 7.928020 8.344273 7.525576 6.835657 8.995732 8.871082 7.821745 8.346335
6 4.900010 4.508470 2.532844 4.197058 4.701989 2.878392 4.939342 2.913649
7 4.378253 7.150406 5.205363 4.872032 5.460943 5.901771 4.986221 4.985459
8 5.697526 4.945624 3.855982 4.460146 4.661727 5.434473 6.283447 5.795619
9 2.772407 2.823246 2.693566 3.019434 2.303198 3.340165 3.624334 2.974910
10 5.063280 6.384693 6.015056 7.935250 5.564225 7.541611 7.154139 8.031693
11 3.202468 1.877924 1.716837 2.286963 2.152085 3.096542 3.026434 4.369565
12 4.684218 5.391280 5.360308 4.675992 4.500356 3.507834 3.795734 7.622723
The same data matrix is transformed to a data frame ("datamat" in the
RData object):
values ind subject condition
1 2.851582 V1 1 1
2 4.521648 V1 2 1
3 3.957449 V1 3 1
4 5.139640 V1 4 1
5 7.928020 V1 5 1
6 4.900010 V1 6 1
7 4.378253 V1 7 1
8 5.697526 V1 8 1
9 2.772407 V1 9 1
10 5.063280 V1 10 1
11 3.202468 V1 11 1
12 4.684218 V1 12 1
13 3.314363 V2 1 2
14 6.171458 V2 2 2
15 4.743048 V2 3 2
16 5.272191 V2 4 2
17 8.344273 V2 5 2
18 4.508470 V2 6 2
19 7.150406 V2 7 2
20 4.945624 V2 8 2
21 2.823246 V2 9 2
22 6.384693 V2 10 2
23 1.877924 V2 11 2
24 5.391280 V2 12 2
.... up to 8 conditions
**********************
* Applying optim resutls in: *
**********************
fm.optim<-lme(values~condition,data=datamat,random=~1|subject,correlation=corSymm(),control=list(msMaxIter=100,opt
= "optim",msVerbose=TRUE))
initial value 209.185062
iter 10 value 186.050690
iter 20 value 184.942325
iter 30 value 184.724836
iter 40 value 184.121038
iter 50 value 183.918251
final value 183.903769
converged
summary(fm.optim)
Linear mixed-effects model fit by REML
Data: datamat
AIC BIC logLik
299.5351 393.6739 -111.7675
Random effects:
Formula: ~1 | subject
(Intercept) Residual
StdDev: 1.110893 1.176853
Correlation Structure: General
Formula: ~1 | subject
Parameter estimate(s):
Correlation:
1 2 3 4 5 6 7
2 -0.056
3 -0.432 0.742
4 -0.213 0.545 0.586
5 0.729 0.593 0.281 0.235
6 0.101 0.474 0.412 0.657 0.430
7 0.625 0.054 -0.212 0.440 0.509 0.635
8 -0.087 0.261 0.557 0.417 0.257 0.354 0.005
Fixed effects: values ~ condition
Value Std.Error DF t-value p-value
(Intercept) 4.591375 0.4671782 77 9.827888 0.0000
condition2 0.485873 0.4936256 77 0.984295 0.3281
condition3 -0.475401 0.5748486 77 -0.827002 0.4108
condition4 -0.280606 0.5290556 77 -0.530391 0.5974
condition5 0.020759 0.2498983 77 0.083068 0.9340
condition6 -0.089050 0.4556203 77 -0.195448 0.8456
condition7 0.000178 0.2941852 77 0.000604 0.9995
condition8 0.638067 0.5009446 77 1.273727 0.2066
Correlation:
(Intr) cndtn2 cndtn3 cndtn4 cndtn5 cndtn6 cndtn7
condition2 -0.528
condition3 -0.615 0.907
condition4 -0.566 0.801 0.846
condition5 -0.267 0.860 0.790 0.627
condition6 -0.488 0.734 0.768 0.847 0.608
condition7 -0.315 0.385 0.406 0.762 0.243 0.783
condition8 -0.536 0.655 0.832 0.748 0.566 0.678 0.366
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.9205059 -0.7430723 -0.1023267 0.7054844 2.1579685
Number of Observations: 96
Number of Groups: 12
***********************
* Applying nlminb results in: *
***********************
fm.nlminb<-lme(values~condition,data=datamat,random=~1|subject,correlation=corSymm(),control=list(msMaxIter=100,opt
= "nlminb",msVerbose=TRUE))
....
98: 183.73454: 0.127036 0.0579634 0.506547 -1.36663 0.306432
-0.838546 -0.395783 -1.18434 -2.01708 -1.91615 0.0747060 -0.199496
-0.740317 -0.398701 -0.835580 0.623620 -0.979961 -0.232816 -0.0991876
-1.63724 2.08835 -12.9687 0.147038 -0.361111 -0.964921 -0.271664
-0.512183 -0.409421 -2.67278
99: 183.73454: 0.127040 0.0579607 0.506550 -1.36663 0.306427
-0.838537 -0.395789 -1.18434 -2.01709 -1.91611 0.0747092 -0.199503
-0.740320 -0.398707 -0.835577 0.623649 -0.979960 -0.232836 -0.0991891
-1.63723 2.08842 -12.9657 0.147034 -0.361111 -0.964925 -0.271658
-0.512177 -0.409425 -2.67234
100: 183.73454: 0.127046 0.0579570 0.506540 -1.36663 0.306412
-0.838550 -0.395794 -1.18434 -2.01711 -1.91609 0.0747311 -0.199514
-0.740328 -0.398712 -0.835571 0.623651 -0.979969 -0.232830 -0.0992005
-1.63723 2.08843 -12.9642 0.147026 -0.361116 -0.964933 -0.271660
-0.512178 -0.409420 -2.67211
100: 183.73454: 0.127046 0.0579570 0.506540 -1.36663 0.306412
-0.838550 -0.395794 -1.18434 -2.01711 -1.91609 0.0747311 -0.199514
-0.740328 -0.398712 -0.835571 0.623651 -0.979969 -0.232830 -0.0992005
-1.63723 2.08843 -12.9642 0.147026 -0.361116 -0.964933 -0.271660
-0.512178 -0.409420 -2.67211
Error en lme.formula(values ~ condition, data = datamat, random = ~1 | :
nlminb problem, convergence error code = 1
message = iteration limit reached without convergence (9)
using more iterations:
fm.nlminb2<-lme(values~condition,data=datamat,random=~1|subject,correlation=corSymm(),control=list(msMaxIter=200,opt
="nlminb",msVerbose=TRUE))
....
98: 183.73454: 0.127036 0.0579634 0.506547 -1.36663 0.306432
-0.838546 -0.395783 -1.18434 -2.01708 -1.91615 0.0747060 -0.199496
-0.740317 -0.398701 -0.835580 0.623620 -0.979961 -0.232816 -0.0991876
-1.63724 2.08835 -12.9687 0.147038 -0.361111 -0.964921 -0.271664
-0.512183 -0.409421 -2.67278
99: 183.73454: 0.127040 0.0579607 0.506550 -1.36663 0.306427
-0.838537 -0.395789 -1.18434 -2.01709 -1.91611 0.0747092 -0.199503
-0.740320 -0.398707 -0.835577 0.623649 -0.979960 -0.232836 -0.0991891
-1.63723 2.08842 -12.9657 0.147034 -0.361111 -0.964925 -0.271658
-0.512177 -0.409425 -2.67234
100: 183.73454: 0.127046 0.0579570 0.506540 -1.36663 0.306412
-0.838550 -0.395794 -1.18434 -2.01711 -1.91609 0.0747311 -0.199514
-0.740328 -0.398712 -0.835571 0.623651 -0.979969 -0.232830 -0.0992005
-1.63723 2.08843 -12.9642 0.147026 -0.361116 -0.964933 -0.271660
-0.512178 -0.409420 -2.67211
101: 183.73454: 0.127045 0.0579574 0.506540 -1.36663 0.306412
-0.838549 -0.395794 -1.18434 -2.01711 -1.91609 0.0747313 -0.199515
-0.740328 -0.398712 -0.835570 0.623651 -0.979969 -0.232834 -0.0991995
-1.63723 2.08843 -12.9642 0.147026 -0.361116 -0.964934 -0.271660
-0.512178 -0.409419 -2.67211
102: 183.73454: 0.127045 0.0579574 0.506540 -1.36663 0.306412
-0.838549 -0.395794 -1.18434 -2.01711 -1.91609 0.0747313 -0.199515
-0.740328 -0.398712 -0.835570 0.623651 -0.979969 -0.232834 -0.0991995
-1.63723 2.08843 -12.9642 0.147026 -0.361116 -0.964934 -0.271660
-0.512178 -0.409419 -2.67211
Error en lme.formula(values ~ condition, data = datamat, random = ~1 | :
nlminb problem, convergence error code = 1
message = false convergence (8)
However, if I log transform the data, both optimization methods converge
using a maxIter of 200. All objects could be found in the Rdata object.
Not to blow up the mail, I did not print the summary of the two fitted
models as they are in the RData object. The two methods converge at
similar values using the log transformed data. I don't know why "nlminb"
does not converge with the not transformed data.
Best,
Stephan
Douglas Bates escribió:
> On 4/29/08, Doran, Harold <HDoran at air.org> wrote:
>
>> I actually trust nlminb more than I trust optim. In some of the functions I have written in the MiscPsycho package (e.g., irt.ability) we get multimodal distributions for some of the item response theory models.
>>
>> When I use optim for optimization, you can get very funny results with difficult to maximize distributions. nlminb, however, gave reasonable results and the maximum was always confirmed to be "correct" via a visual examination of the likelihoods.
>>
>> In other words, I can plot the likelihood function and get a sense of where the max is. optim would more often than not give results that are very (very) far away from the max with multimodal distributions whereas nlminb always gave back the max that seemed consistent with the visual plot of the likelihood.
>>
>> With all of this said, you should switch to lme4.
>>
>
> As Harold indicated, switching the optimizer used in the nlme package
> from optim to nlminb was intentional. I had seen cases of spurious
> convergence for optim and found it easier to introduce the nlminb
> optimizer.
>
> Could you give more details of the cases where functions from the nlme
> package using nlminb are not converging? Are you using the lme
> function or the nlme function? Does your model have random effects
> and fixed effects only or does it also use parameterized correlation
> structures?
>
> As Harold indicated, the lme4 package is now the preferred way to fit
> linear or nonlinear mixed-effects models that have only fixed effects
> and random effects. If you can show us one of your calls to lme or
> nlme we can tell you how to translate it to a call to lmer or nlmer.
>
>
>> > -----Original Message-----
>> > From: r-sig-mixed-models-bounces at r-project.org
>> > [mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf
>> > Of Stephan Moratti
>> > Sent: Tuesday, April 29, 2008 7:05 AM
>> > To: r-sig-mixed-models at r-project.org
>> > Subject: [R-sig-ME] nlme optim control option
>> >
>> > Hello,
>> >
>> > I am using the nlme package for quite a while, but recently I
>> > have discovered the "optim" option in the lmeControl function
>> > and I have read that until the R verison 2.2.0 the
>> > optimization function is "nlminb" and before it was "optim"
>> > by default. As I am a "user" of R I am not into the
>> > optimization functions. However, I have a case now where with
>> > optim = "nlminb" the model does not converge, whereas with
>> > optim = "optim" the model converges.
>> >
>> > Can somebody explain to me in easy words what the difference
>> > is ? If nlminb does not converge, but optim does, is the
>> > result less trustable ?
>> >
>> > Thanks,
>> >
>> > Stephan
>> >
>> >
>> > --
>> > *Stephan Moratti, PhD/
>> > /**/see also: http://web.mac.com/smoratti/ /*Centro de
>> > Tecnología Biomédica CBT, Universidad Politécnica de Madrid,
>> >
>> > en la actualidad (currently at) en el
>> > Centro de Magnetoencefalografía Dr. Perez Modrego,
>> > Universidad Complutense de Madrid, Faculdad de Medicina,
>> > Pabellón 8, Avda. Complutense, s/n, 28040 Madrid, Spain,
>> >
>> > email: moratti at gbt.tfo.upm.es <mailto:moratti at gbt.tfo.upm.es>
>> > moratti at med.ucm.es
>> > Tel.: +34 91 394 2292
>> > Fax.: +34 91 394 2294
>> > */
>> > /*
>> >
>> > _______________________________________________
>> > R-sig-mixed-models at r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> >
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>>
>
>
--
*Stephan Moratti, PhD/
/**/see also: http://web.mac.com/smoratti/
/*Centro de Tecnología Biomédica CBT,
Universidad Politécnica de Madrid,
en la actualidad (currently at) en el
Centro de Magnetoencefalografía Dr. Perez Modrego,
Universidad Complutense de Madrid,
Faculdad de Medicina,
Pabellón 8,
Avda. Complutense, s/n,
28040 Madrid,
Spain,
email: moratti at gbt.tfo.upm.es <mailto:moratti at gbt.tfo.upm.es>
moratti at med.ucm.es
Tel.: +34 91 394 2292
Fax.: +34 91 394 2294
*/
/*
More information about the R-sig-mixed-models
mailing list