[R-sig-ME] correlation structure in nlme
Juliette Chamagne
juliettechamagne at gmail.com
Wed Jun 29 11:00:45 CEST 2011
That's weird, no attachment. Let me try again
-------------- next part --------------
A non-text attachment was scrubbed...
Name: individual_fit.pdf
Type: application/pdf
Size: 121233 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110629/50ddf31c/attachment-0004.pdf>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: data_subset.txt
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20110629/50ddf31c/attachment-0004.txt>
-------------- next part --------------
Le Jun 28, 2011 ? 5:25 PM, Juliette Chamagne a ?crit :
> Hi all,
>
>
> I have a dataset of dbh increment (Cum.TRW) against Age. Here is attached a subset of my data, with only 6 trees.
> I performed a very simple nlme (no fixed or random effects), estimating the parameters of a monomolecular growth model based on log(Cum.TRW), for every individual.
> I then looked at the residuals of this model, and plotted acf and pacf on those. There is a strong autocorrelation at lag 1, for every individual.
> I thought I would put a correlation structure into my model. So I made the exact same model, but with the corAR1(). The AIC is much better, but when I look at the fit (I plot real data super imposed with both models: without correlation structure in green and with in red).
>
> When I calculate myself the autocorrelation coefficients for the 6 individuals and take an average, and then force it into CorAR1(), R crashes (see code below).
>
> Any ideas why is the fit with autocorrelation better based on AIC? And why does R crash when I fix the correlation parameter?
>
>
>
>
>
>
> library(nlme)
>
> sub <- read.table("data_subset.txt", h=T)
>
> model.data <- groupedData(log(Cum.TRW) ~ Age | IndID, data=sub)
>
> fit.nlme.NoCor <- nlme(log(Cum.TRW) ~ SSasymp (Age, Asym, M0, lrc), fixed=Asym + lrc + M0 ~ 1, data = model.data, random=(Asym + lrc ~ 1), verbose=T, start=c(4.7, -3, 2.5))
>
> fit.nlme.CorAR <- nlme(log(Cum.TRW) ~ SSasymp (Age, Asym, M0, lrc), fixed=Asym + lrc + M0 ~ 1, data = model.data, random=(Asym + lrc ~ 1), verbose=T, start=c(4.7, -3, 2.5), correlation=corAR1())
>
> fit.nlme.FixCorAR <- nlme(log(Cum.TRW) ~ SSasymp (Age, Asym, M0, lrc), fixed=Asym + lrc + M0 ~ 1, data = model.data, random=(Asym + lrc ~ 1), verbose=T, start=c(4.7, -3, 2.5), correlation=corAR1(0.9398526, fixed=TRUE))
>
>
>
> Thanks a lot,
>
> Juliette Chamagne
>
>
> --
> PhD student
> Institute for Evolution Biology and Environmental Studies
> University of Z?rich
> Winterthurerstrasse 190
> 8057 Z?rich, Switzerland
> office: +41 (0)44 635 61 21
>
> _______________________________________________
> 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