convergence error code in mixed effects models
On Dec 14, 2007 9:26 PM, Mark Difford <mark_difford at yahoo.co.uk> wrote:
Is there a solution for this problem?
If there is, then Professor Bates (the gentleman who replied to your question) will have tried to find it, and fix it, for you.
What I was trying to indicate in my replies is that "this problem" is the user attempting to fit a model that is not appropriate for the data. I feel, and I hope that Jos? Pinheiro and I conveyed in our book, that analysis of longitudinal data should always begin with plots of the response versus time by experimental unit. Much of the wonderful work that Deepayan Sarkar did in developing the lattice package was motivated by the desire to produce exactly those plots. Because the purpose of the analysis is to look at the behavior within units over time and see how these patterns differ between units, it is clear that you should always begin with such a plot. It is disappointing to have users feel that the software is somehow inferior because it doesn't mindlessly produce estimates for inappropriate models when it is so simple for the user to check the patterns in the data and decide if the model is appropriate. As I said, I think that SAS and SPSS are better tools for fitting "the" repeated measures model or "the" longitudinal data model when you don't want to be bothered with actually looking at your data and thinking about the model. (And yes, I am being a trifle sarcastic in saying that. Please don't quote me as having endorsed the SAS and SPSS mixed model software.) By the way, I received a personal reply from a friend who asked what I meant when I said that the model fit by SAS PROC MIXED to these data may not make sense. I haven't used SAS PROC MIXED in a long time (I can never get past the "CARDS;" statement and still take the software seriously) but I strongly suspect that it would produce one of those mixed models that occurs in SAS computation and nowhere else. If you allow correlated random effects for the intercept and slope on these data you will probably get the ML or REML estimate of the variance of the intercept random effects being driven to zero while the estimate of the covariance of the intercept and slope random effects stays decidedly nonzero. Apparently mathematical impossibility is not an impediment to parameter estimation in such cases. In fact, Singer and Willett claim on p. 154 of their 2003 book "Applied Longitudinal Data Analysis" that one can go further and invoke an option to allow for negative variance estimates. The mind boggles.
Professor Bates wrote/co-wrote the software package (nlme) you are using. And while I have nothing against Crawley's book, you are usually much better off going to primary sources first, to solve this kind of problem (which, of course you have done, though may not have been aware of it ;) Mixed-Effects Models in S and S-PLUS, by: Pinheiro, Jos?, Bates, Douglas http://www.springer.com/west/home/statistics/computational?SGWID=4-10130-22-2102822-0 Hope this speeds you on your way... Regards, Mark. Ilona Leyer wrote:
Here an simple example: rep treat heightfra leaffra leafvim week ID1 pHf 1.54 4 4 4 ID2 pHf 1.49 4 4 4 ID3 pHf 1.57 4 5 4 ID4 pHf 1.48 4 4 4 ID5 pHf 1.57 4 4 4 ID6 pHs 1.29 4 5 4 ID7 pHs 0.97 4 5 4 ID8 pHs 2.06 4 4 4 ID9 pHs 0.88 4 4 4 ID10 pHs 1.47 4 4 4 ID1 pHf 3.53 5 6 6 ID2 pHf 4.08 6 6 6 ID3 pHf 3.89 6 6 6 ID4 pHf 3.78 5 6 6 ID5 pHf 3.92 6 6 6 ID6 pHs 2.76 5 5 6 ID7 pHs 3.31 6 7 6 ID8 pHs 4.46 6 7 6 ID9 pHs 2.19 5 5 6 ID10 pHs 3.83 5 5 6 ID1 pHf 5.07 7 7 9 ID2 pHf 6.42 7 8 9 ID3 pHf 5.43 6 8 9 ID4 pHf 6.83 6 8 9 ID5 pHf 6.26 6 8 9 ID6 pHs 4.57 6 9 9 ID7 pHs 5.05 6 7 9 ID8 pHs 6.27 6 8 9 ID9 pHs 3.37 5 7 9 ID10 pHs 5.38 6 8 9 ID1 pHf 5.58 7 9 12 ID2 pHf 7.43 8 9 12 ID3 pHf 6.18 8 10 12 ID4 pHf 6.91 7 10 12 ID5 pHf 6.78 7 10 12 ID6 pHs 4.99 6 13 12 ID7 pHs 5.50 7 8 12 ID8 pHs 6.56 7 10 12 ID9 pHs 3.72 6 10 12 ID10 pHs 5.94 6 10 12 I used the procedure described in Crawley?s new R Book. For two of the tree response variables (heightfra,leaffra) it doesn?t work, while it worked with leafvim (but in another R session, yesterday, leaffra worked as well...). Here the commands:
attach(test) names(test)
[1] "week" "rep" "treat" "heightfra" "leaffra" "leafvim"
library(nlme)
test<-groupedData(heightfra~week|rep,outer=~treat,test)
model1<-lme(heightfra~treat,random=~week|rep)
Error in lme.formula(heightfra ~ treat, random = ~week
| rep) :
nlminb problem, convergence error code = 1;
message = iteration limit reached without convergence
(9)
test<-groupedData(leaffra~week|rep,outer=~treat,test)
model2<-lme(leaffra~treat,random=~week|rep)
Error in lme.formula(leaffra ~ treat, random = ~week |
rep) :
nlminb problem, convergence error code = 1;
message = iteration limit reached without convergence
(9)
test<-groupedData(leafvim~week|rep,outer=~treat,test)
model3<-lme(leafvim~treat,random=~week|rep) summary(model)
Error in summary(model) : object "model" not found
summary(model3)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
129.6743 139.4999 -58.83717
Random effects:
Formula: ~week | rep
Structure: General positive-definite, Log-Cholesky
parametrization
StdDev Corr
(Intercept) 4.4110478 (Intr)
week 0.7057311 -0.999
Residual 0.5976143
Fixed effects: leafvim ~ treat
Value Std.Error DF t-value p-value
(Intercept) 5.924659 0.1653596 30 35.82893 0.0000
treatpHs 0.063704 0.2338538 8 0.27241 0.7922
Correlation:
(Intr)
treatpHs -0.707
Standardized Within-Group Residuals:
Min Q1 Med Q3
Max
-1.34714254 -0.53042878 -0.01769195 0.40644540
2.29301560
Number of Observations: 40
Number of Groups: 10
Is there a solution for this problem?
Thanks!!!
Ilona
--- Douglas Bates <bates at stat.wisc.edu> schrieb:
On Dec 13, 2007 4:15 PM, Ilona Leyer <ileyer at yahoo.de> wrote:
Dear All, I want to analyse treatment effects with time
series
data: I measured e.g. leaf number (five replicate plants) in relation to two soil pH - after 2,4,6,8 weeks. I used mixed effects models, but some
analyses
didn?t work. It seems for me as if this is a
randomly
occurring problem since sometimes the same model
works
sometimes not. An example:
names(test)
[1] "rep" "treat" "leaf" "week"
library (lattice) library (nlme)
test<-groupedData(leaf~week|rep,outer=~treat,test)
model<-lme(leaf~treat,random=~leaf|rep)
Error in lme.formula(leaf~ treat, random =
~week|rep) Really!? You gave lme a model with random = ~ leaf | rep (and no data specification) and it tried to fit a model with random = ~ week | rep? Are you sure that is an exact transcript?
:
nlminb problem, convergence error code =
1;
message = iteration limit reached without
convergence
(9)
Has anybody an idea to solve this problem?
Oh, I have lots of ideas but without a reproducible example I can't hope to decide what might be the problem. It appears that the model may be over-parameterized. Assuming that there are 4 different values of week then ~ week | rep requires fitting 10 variance-covariance parameters. That's a lot. The error code indicates that the optimizer is taking
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
-- View this message in context: http://www.nabble.com/convergence-error-code-in-mixed-effects-models-tp14325990p14340592.html Sent from the R help mailing list archive at Nabble.com.
______________________________________________ R-help at r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.