[R-sig-ME] Problem with convergence, please help

Kingsford Jones kingsfordjones at gmail.com
Fri Jul 17 00:17:01 CEST 2009


It's a complex model with algorithms trying to optimize many estimates
together:  random intercept and slope, residual variances within each
plot and AR1 correlation parameter for points across years, as well as
the 8? fixed effects.  Fitting residual covariance structures is
particularly tricky, and you're structuring diagonals and
off-diagonals, so it's not too surprising you've run into a
convergence problem.

As for solutions, take a look at the 'control' argument to lme, and
the lmeControl function.  At a minimum, try increasing the number of
iterations (maxIter=<something more than 50>) and set msVerbose=TRUE
to view progress during iterations.  If the algorithm is having
trouble with the varFunc or corStruct classes, you can use the 'value'
argument to specify starting values (e.g. correlation = corAR1(value =
.6, form = ~Year)).  You can also try simpler structures; e.g.,
weights = varPower(form =~fitted(.)) for the situation where variance
increases with the predicted response.

And of course P&B is the place to look for details...

hth,
Kingsford

On Thu, Jul 16, 2009 at 2:43 AM, Djibril
Dayamba<Djibril.Dayamba at ess.slu.se> wrote:
> Hello,
> Thanks a lot Kingsford for all your help. I went through the process step by step and my model kept improving. Fitting the Model with correlation structure corAR1(form=~Year) gave the smaller AIC but when checking the model there was a trend in the spread of the residuals. I then tried to account for the heterogeneity using varIdent as specified in the model below.
>
> lme(BA.130~Grazing*Fire*Year,random=~Year|Plot,correlation=corAR1(form=~Year), weights=varIdent(form=~1|Plot),na.action=na.omit)
>
> But this model could not work and I got an error message about convergence (see below).
>
> Error in lme.formula(BA.130 ~ Grazing * Fire * Year, random = ~Year |  :
>  nlminb problem, convergence error code = 1
>  message = iteration limit reached without convergence (9)
>
> Since then I have been reading but could not get the point. Could anyone tell me something about this problem and how it could be solved? Thanks in advance.
>
> With regards,
>
> Djibril.
>
>
>
> -----Original Message-----
> From: Kingsford Jones [mailto:kingsfordjones at gmail.com]
> Sent: den 12 juli 2009 22:01
> To: Djibril Dayamba
> Cc: r-sig-mixed-models at r-project.org
> Subject: Re: [R-sig-ME] Need some advice on my model specification
>
> Hi Djibril,
>
> The model is overspecified.  I strongly suspect a call of
> intervals(Model) will throw an error, or possibly produce intervals of
> infinite width.  A more likely candidate is:
>
> Model <- lme(BA.B ~ Grazing*Fire + Year, data = <your data>,
>    random = ~<random intercept and possible linear and quadratic slopes>|Plot,
>    correlation = corAR1(form=~Year|Plot), weights = <possible
> heterogeneity structure>)
>
> Without being familiar with the data this is of course only a suggestion.
>
> One strategy is to start with the fit with no correlation structure
> and explore residual temporal correlation within plots using the ACF
> function and its plot method, choosing a structure, fitting it, and
> repeating the exploration.  The same can be done for variance
> structures with the weights argument.  The following link (and P&B
> 2000) will suggest plots for exploration:
>
> http://bm2.genes.nig.ac.jp/RGM2/index.php?query=nlme
>
> As I mentioned in reply to your r-help post there are a variety of
> tricky issues surrounding inference; nonetheless, if you carefully
> build your linear model, checking that you are reasonably meeting
> normality assumptions (for errors and random effects), and using the
> lme arguments to structure error covariance so as to account for lack
> of independence and heterogeneity of variance, you will have results
> that should be quite defensible to reviewers.
>
> best of luck,
>
> Kingsford
>
>
> On Sun, Jul 12, 2009 at 11:43 AM, Djibril
> Dayamba<Djibril.Dayamba at ess.slu.se> wrote:
>> Hello,
>> I previously wrote to R-help, I got some advices and was kindly redirected by Kingsford to this specific Help-list (mixed models) for further questions. Since then I have moved a bit but I still have some points where I would appreciate having clarification.
>>
>> I have a factorial experiment ( with 4 repetitions for each treatment combination) to study the effects of Grazing and Fire on Forest biomass production. The experimental unit (to which the treatment combinations are applied) are PLOTs. The measures were made repeatedly for 13 years. Below is how I organized my data; Plot is the plot naming in the field; BA.B is the response variable (Basal area) I am using to express myself here
>>
>>
>> Grazing             Fire                     Plot                     Year                    BA.B
>> Ungrazed         Unburnt            102                      1                           398.13
>> Ungrazed         Unburnt            102                      2                           4728.54
>> Ungrazed         Unburnt            102                      3                           2092.05
>> Ungrazed         Unburnt            102                      4                           3076.70
>> Ungrazed         Unburnt            102                      5                           2578.54
>> Ungrazed         Unburnt            102                      6                           2541.07
>> Ungrazed         Unburnt            102                      7                           3191.61
>> Ungrazed         Unburnt            102                      8                           2526.75
>> Ungrazed         Unburnt            102                      9                           3665.42
>> Ungrazed         Unburnt            102                      10                        3077.42
>> Ungrazed         Unburnt            102                      11                        3911.63
>> Ungrazed         Unburnt            102                      12                        4067.28
>> Ungrazed         Unburnt            102                      13                        4457.94
>> Ungrazed         Unburnt            108                      1                           370.99
>> Ungrazed         Unburnt            108                      2                           2184.39
>> Ungrazed         Unburnt            108                      3                           2008.66
>> .
>> .
>> .
>> .
>> .
>>
>> I fitted the below model to account for the temporal autocorrelation and the variance heterogeneity. I also have a missing value.
>>
>> Model<-lme(BA.B~Grazing*Fire*Year, random=~1|Year/Plot/Fire/Grazing, correlation=corAR1(form=~Year), weights=varIdent(form=~1|Grazing*Fire*Year), na.action=na.omit)
>>
>> For the random effect I got something like this
>>
>> Random effects:
>>  Formula: ~1 | Year
>>        (Intercept)
>> StdDev:    234.6285
>>
>>  Formula: ~1 | Plot %in% Year
>>        (Intercept)
>> StdDev:    67.01272
>>
>>  Formula: ~1 | Fire %in% Plot %in% Year
>>        (Intercept)
>> StdDev:    66.80442
>>
>>  Formula: ~1 | Grazing %in% Fire %in% Plot %in% Year
>>        (Intercept) Residual
>> StdDev:    66.83408           153.0221
>>
>>
>> For the correlation structure, the Phi value is 0 (zero)
>>
>> Correlation Structure: AR(1)
>>  Formula: ~Year | Year/Plot/Fire/Grazing
>>  Parameter estimate(s):
>> Phi
>>  0
>>
>> For the fixed effects, please note the identical degree of freedom for every term (except Year)
>>
>> Fixed effects: BA.B ~ Grazing * Fire * Year
>>                                     Value                   Std.Error       DF                 t-value           p-value
>> (Intercept)                       461.4655          178.15612     396                2.590231       0.0099
>> GrazingUngrazed                  -180.1041   104.88596     396                -1.717143      0.0867
>> FireUnburnt                      -236.3691       124.54654     396                -1.897837      0.0584
>> Year                              410.9426             31.68234       11                  12.970718     0.0000
>> GrazingUngrazed:FireUnburnt 224.5099 153.17048    396                1.465752       0.1435
>> GrazingUngrazed:Year         -104.4984  28.99989       396                -3.603406      0.0004
>> FireUnburnt:Year                  -21.3375     38.52046       396                -0.553927      0.5799
>> GrazingUngrazed:FireUnburnt:Year 85.2475  44.47663   396           1.916680       0.0560
>>
>>
>> My concern is:
>> Does my model specification looks ok?
>> Is Phi = 0 real and what does this mean (in all case studies I have seen, it had value different from zero)?
>> I am also wondering about the identical degree of freedom for all term (except for Year).
>> Thanks in advance
>>
>>
>> With regards,
>>
>> Djibril.
>>
>>
>>
>>
>>        [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>
>




More information about the R-sig-mixed-models mailing list