adonis in vegan: order of variables -- non nested with one degree of freedom for each

Question

I am trying to understand why the order of variables seems to effect p values for both of those variables for the adonis() (Permutational Multivariate Analysis of Variance Using Distance Matrices) test in vegan in R.

If I have a community data similarity matrix and two numeric predictors:

library(vegan)
data(dune)
data(dune.env)
dune.env2 <- dune.env
dune.env2$Moisture = as.numeric(dune.env$Moisture)

I can run adonis with the two variables trying to predict environmental similarity, with either varaible leading, like so:

adonis(dune ~ A1 + Moisture, data=dune.env2, permutations=9999)

Call:
adonis(formula = dune ~ A1 + Moisture, data = dune.env2, permutations = 9999) 

Permutation: free
Number of permutations: 9999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
A1         1    0.7230 0.72295  4.5132 0.16817  7e-04 ***
Moisture   1    0.8529 0.85292  5.3245 0.19840  2e-04 ***
Residuals 17    2.7232 0.16019         0.63344           
Total     19    4.2990                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

adonis(dune ~ Moisture + A1, data=dune.env2, permutations=9999)

Call:
adonis(formula = dune ~ Moisture + A1, data = dune.env2, permutations = 9999) 

Permutation: free
Number of permutations: 9999

Terms added sequentially (first to last)

          Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
Moisture   1    1.3782 1.37822  8.6039 0.32059 0.0001 ***
A1         1    0.1976 0.19765  1.2339 0.04598 0.2804    
Residuals 17    2.7232 0.16019         0.63344           
Total     19    4.2990                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

In the first case, both variables seem to have a statistically significatnt relationship to community structure. Meanwhile in the second case, only the Moisture variable appears to have a statistically significant relationship to community structure. In both cases, the residuals account for about 63% of the variance.

My working interpretation here is that in the first case, A1 ends up describing about 17% of the variance, and then Moisture predicts about 20% of what is left. Meanwhile in the second case Moisture describes 32% of the variance, and then there is not as much unexplained variance left for A1 to describe and so only Moisture ends up being statistically significant.

I find this surprising that adonis works this way, since with OLS regression, predictor variable order doesn't seem to matter this way.

summary(lm(log10(dune$Poaprat + 1) ~ dune.env2$A1 + dune.env2$Moisture))

Call:
lm(formula = log10(dune$Poaprat + 1) ~ dune.env2$A1 + dune.env2$Moisture)
>     
>     Residuals:
>          Min       1Q   Median       3Q      Max 
>     -0.32734 -0.11272 -0.02233  0.13877  0.42511 
>     
>     Coefficients:
>                        Estimate Std. Error t value Pr(>|t|)    
>     (Intercept)         0.97846    0.13337   7.337 1.16e-06 ***
>     dune.env2$A1       -0.04629    0.02677  -1.729   0.1019    
dune.env2$Moisture -0.12754    0.04431  -2.879   0.0104 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2239 on 17 degrees of freedom
Multiple R-squared:  0.5483,  Adjusted R-squared:  0.4951 
F-statistic: 10.32 on 2 and 17 DF,  p-value: 0.001165

 summary(lm(log10(dune$Poaprat + 1) ~ dune.env2$Moisture + dune.env2$A1 ))

Call:
lm(formula = log10(dune$Poaprat + 1) ~ dune.env2$Moisture + dune.env2$A1)
>     
>     Residuals:
>          Min       1Q   Median       3Q      Max 
>     -0.32734 -0.11272 -0.02233  0.13877  0.42511 
>     
>     Coefficients:
>                        Estimate Std. Error t value Pr(>|t|)    
>     (Intercept)         0.97846    0.13337   7.337 1.16e-06 ***
>     dune.env2$Moisture -0.12754    0.04431  -2.879   0.0104 *  
dune.env2$A1       -0.04629    0.02677  -1.729   0.1019    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2239 on 17 degrees of freedom
Multiple R-squared:  0.5483,  Adjusted R-squared:  0.4951 
F-statistic: 10.32 on 2 and 17 DF,  p-value: 0.001165

My question is similar to this one adonis in vegan: order of variables or use of strata. In that post, it was suggested that differences in degrees of freedom and nestedness accounted for this difference between models with different variable order. In my case, however, degrees of freedom for each variable equal one, and yet the pattern persists.

Could anyone help me in interpreting differences between adonis tests with different orders of predictor variables, even when each predictor variable has the same number of degrees of freedom? Thanks.

Answer 1

The output states:

Terms added sequentially (first to last)

So ordering is important. The issue is the same as the Type II (sequential) sums of squares of a linear model. While the partial effects don't change in your example, you still get the same order-sensitive behaviour as with adonis() if you compute the Type II sums of squares.

> anova(lm(log10(dune$Poaprat + 1) ~ A1 + Moisture, data = dune.env2))
Analysis of Variance Table

Response: log10(dune$Poaprat + 1)
          Df  Sum Sq Mean Sq F value   Pr(>F)   
A1         1 0.61907 0.61907 12.3462 0.002664 **
Moisture   1 0.41553 0.41553  8.2869 0.010422 * 
Residuals 17 0.85242 0.05014                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> anova(lm(log10(dune$Poaprat + 1) ~ Moisture + A1, data = dune.env2))
Analysis of Variance Table

Response: log10(dune$Poaprat + 1)
          Df  Sum Sq Mean Sq F value    Pr(>F)    
Moisture   1 0.88470 0.88470 17.6437 0.0006011 ***
A1         1 0.14990 0.14990  2.9894 0.1019254    
Residuals 17 0.85242 0.05014                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This occurs because A1 and Moisture are not uncorrelated; as there is some linear dependency, which ones goes into the model first determines how much variation is left to be explained by the second of the pair of covariates.

You can use dbrda(), for distance-based redundancy analysis, as it has better capabilities in its anova() method. Note that we also have adonis2() which provides the same marginal tests as I describe below.

ord <- dbrda(dune ~ A1 + Moisture, data = dune.env2, dist = 'bray')
anova(ord, by = 'margin')

which gives

> anova(ord, by = 'margin')
Permutation test for dbrda under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 999

Model: dbrda(formula = dune ~ A1 + Moisture, data = dune.env2, distance = "bray")
         Df SumOfSqs      F Pr(>F)    
A1        1  0.19765 1.2339  0.277    
Moisture  1  0.85292 5.3245  0.001 ***
Residual 17  2.72315                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

indicating that the marginal effect of Moisture after controlling for A1 is significant, but the reverse is not true.