[R-sig-ME] Temporal correlation structure in mixed model

[Worthington, Thomas A]

Dear Ben 

Thank you for your suggestions  in answer to your questions

A trial is multiple fish but it is impossible to tell individuals apart, therefore the dependent variable is for the group. 

Yes the ten levels relates to an overall trial length of 10 minutes and yes there are a total of 200 observations

The typical number of approaches varies from 0.2 - 0.5 for setup 0 and 0.2 - 1.7 for setup 1, this value is corrected for fish that past the barrier. >From examination of the data the number of approaches is only different between the two setups for minutes 1, 2 and possibly 3.  The rational for treating minute as a categorical variable  was determine if there was an attraction to the barrier due to setup in those first minutes, which I was going to test using planned contrast once the overall model had been specified. WE had toyed with a simpler approach of multiple t-tests the correcting for multiple comparisons

In terms of the model you suggested maybe it would be preferable to set the variance structure to weights = varIdent(form=~ 1| Minute) as the difference in variance is mostly between minutes? I avoided transforming the data as it was suggested that it might be preferable to model the variance structure in Zuur et al 2009 Mixed Effects Models and Extensions in Ecology with R. 
  
test4b<-lme(Approaches ~ Setup * Minute, data=data, random = ~ 1 | Trial, weights = varIdent(form=~ 1| Minute), correlation = corAR1(form =~ Minute|Trial), control=list(opt="optim"))

I have compared the 5 models with AIC and BIC and 4b seems to fit the best (Test 4a being the one you put forth as a suggestion)

       df       AIC
test1  25 -40.44546
test2  27 -32.25960
test3  28 -38.20182
test4a 26 -42.20184
test4b 16 -20.23438

and 

       df      BIC
test1  25 41.50741
test2  27 56.24950
test3  28 53.58539
test4a 26 43.02914
test4b 16 32.21545

Best wishes

Tom 

-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com] 
Sent: Thursday, March 27, 2014 12:58 PM
To: Worthington, Thomas A; r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Temporal correlation structure in mixed model

On 14-03-27 12:23 PM, Worthington, Thomas A wrote:
> Dear All

> I'm trying to run a mixed model on ecology data. The setup is as
  follows, a group of fish are placed in a flume and attempt to pass a
  barrier. The dependent variable is the number of approaches to pass
  the barrier per minute of the trial, with each trial lasting 10
  minutes. There are two configurations of the barrier.  There are 20
  trials in total, 10 for each setup. The aim of the study is to see
  whether the number of approaches in minute 1, minute 2, minute 3
  etc.  is different between the two setups. Therefore my model is of
  the form Approaches ~ Setup + Minute + Setup * Minute, where minute
  (10 levels) and setup (2 levels) are both factors.

The dependent variable sounds more like a count than a continuous variable, so you might consider a Poisson GLMM, _but_ it is harder to incorporate temporal autocorrelation in that case, so if it seems reasonable you might try to stick with the linear mixed model (i.e.
Gaussian/Normal responses).  How big is your data set (is a 'trial'
equivalent to a single fish?  10 levels=10 minutes?  Does that you mean you have a total of 200 observations?)  What is a typical number/range of attempts per minute?


> I have three potential models and am slightly confused over the
  model specification, I believe I need to include a correlation
  structure to deal with temporal autocorrelation but am unclear over
  whether a random effect is needed instead or both are required. I
  have included a variance structure to deal with the issue of
  homogeneity of variance in the data.

A Poisson GLMM might help here.

> The three models are as follows



> test1 <- gls(Approaches ~ Setup * Minute, data=data, correlation =
  corAR1(form =~ Minute|Trial), weights = varIdent(form=~ 1| Minute *
  Setup), control=list(opt="optim"))

This will estimate correlation, but will assume that the average number of attempts is the same for every trial within a setup category.


> test2 <- lme(Approaches ~ Setup * Minute, data=data, random = ~
  Minute | Trial, weights = varIdent(form=~ 1| Minute * Setup),
  control=list(opt="optim"))

  This looks overspecified (does it really work??) -- it is allowing the differences among every pair of minutes to be characterized for every trial.  Even ~1|Trial/Minute is overspecified, since there is only one observation per minute per trial (thus the among-minute-within-trial variance would be confounded with the residual variance).

> test3 <- lme(Approaches ~ Setup * Minute, data=data, random = ~
  Minute | Trial, weights = varIdent(form=~ 1| Minute * Setup),
  correlation = corAR1(form =~ Minute|Trial),
  control=list(opt="optim"))

This looks even more overspecified.  I would suggest

lme(Approaches ~ Setup * Minute, data=data, random = ~ 1 | Trial,
   weights = varIdent(form=~ 1| Minute * Setup),
  correlation = corAR1(form =~ Minute|Trial),
  control=list(opt="optim"))

 which allows for variation among trials (1|Trial) and autocorrelation within Trial -- the Minute:Trial interaction will show up in the residual variance.

  This still looks ambitious if you have a total of 200 data points, as the varIdent() specification will estimate variances for all 20 of the minute*setup combinations.
  * Would log- or square-root transforming your data take care of the heteroscedasticity adequately?
  * Are you sure you want to treat minute as categorical, estimating a separate response for every minute?  Or would a smooth trend (linear, quadratic, or spline) describe the data adequately and more parsimoniously?  (This applies to both the mean and variance
description)


> I believe (maybe wrongly) that Test 1 and 2 would be broadly
  similar, but if some could suggest the most sensible model I would
  greatly appreciate it

> Best wishes
> 
> Tom
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>