From: Douglas Bates <bates_at_stat.wisc.edu>

Date: Sun, 16 Dec 2007 12:20:51 -0600

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. Received on Sun 16 Dec 2007 - 18:25:12 GMT

Date: Sun, 16 Dec 2007 12:20:51 -0600

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.
**>
*

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. Received on Sun 16 Dec 2007 - 18:25:12 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sun 16 Dec 2007 - 20:30:20 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*