[R] convergence error code in mixed effects models

Douglas Bates bates at stat.wisc.edu
Sun Dec 16 19:20:51 CET 2007


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



More information about the R-help mailing list