[R-sig-ME] specify random term and autocorrelation and plot lme

[Mathew Vickers]

Hi,

I am modeling autocorrelation in a an individual's walk. For a given
walker, I estimate the autocorrelation decay as lag increases.

Dependent
corel (autocorrelation)

Factors:
2 species (a,b)
2 Tanks (t1, t2 - these are two experimental arenas)

random:
subject (there was 10 individuals in each species)

I want to ask whether there is a significant difference in "corel" as a
function of "lag" between species in each Tank, and between the two Tanks
per species


The data look not dissimilar to this:


# data for SIG question

set.seed(100)
dataz <- expand.grid(subject = 1:10, sp = c("a", "b"), Tank = c("t1",
"t2"), lag = 1:100)
dataz$subject <- as.factor(paste(dataz$subject, dataz$sp, (dataz$Tank),
sep=""))
slope <- c(rnorm(10, 40, 2), rnorm(10, 5, 2)); slope <- rep(slope, 200)
dataz$corel <- dataz$lag*slope + rnorm(4000, 0, 1)*dataz$lag
dataz$corel[dataz$Tank=="t2" & dataz$sp=="a"] <-
dataz$corel[dataz$Tank=="t2"& dataz$sp=="a"]*2
dataz$corel[dataz$Tank=="t1" & dataz$sp=="b"] <-
dataz$corel[dataz$Tank=="t1"& dataz$sp=="b"]*1.8
dataz$corel <- -1*dataz$corel
plot(corel~lag, data=dataz, col=interaction(sp, Tank),
pch=c(1,2)[as.numeric(Tank)], cex=0.5)
plot(corel~lag, data=dataz, col=as.numeric(subject),
pch=c(1,2)[as.numeric(Tank)], cex=0.5)


# a dataset with:
# corel = correlation score
# sp = species (a,b)
# Tank = experimental setup (t1, t2)
# lag = indepentdent variable. A time-step (1:100)
# subject = individual
# in this case, I want to force the intercept to 0. The response variable,
"corel" is the autocorrelation between points with lag "lag". This means
point 1 (lag 0) = 0

library(ggplot)
library(nlme)

lma <- gls(corel~0+lag, data=dataz)
ggplot(dataz, aes(x=lag, y=corel)) + geom_point()+
 stat_smooth(method="lm", size=1, se = T)
# this is a test model


lm0 <- gls(corel~0+sp*Tank*lag, data=dataz)
anova(lm0)
plot(lm0)

# add autocorrelation strucuture, allow it to vary by subject?
lm1 <- gls(corel~0+sp*Tank*lag, correlation = corAR1(form = ~ lag |
subject), data=dataz)
anova(lm0, lm1)
plot(lm1)  # i am not convinced autocor structure has helped
summary(lm1)

# is this a plot of model  lm1   ? There is no allowance for
autocorrelation:
ggplot(dataz, aes(x=lag, y=corel, group=interaction(sp, Tank))) +
  geom_point(aes(colour=interaction(sp, Tank)))+
  stat_smooth(method="lm", size=1, se = T)


# add random effect.
# i am not sure if the random effect is specified correctly: repeated
measured of corel per subject, subject nested in sp (species)?
# if random=~subject|sp, the autocorrelation and random terms are
incompatible
lm2 <- lme(corel~0+sp*Tank*lag, random=~1|subject, correlation =
corAR1(form = ~ lag | subject), data=dataz)
plot(lm2)
anova(lm2)
anova(lm0, lm1, lm2)  # seems like lm2 is the best

#how do I plot this?
# I want a plot like the ggplot above, but I am not sure how to do it.




-- 
Mathew Vickers
Post Doc
OuLaLab
CNRS
Moulis, France

	[[alternative HTML version deleted]]