[R] cca with repeated measures

Gavin Simpson gavin.simpson at ucl.ac.uk
Fri Nov 18 15:47:20 CET 2011


On Fri, 2011-11-18 at 15:17 +0100, René Mayer wrote:
> Thanks a lot Gavin!,
> this was what I was looking for.
> Have I got this right that with no 'cyclic shifts *within* strata' you  
> mean that I cannot define a nesting within animal, e.g.,  
> animal/year/season (speaking in regression-terms  random-effects for  
> the animal-specific season and year variation).

In the permutation framework, you condition the permutation on animal
(`strata = animal` in the anova() method), but that alone would say that
observations are freely exchangeable within each animal, but not
exchangeable between animals. As you have time series data then the
observations are not really freely exchangeable within animal either. We
can try to maintain the within-animal temporal dependence structure by
using cyclic shift permutations:

> require(permute)
> shuffleSeries(1:10)
 [1]  4  5  6  7  8  9 10  1  2  3
> t(replicate(10, shuffleSeries(1:10)))
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   10    1    2    3    4    5    6    7    8     9
 [2,]    1    2    3    4    5    6    7    8    9    10
 [3,]    3    4    5    6    7    8    9   10    1     2
 [4,]   10    1    2    3    4    5    6    7    8     9
 [5,]    7    8    9   10    1    2    3    4    5     6
 [6,]    8    9   10    1    2    3    4    5    6     7
 [7,]    9   10    1    2    3    4    5    6    7     8
 [8,]    3    4    5    6    7    8    9   10    1     2
 [9,]    4    5    6    7    8    9   10    1    2     3
[10,]    5    6    7    8    9   10    1    2    3     4

cyclic shifts basically take a random starting point for the permuted
timeseries and then keep sample in the order they appear in the data,
bending the ends of the series round together into a circle. The above
example:

> shuffleSeries(1:10)
 [1]  4  5  6  7  8  9 10  1  2  3

shows a single cyclic shift of the observations 1:10 in an equally
spaced time series. Here the random starting point was observation 4 and
note that we join the ends of the timeseries so after observation 10, we
have observation 1.

There are problems with this if there is a trend because you'll get a
discontinuity when you join the two ends.

As you only have two years of data, and include season as a fixed effect
(main term in the model), you could assume that the residuals within
animal might be free of temporal dependence. You can check that as I
said by fitting the model and looking at the residuals within animal
over time.

Unless you are prepared to do some coding, the above discussion is moot
as vegan doesn't possess the code to use these restricted permutations
yet. If you want to do it yourself by hacking the anova.ccabyterm()
function, take a look at the vignette that comes with the permute
package for ideas on how to specify the permutation design you want and
then generate the permutations using shuffle() or shuffleSet() instead
of via permuted.index() that is used in vegan.

Basically, with:

mod <- cca(food ~ season*sex, data)
anova(mod, by = "term", strata = data$animal)

I guess you are allowing for a random intercept in animal only. For
animal specific season and year effects you'll need to include animal
and year in the model with appropriate interactions - we don't, as far
as I know, have ways of dealing with more complex dependencies in the
residuals via permutations built into vegan other than permuting within
animal, or maintaining temporal structure within animal (if you hack the
code yourself using functions from permute).

HTH

G

> thanks,
> René
> 
> 

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%



More information about the R-help mailing list