Hi Thierry,
thank you for your answer.
What I am looking for is an interaction between sp and Tank in slope
(perhaps this is not what I modelled, but it was what I was trying to
model).
I am OK if the strength of the correlation is equal among groups (in this
case, subjects), I thought that including a correlation structure was to
correct for non-independent residuals? (am I incorrect?)
mat.
On Tue, Apr 12, 2016 at 9:37 AM, Thierry Onkelinx <
thierry.onkelinx at inbo.be> wrote:
Dear Mathew,
AFAIK you can't model different correlation parameters with nlme. The
correlation structures only allow to specify the grouping and the time
covariate. The strength of the correlation is assumed to be equal among
groups.
In your case you are looking for correlation structures with different
strength among the groups. So I would suggest that you look for a technique
which allows you to model that. Then you can compare a model with different
correlation strengths with a model with equal correlation strength.
Best regards,
Thierry
ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature
and Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium
To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey
2016-04-11 15:30 GMT+02:00 Mathew Vickers <vickers.mathew at gmail.com>:
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]]