[R-sig-ME] Just to say thanks to Kingsford
Djibril Dayamba
Djibril.Dayamba at ess.slu.se
Sun Aug 2 16:02:57 CEST 2009
Hello Kingsford,
I just want to say thanks to you for the guidance I have got in linear mixed effect modeling with my data. I believe there is much left to learn but at this point, I have made a lot of progress and I am happy with it.
Best wishes,
Djibril.
-----Original Message-----
From: Kingsford Jones [mailto:kingsfordjones at gmail.com]
Sent: den 17 juli 2009 00:17
To: Djibril Dayamba; R-Sig Mixed-models
Subject: Re: Problem with convergence, please help
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