convergence error code in mixed effects models
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