Restricted permutations

Question

I have some multivariate data and want to investigate the effect of some environmental gradient. I want to use capscale but I don´t know how to deal with the permutation scheme. I have made up some artificial data, with 20 sites along a gradient ("env"):

######### create some species data along a gradient
df <- structure(list(site = 1:20, 
                     sp1 = c(7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp2 = c(1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp3 = c(0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp4 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L, 7L, 4L, 2L, 1L, 0L), 
                     sp5 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 2L, 4L, 7L, 8L), 
                     sp6 = c(0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp7 = c(0L, 0L, 0L, 0L, 0L, 0L, 2L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp8 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 3L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
                     sp9 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 4L, 1L, 0L, 0L)), 
                .Names = c("site", "sp1", "sp2", "sp3", "sp4", "sp5", "sp6", "sp7", "sp8", "sp9"), 
                class = "data.frame", row.names = c(NA, -20L))

# add some linear responses
df$sp6 <- round(seq(1, 8, 7/19), digits = 2) # linear response
df$sp7 <- round(seq(1, 4, 3/19), digits = 2) # no so strong linear response
df$sp9 <- round(seq(1, 6, 5/19), digits = 2)

# gradient
df$env <- 1:20

If I sampled only once I would do something like this:

# db-RDA sampled at one time
require(vegan)
mod <- capscale(df[ ,-c(1, 10)] ~ env, data = df, distance = "bray")
anova(mod, by = "terms", step = 999)  # assess the "significance" of contraining variable
plot(mod)

Now imagine I sampled the same data trice, but in three different months:

# now we replate exactly the same data 2 more times 
repdf <- rbind(df, df, df)  
repdf$time <- rep(1:3, each = nrow(df))
repdf$site <- factor(repdf$site)

If I would use unrestricted permutations, then this won't capture the repeated measures and the p-values would be to low.

I could restrict the permutations within each sites (using strata = site in vegan), but this destroys only the temporal effect and yields to a p of 1 (because every permuation is the same):

repmod <- capscale(repdf[ ,-c(1, 11, 12)] ~ env, data = repdf, distance = "bray") # db-RDA
anova(repmod, by = "terms", strata = repdf$site, step = 999)  

My question: How should I restrict the permutations assessing the effect of the gradient taking this temporal correlation into account? What permutationscheme should I use?

Some ideas: a) Permute the strata ( = sites), but not within the strata. This will destroy the env gradient, so the p-value is only determined by this.

b) Include "time" into the model (with interactions) and the run for time-effect a different permutation-scheme (permute within sites) than for env (permute sites, but not within sites).

I know about the permute-package and can incorporate it into permutest.cca, so my question is more theoretical.

Answer 1

If env is the variable you are interested in then we need to test the "effect" of this variable in terms of the variance explained by env. To do the test using a permutation we need to think about what is and isn't exchangeable under the Null hypothesis we are testing. In this case, the Null is no effect of env so under this the sites are freely exchangeable. Within site, the samples could be permuted but we should use a cyclic shift permutation rather than free permutation to preserve the autocorrelations. If the three time points are the same across all sites we might also wish to restrict the permutation so that the same permutation is used within all sites, thus preserving those correlations across site.

If time is a nuisance variable that you are not interested in, I would do a partial analysis, conditioning on time to remove it and look at the effect of env. To do this, include + Condition(time) in the model formula.

> mod <- capscale(repdf[, -c(1,11:12)] ~ env + Condition(time), data = repdf)
> mod
Call: capscale(formula = repdf[, -c(1, 11:12)] ~ env +
Condition(time), data = repdf)

              Inertia Proportion Rank
Total         38.1148     1.0000     
Conditional    0.0000     0.0000    1
Constrained   16.5286     0.4337    1
Unconstrained 21.5862     0.5663    9
Inertia is mean squared Euclidean distance 

Eigenvalues for constrained axes:
 CAP1 
16.53 

Eigenvalues for unconstrained axes:
     MDS1      MDS2      MDS3      MDS4      MDS5      MDS6 
8.890e+00 7.312e+00 4.611e+00 4.793e-01 1.560e-01 1.370e-01 
     MDS7      MDS8      MDS9 
7.542e-06 5.719e-06 3.341e-06

To permute strata (here site), you can use permute. To shuffle just the strata as blocks of three samples, a suitable control object would be:

ctrl <- with(repdf, permControl(strata = site,
                                within = Within(type = "none"),
                                blocks = Blocks(type = "free")))

If you wanted to permute within site we would use this

ctrl <- with(repdf, permControl(strata = site,
                                within = Within(type = "series",
                                                constant = TRUE),
                                blocks = Blocks(type = "free")))

The constant = TRUE forces the same random cyclic shift to be applied across each site.

Even controlling for time (by conditioning on it in partial analysis) this only removes a linear effect of time so there could be some remaining autocorrelation in the residuals that are then used as the response data in the model testing the effect of env, hence the reason for using cyclic permutations even though we controlled for time in the model.

Something to be aware of is that when permuting strata, you have to have a balanced design, as here, with the same number of samples within each level of strata. The code that generates permutations from control objects does not enforce this as it is designed to be very general purpose. There is a function permCheck() that exists to check a given control object and balance is one of the checks performed. When we interface permute with vegan we will enforce balance using permCheck() but I don't think permCheck() works correctly in the current version of permute.

Of course, as we (well, I) haven't yet interfaced permute and vegan yet, you'll need to hack the anova() method to generate the relevant permutations.

---

[Apologies for not responding sooner despite being alerted to this by @chlalanne on Twitter - my reply has nothing to do with you awarding a bounty, I had a poorly 11-month-old asleep on my chest most of the past few days that curtailed my typing ;-)]