Something is off with adonis and adonis2 output with strata/perm constraints

[CarlyMuletzWolz]

Hi,

I am going back through code that I ran some months ago, and getting wildly different p-values, but all the same summary stats.

I ran the code prior to strata deprecation in adonis2 and wondering if something is going on with that code that is miscalculating. Visually, it looks like my prior results were correct and something is going on with new p-value calculations that is incorrect.

Here is the old model set-up

b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, strata = Animal.Name, data=df.rep)

Output pasted from saved excel file

 | Df | SumOfSqs | R2 | F | Pr(>F)

SuperOrder | 2 | 2.23203697 | 0.1294166 | 2.64512209 | 0.001
Lactation.Stage | 1 | 0.28383607 | 0.01645721 | 0.67273174 | 0.998
DietType | 1 | 0.93520187 | 0.0542243 | 2.21656106 | 0.001
Habitat | 1 | 0.73570578 | 0.04265724 | 1.74372704 | 0.005
SuperOrder:Lactation.Stage | 2 | 0.53151843 | 0.03081818 | 0.62988704 | 1
Lactation.Stage:DietType | 1 | 0.29306001 | 0.01699203 | 0.69459379 | 0.977
Residual | 29 | 12.2355547 | 0.70943444 | NA | NA
Total | 37 | 17.2469138 | 1 | NA | NA

Here is the new model set-up:

perm <- how(nperm = 1000)
setBlocks(perm) <- with(df.rep, Animal.Name)

b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, data=df.rep, permutations = perm)
b2

And output

> b2
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(df.rep, Animal.Name) 
Permutation: free
Number of permutations: 1000

adonis2(formula = brayRepProp ~ SuperOrder * Lactation.Stage + DietType * Lactation.Stage + Habitat, data = df.rep, permutations = perm)
                           Df SumOfSqs      R2      F Pr(>F)
SuperOrder                  2   2.2320 0.12942 2.6451 0.8661
Lactation.Stage             1   0.2838 0.01646 0.6727 0.7742
DietType                    1   0.9352 0.05422 2.2166 0.8661
Habitat                     1   0.7357 0.04266 1.7437 0.8661
SuperOrder:Lactation.Stage  2   0.5315 0.03082 0.6299 0.7842
Lactation.Stage:DietType    1   0.2931 0.01699 0.6946 0.7872
Residual                   29  12.2356 0.70943              
Total                      37  17.2469 1.00000      

Also, when I use adonis and strata it gives similar as new adonis2

b2 <- adonis(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, strata = df.rep$Animal.Name, data=df.rep)

Call:
adonis(formula = brayRepProp ~ SuperOrder * Lactation.Stage +      DietType * Lactation.Stage + Habitat, data = df.rep, strata = df.rep$Animal.Name) 

Blocks:  strata 
Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

                           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
SuperOrder                  2    2.2320 1.11602 2.64512 0.12942  0.861
Lactation.Stage             1    0.2838 0.28384 0.67273 0.01646  0.753
DietType                    1    0.9352 0.93520 2.21656 0.05422  0.861
Habitat                     1    0.7357 0.73571 1.74373 0.04266  0.861
SuperOrder:Lactation.Stage  2    0.5315 0.26576 0.62989 0.03082  0.789
Lactation.Stage:DietType    1    0.2931 0.29306 0.69459 0.01699  0.767
Residuals                  29   12.2356 0.42192         0.70943       
Total                      37   17.2469                 1.00000       

Notice all same sums of squares, f-values, dfs, but different p-values for all three varieties

When I just account for animal.name as a covariate, I get similar to what I would expect based on plotting the data

b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat + Animal.Name, data=df.rep)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = brayRepProp ~ SuperOrder * Lactation.Stage + DietType * Lactation.Stage + Habitat + Animal.Name, data = df.rep)
                           Df SumOfSqs      R2      F Pr(>F)    
SuperOrder                  2   2.2320 0.12942 3.6598  0.001 ***
Lactation.Stage             1   0.2838 0.01646 0.9308  0.607    
DietType                    1   0.9352 0.05422 3.0668  0.001 ***
Habitat                     1   0.7357 0.04266 2.4126  0.001 ***
Animal.Name                14   7.6614 0.44422 1.7946  0.001 ***
SuperOrder:Lactation.Stage  2   0.5315 0.03082 0.8715  0.872    
Lactation.Stage:DietType    1   0.2931 0.01699 0.9610  0.519    
Residual                   15   4.5741 0.26521                  
Total                      37  17.2469 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Here is one figure showing how super order and diet type indeed look significant

But, lactation stage does not

Thoughts on what is going on here?

[jarioksa]

Could you please concisely specify what is the problem: I fail to find it.

[gavinsimpson]

@jarioksa the issue seems to be that with strata = fct they got one set of results and since deprecation of strata adding the equivalent code to do a restricted permutation setting blocks to be fct they get some very different results.

@CarlyMuletzWolz Can you check with the simpler (in terms of the underlying code):

perms <- with(df.rep, how(nperm = 1000, blocks = Animal.Name))

as this could signal an issue with whether the blocks are being set on the permutations design post hoc.

[gavinsimpson]

@CarlyMuletzWolz I'll add, what I'm seeing with the example in ?adonis2 doesn't suggest there is any issue in permuting within blocks

r$> shuffleSet(nset = 10, control = perm, n = nrow(dat))                                                
No. of Permutations: 10
No. of Samples: 12 (Nested in: blocks; Randomised)
Restricted by Blocks: with(dat, field) (3 blocks)

    1 2 3 4 5 6 7 8  9 10 11 12
p1  1 4 2 3 7 5 8 6  9 10 12 11
p2  4 1 3 2 8 5 6 7 11 10 12  9
p3  4 2 3 1 5 6 8 7 12  9 11 10
p4  3 1 2 4 6 5 8 7 12  9 10 11
p5  3 2 4 1 8 5 7 6  9 10 11 12
p6  3 2 1 4 6 7 8 5 12  9 10 11
p7  2 1 3 4 7 5 8 6 12 11  9 10
p8  1 3 4 2 7 6 5 8  9 10 11 12
p9  1 4 3 2 7 5 6 8 12 11  9 10
p10 1 3 4 2 5 7 8 6 11  9 10 12

where there are 12 samples, 3 blocks, 4 samples each. You can see that the numbers 1:4 only ever occur in the first 4 positions of the permutation, and same for the other blocks of 4 samples.

You should probably confirm that you really did store the correct output from the previous run by reinstalling an older version of vegan with strata working and confirm that the results you obtain with that version match those you saved in the Excel file and are different to the ones you are getting now

[jarioksa]

It indeed is possible in adonis2 that there is a confusion with strata: these were never implemented nor documented to work with adonis2, but they were silently ignored in vegan pre 2.5-7, and non-stratified result was returned. If this was the point. In 2.5-7 we issue an error if strata are used with adonis2.

[CarlyMuletzWolz]

Hi all,

Indeed verified if I go back to version 2.5.6 I get completely different p-values. Again, note that all of the other summary statistics are the same, so something is off when it gets to p-value calculations.

##vegan 2.5.6 

b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, strata = Animal.Name, data=df.rep, REML = F)
> b2
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = brayRepProp ~ SuperOrder * Lactation.Stage + DietType * Lactation.Stage + Habitat, data = df.rep, strata = Animal.Name, REML = F)
                           Df SumOfSqs      R2      F Pr(>F)    
SuperOrder                  2   2.0681 0.10754 2.4509  0.001 ***
Lactation.Stage             1   0.2755 0.01432 0.6529  0.996    
DietType                    3   2.6621 0.13842 2.1032  0.001 ***
Habitat                     1   0.7357 0.03826 1.7437  0.003 ** 
SuperOrder:Lactation.Stage  2   0.5191 0.02699 0.6151  1.000    
Lactation.Stage:DietType    3   0.7353 0.03823 0.5809  1.000    
Residual                   29  12.2356 0.63623                  
Total                      41  19.2313 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> packageVersion("vegan")
[1] ‘2.5.6


## vegan2.5.7
> perm <- how(nperm = 1000)
> setBlocks(perm) <- with(df.rep, Animal.Name)

> b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, data=df.rep, permutations = perm)
> b2
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  with(df.rep, Animal.Name) 
Permutation: free
Number of permutations: 1000

adonis2(formula = brayRepProp ~ SuperOrder * Lactation.Stage + DietType * Lactation.Stage + Habitat, data = df.rep, permutations = perm)
                           Df SumOfSqs      R2      F Pr(>F)
SuperOrder                  2   2.0681 0.10754 2.4509 0.8452
Lactation.Stage             1   0.2755 0.01432 0.6529 0.7872
DietType                    3   2.6621 0.13842 2.1032 0.8452
Habitat                     1   0.7357 0.03826 1.7437 0.8452
SuperOrder:Lactation.Stage  2   0.5191 0.02699 0.6151 0.7592
Lactation.Stage:DietType    3   0.7353 0.03823 0.5809 0.8402
Residual                   29  12.2356 0.63623              
Total                      41  19.2313 1.00000              
> packageVersion("vegan")
[1] ‘2.5.7’

[CarlyMuletzWolz]

Hi Gavin,

And, it is the same with using the simpler line of code

> perm <- with(df.rep, how(nperm = 1000, blocks = Animal.Name))
> b2 <- adonis2(brayRepProp ~ SuperOrder*Lactation.Stage + DietType*Lactation.Stage + Habitat, data=df.rep, permutations = perm)
> b2
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  Animal.Name 
Permutation: free
Number of permutations: 1000

adonis2(formula = brayRepProp ~ SuperOrder * Lactation.Stage + DietType * Lactation.Stage + Habitat, data = df.rep, permutations = perm)
                           Df SumOfSqs      R2      F Pr(>F)
SuperOrder                  2   2.0681 0.10754 2.4509 0.8661
Lactation.Stage             1   0.2755 0.01432 0.6529 0.8012
DietType                    3   2.6621 0.13842 2.1032 0.8661
Habitat                     1   0.7357 0.03826 1.7437 0.8661
SuperOrder:Lactation.Stage  2   0.5191 0.02699 0.6151 0.7882
Lactation.Stage:DietType    3   0.7353 0.03823 0.5809 0.8202
Residual                   29  12.2356 0.63623              
Total                      41  19.2313 1.00000 

[CarlyMuletzWolz]

Ok, and with the example data provided. Completely different p-values between using strata (2.5.6) and the new format (2.5.7)

dat <- expand.grid(rep=gl(2,1), NO3=factor(c(0,10)),field=gl(3,1) )
dat

Agropyron <- with(dat, as.numeric(field) + as.numeric(NO3)+2) +rnorm(12)/2
Schizachyrium <- with(dat, as.numeric(field) - as.numeric(NO3)+2) +rnorm(12)/2
total <- Agropyron + Schizachyrium


Y <- data.frame(Agropyron, Schizachyrium)

## For vegan 2.5.6, with REML = F and T

adonis2(Y ~ NO3, strata = field, data=dat)

> adonis2(Y ~ NO3, strata = field, data=dat, REML = F)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = Y ~ NO3, data = dat, strata = field, REML = F)
         Df SumOfSqs      R2      F Pr(>F)
NO3       1 0.022236 0.12014 1.3654  0.237
Residual 10 0.162849 0.87986              
Total    11 0.185085 1.00000              
> adonis2(Y ~ NO3, strata = field, data=dat)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = Y ~ NO3, data = dat, strata = field)
         Df SumOfSqs      R2      F Pr(>F)
NO3       1 0.022236 0.12014 1.3654  0.246
Residual 10 0.162849 0.87986              
Total    11 0.185085 1.00000 


## For vegan 2.5.7

perm <- with(dat, how(nperm = 1000, blocks = field))

adonis2(Y ~ NO3, data = dat, permutations = perm)

> adonis2(Y ~ NO3, data = dat, permutations = perm)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  field 
Permutation: free
Number of permutations: 1000

adonis2(formula = Y ~ NO3, data = dat, permutations = perm)
         Df SumOfSqs      R2      F  Pr(>F)  
NO3       1 0.022236 0.12014 1.3654 0.04995 *
Residual 10 0.162849 0.87986                 
Total    11 0.185085 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


## Why such different results?

## And, if I run strata in 2.5.7, I get something similar to 2.5.6 results

adonis2(Y ~ NO3, strata = field, data=dat)

> adonis2(Y ~ NO3, strata = field, data=dat)
Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = Y ~ NO3, data = dat, strata = field)
         Df SumOfSqs      R2      F Pr(>F)
NO3       1 0.022236 0.12014 1.3654  0.241
Residual 10 0.162849 0.87986              
Total    11 0.185085 1.00000 

[jarioksa]

I think I have already answered to this long and winding question. Here a recap:

1. adonis2 has no argument strata and never had.
2. In vegan 2.5-6 and earlier, adonis2 silently ignored non-existing strata argument and performed non-stratified analysis.
3. In vegan 2.5-7 adonis2 stops with error if you try to use non-existing strata argument. Your last example is not from vegan 2.5-7 although it is so labelled.
4. The lack of strata in adonis2 is documented, but not emphasized: we now issue an error to make this clear. However, output also tells that there are no blocks of permutation: you only get Permutation: free in the output without reference to strata or blocks. I apologize for not catching the user error before 2.5-7.
5. The Examples section of the manual tells how to implement strata as blocks: this should be followed if you want to have strata in adonis2. The same procedure must be followed both in vegan 2.5-6 and 2.5-7: adonis2 has never used strata.
6. Minor issue: there is no argument REML in adonis or adonis2 and it has no effect.

[CarlyMuletzWolz]

Hi Jari,

I was trying to provide information that would potentially be pertinent to the issue as I did not know where the issue stemmed from. Thank you for looking into this and providing the re-cap. It is now clear - strata was doing nothing in 2.5-6 version and that is why the p-values are so different. Thank you for now issuing a warning.