[R-sig-ME] correlation structure in nlme

[Juliette Chamagne]

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