Introduction

The purpose of this chapter is to illustrate how to control permutations.  This is important primarily so that all aspects of an analysis – including the P-value – are replicable.  We’ve already seen multiple examples of analyses in which the test statistic is identical but the statistical significance varies.  Here, we want to consider ways to control which set of permutations is used.  Doing so means that the statistical significance will also be unchanged from one run to another, even if those runs are at different times or even on different machines.

These notes are illustrated with analyses using vegan::adonis2().  However, the permutation strategies that are developed here can be applied to any of the functions from the vegan package that use the same permutations argument.

RRPP is coded differently, and has a built-in ability to control permutations.  In RRPP::lm.rrpp(), the set of permutations is identified by the number of permutations requested.  As long as you request the same number of permutations, the statistical significance will be the same from run to run.  If you change the number of permutations, a different set of permutations is used and the statistical significance may differ.

Grazing Example

We will illustrate the techniques in this section using our oak dataset.

source("scripts/load.oak.data.R")

The oak plant community dataset contains factors indicating whether stands were grazed in the past and/or grazed at the time the data were collected (GrazPast and GrazCurr, respectively).  As we’ve done before, we’ll combine these factors together into a new Grazing factor (capitalized):

Grazing <- factor(with(Oak, paste(GrazPast, GrazCurr, sep = "_")))

summary(Grazing)

  No_No  Yes_No Yes_Yes 
     24       6      17

How to Control Permutations

Permutations are conducted by drawing on a random number generator within R (we won’t worry about how it uses this random number generator).

It is sometimes helpful to know exactly which permutations(s) you used so that, for example, you can exactly replicate an analysis.  Here are three ways to do so.

set.seed()

The sample() function was introduced in the chapter about permutations tests.  This function simply rearranges the data – it doesn’t change the data itself.  For example, here is a rearrangement of the numbers from 1 to 10:

sample(10)

6  4 10  1  9  7  5  2  8  3

If we re-run it, we get a different rearrangement:

sample(10)

2 10  3  8  6  4  1  7  9  5

To exactly replicate a permutation test, we need to use the same set of permutations and therefore have to start at the same value in the random number generator.  The set.seed() function allows you to do so.  If you set this, run a permutation test, re-set it, and re-run the permutation test, you will get the same permutation(s) each time and thus exactly the same results.

set.seed(42); sample(10)

10  9  3  6  4  8  5  1  2  7

set.seed(42); sample(10) #identical to previous

10  9  3  6  4  8  5  1  2  7

set.seed(40); sample(10) #different

7  8  6  1  2  3 10  9  4  5

set.seed(NULL) #reset random number generator
sample(10)

8  1  6  7  4 10  3  9  5  2

Each time you call the random number generator you ‘use’ a different set of numbers.  This means that the set.seed() function should be reset before each use.  For example, setting it once and then generating two sets results in two different sets:

set.seed(42); sample(10); sample(10)

 [1]  1  5 10  8  2  4  6  9  7  3
 [1]  8  7  4  1  5 10  2  6  9  3

Resetting it before each set ensures that they’re the same:

set.seed(42); sample(10); set.seed(42); sample(10)

 [1]  1  5 10  8  2  4  6  9  7  3
 [1]  1  5 10  8  2  4  6  9  7  3

In older versions of R (pre-3.6.0), this function gave identical results regardless of platform and in different sessions.  Newer versions of R use a different random number generator by default and therefore do not necessarily give identical results unless you change the type of random number generator that is being used.  An explanation of this issue is here.  You can change the random number generator for an individual action or for an entire session.

Changing one action:

set.seed(40, sample.kind = "Rounding"); sample(10)

7  8  6  1  2  3 10  9  4  5

This value should match among all of our machines and runs.

Changing the random number generator for the rest of an R session:

RNGkind(sample.kind = "Rounding")

This function would need to be run once per session, before using set.seed().

Resetting to default random number generator:

RNGkind(sample.kind = "default")

sample() and permute::shuffle()

You can create a permutation using the sample() function, which we’ve seen earlier, or the shuffle() function.  The shuffle() function has the following usage:

shuffle(
n,
control = how()
)

The main arguments are:

• n – how many permuted values to return.  Usually set to the number of samples in the object.
• control – how the sample units are to be permuted.  See chapter about restricting permutations for more information about the how() function.

The result is a permutation of the integers from 1:n, just like we saw with sample() above.  For example, here is a permutation of the sample units in our oak dataset:

perm1 <- shuffle(n = length(Grazing))

 [1] 44 36 15 41 35 21 28 26 32  4 47 23 25 17  6  1 46 24  5  8 40  3 42
[24] 37 29  7 31 39 13 20 12 34  9 14  2 45 38 10 43 22 11 18 30 27 16 19
[47] 33

Note how the integers from 1 to 47 have been permuted.  Each integer refers to a position in the original object.  For example, the first value (shown in bold) is now in the 16th position.

This permutation can be used to index other objects as occurs in a permutation test.  For example, let’s use this permutation to index the species richness of the stands.  Here’s the original species richness data:

Oak_explan[ , "SppRich"]

 [1] 32 41 35 33 32 51 37 30 29 38 35 33 38 45 41 34 23 22 22 32 29 18 51
[24] 20 32 36 40 26 34 35 32 18 46 16 38 49 40 33 41 28 38 33 61 44 27 38
[47] 33

It can be confusing to distinguish sets of numbers.  These aren’t the integers that we saw above – these values don’t get below 16 and get as high as 61.  The first stand has 32 species (shown in bold).

Here’s the permuted species richness data:

Oak_explan[perm1, "SppRich"]

 [1] 44 49 41 38 38 29 26 36 18 33 33 51 32 23 51 32 38 20 32 30 28 35 33
[24] 40 34 37 32 41 38 32 33 16 29 45 41 27 33 38 61 18 35 22 35 40 34 22
[47] 46

By comparing these sets of data, you can verify that the order of the species richness values has been adjusted based on the ordering in perm1.  For example, the integer 1 is in the 16th position in perm1 (shown in bold).   In the permuted species richness data, the richness from stand 1 is now in the 16th position.

Can you calculate the average richness in each grazing treatment, using both the real data and the ordering in perm1?

Start with the real data:

Oak_explan |>
mutate(Grazing = Grazing) |>
  group_by(Grazing) |>
reframe(SppRich = mean(SppRich))

# A tibble: 3 × 2
  Grazing SppRich
  <fct>     <dbl>
1 No_No      30.4
2 Yes_No     36.2
3 Yes_Yes    39.6

These are the average number of species in each grazing treatment – on average, stands in the Yes_Yes treatment contain the most species and those in the No_No treatment contain the least species.

We could use a set of permutations to assess statistical significance, as we’ve seen throughout this section of the course.  Here are the average values based on one permutation:

Oak_explan |>
mutate(Grazing = Grazing[perm1]) |>
group_by(Grazing) |>
reframe(SppRich = mean(SppRich))

# A tibble: 3 × 2
  Grazing SppRich
  <fct>     <dbl>
1 No_No      35.1
2 Yes_No     35.8
3 Yes_Yes    33

As expected, the average richness values are all very similar to one another once the data are permuted.

We could have permuted either Grazing or SppRich but we would not have wanted to do both … do you see why?

permute::shuffleSet()

Since shuffle() creates one permutation, you have to call it each time you want to generate a new permutation.  The shuffleSet() function works just like shuffle(), except that it includes an additional argument (nset) that allows you to specify how many permutations to create at once.

perms <- shuffleSet(n = length(Grazing), nset = 5); perms

Note the dimensionality of this object: one row per permutation and one column per sample unit.

This set of permutations could be used in a for() loop or a matrix algebra calculation to calculate the results of all permutations.  See the ‘Restricting Permutations‘ chapter for an example.

For even more control over this set of permutations, it can be proceeded by set.seed().  This would ensure that the entire set of permutations is replicable.