[R] convergence error code in mixed effects models
Douglas Bates
bates at stat.wisc.edu
Fri Dec 14 21:27:18 CET 2007
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
> xtabs(~ rep + treat, leaf)
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:
>
> 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
> >
>
>
-------------- 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
More information about the R-help
mailing list