Skip to content

convergence error code in mixed effects models

7 messages · Ilona Leyer, Douglas Bates, Mark Difford

Ilona Leyer
# 
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:
[1] "rep"    "treat"  "leaf"   "week"
Error in lme.formula(leaf~ treat, random = ~week|rep)
: 
        nlminb problem, convergence error code = 1;
message = iteration limit reached without convergence
(9)

Has anybody an idea to solve this problem?
Thanks

Ilona

Ilona Leyer
Conservation Biology
University of Marburg
Germany
e-mail: ileyer at yahoo.de
Douglas Bates
#
On Dec 13, 2007 4:15 PM, Ilona Leyer <ileyer at yahoo.de> wrote:
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?
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
Ilona Leyer
# 
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:
[1] "week"      "rep"       "treat"     "heightfra"
"leaffra"   "leafvim"
test<-groupedData(heightfra~week|rep,outer=~treat,test)
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)
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)
Error in summary(model) : object "model" not found
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:
Douglas Bates
# 
Thank you for sending the data.  It is very helpful in understanding
the nature of the problem.

For example, in your original description of your study you referred
to week as a factor, which is a completely reasonable term, but I
mistakenly thought that you meant an object of class "factor" and that
was why I replied that you would be estimating too many variances and
covariances.

I can tell you why you are having problems fitting a mixed-effects
model.  Strangely it is because there is too little variability in the
patterns across the replicates, especially at the early times.  The
leaf number is discrete with a small range (all the leaffra
observations in this example are 4, 5, 6, 7 or 8) and non-decreasing
over time.  (I assume the nature of the experiment is such that the
leaf number is necessarily non-decreasing.)  That doesn't allow for
many patterns.  I'm sure some clever person reading this will be able
to tell us exactly how many different such patterns you could get but
I will simply say "not many".

Notice that the first 10 lines show that the leaffra is 4 at week 4 in
*every* replicate.  There isn't a whole lot of variation here for the
random effects to model.

The best place to start is with a plot of the data.  I changed the
levels of the rep factor to "f1"-"f5" and "s1"-"s5" to indicate that
each rep is at one level of the treat. (Those who are playing along at
home should be careful of the ordering of the original levels because
ID10, which I now call s5, occurs between ID1 and ID2.)

With this set of labels the cross-tabulation and treat should be
treat
rep  pHf pHs
  f1   4   0
  f2   4   0
  f3   4   0
  f4   4   0
  f5   4   0
  s1   0   4
  s2   0   4
  s3   0   4
  s4   0   4
  s5   0   4

Now look at lattice plots such as

library(lattice)
xyplot(heightfra ~ week | rep, leaf, type = c("g", "p", "r"), layout =
c(5,2), aspect = 'xy', groups = treat)

(I enclose PDF files of these plots for each of the three responses.)
First you can see that there is very little variation at the low end.
Strangely enough, this causes a problem in fitting mixed-effects
models because the mle's of the variances of  the random effects for
the intercept will be zero.  The lme function does not handle this
gracefully.  The lmer function from the lme4 package does a better job
on this type of model.

Also, note that the pattern of heightfra over time is not linear.  It
is consistently concave down.  Thuso a mixed-effects model that is
linear in week will miss much of the structure in the data.

The point of R is to encourage you to explore your data rather than
subjecting it to a "canned" analysis.  You could try fitting a
mixed-effects model to these data in SAS PROC MIXED or SPSS MIXED and
I have no doubt that those packages would give you estimates (not to
mention p-values, something that the author of lmer has been woefully
negligent in not providing :-) but you probably won't get much of a
hint that the model doesn't make sense.  I would prefer to start with
the plot and see what the data have to say.


The technical problem with convergence in lme is that the mle of the
variance of the intercept term is zero.  You can see that if you use
lmer from the lme4 package instead to fit the model.
the random effect for the intercept is estimate
On Dec 14, 2007 4:40 AM, Ilona Leyer <ileyer at yahoo.de> wrote:
-------------- next part --------------
A non-text attachment was scrubbed...
Name: hf.pdf
Type: application/pdf
Size: 22554 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20071214/e8cff557/attachment.pdf 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lf.pdf
Type: application/pdf
Size: 23715 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20071214/e8cff557/attachment-0001.pdf 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: lv.pdf
Type: application/pdf
Size: 23776 bytes
Desc: not available
Url : https://stat.ethz.ch/pipermail/r-help/attachments/20071214/e8cff557/attachment-0002.pdf
Mark Difford
# 
Hi Ilona,
If there is, then Professor Bates (the gentleman who replied to your
question) will have tried to find it, and fix it, for you.

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:

  
    

      
      
    

  


  


      

    

    
1 day later

    

      
      

        


  

        

      Douglas Bates
    

    

      
      #
    

  


  

      
On Dec 14, 2007 9:26 PM, Mark Difford <mark_difford at yahoo.co.uk> wrote:

        

            

      
      
      
    

        
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.

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Ilona Leyer
    

    

      
      #
    

  


  

      
Thank you for your helpful answer. The only thing
obout what I still wonder is that some of the analyses
work with R some time and some time not without any
changes in the data or commands.
Your text (from both emails) is cut at the end, so
that  your information about the convergence error
problem is  not complete. Would you send it again?
Than you for your help!

Ilona


--- Douglas Bates <bates at stat.wisc.edu> schrieb:

            

      
      
      
    

        
test<-groupedData(heightfra~week|rep,outer=~treat,test)

            

      
      
      
    

        
test<-groupedData(leaffra~week|rep,outer=~treat,test)

            

      
      
      
    

        
=== message truncated ===