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

Djibril Dayamba Djibril.Dayamba at ess.slu.se
Thu Jul 16 10:43:46 CEST 2009


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