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

[Mathew Vickers]

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]]
>>
>> _______________________________________________
>> R-sig-mixed-models at 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]]