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

[Thierry Onkelinx]

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 op 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]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>

	[[alternative HTML version deleted]]