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

[Thierry Onkelinx]

Dear Mathew,

My point was not to estimate the autocorrelation and model these estimated
autocorrelations. But rather model the random walk with different
autocorrelation structures and compare these models.

Modelling the estimated autocorrelation at different lags inflates the data
and ignores the uncertainty on the autocorrelation measurements.

Best regards,

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-12 14:42 GMT+02:00 Mathew Vickers <vickers.mathew op gmail.com>:

> 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 op 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 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
>>>
>>
>>
>
>
> --
> Mathew Vickers
> Post Doc
> OuLaLab
> CNRS
> Moulis, France
>

	[[alternative HTML version deleted]]