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
correlation structure in nlme
2 messages · Juliette Chamagne
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.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.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