[R-sig-ME] Need some advice on my model specification
Kingsford Jones
kingsfordjones at gmail.com
Sun Jul 12 22:01:15 CEST 2009
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