Section: Ordination analysis

## Variation partitioning (constrained ordination)

### Example 1: Use of varpart function (using Barley field weed community dataset)

Use data from Effect of fertilizer on composition of weed community in barley fields. Authors focused on weed communities in barley fields; weeds are plants accompanying the crop (here barley), not wanted by farmers (because they feel weed is competing with crop), but often containing rare and endangered species). They asked how the species composition of weeds changes after application of fertiliser on the barley field. The problem is that fertiliser may influence the weed community directly, but also indirectly via the increased cover of barley (indeed, barley grows better if fertilised). This example aims to separate variation in species composition of weed community caused by the amount of fertiliser (`dose`

) and by estimated `cover`

of barley at the situation when both variables are correlated (since fertiliser causes the barley to grow better and have higher cover).

First, import the data:

fertil.spe <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/fertil.spe.txt', row.names = 1) fertil.env <- read.delim ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/fertil.env.txt', row.names = 1)

Species data in `fertil.spe`

are represented by the estimated cover of weed community species in a plot, and it is not necessary to transform them. One may first apply DCA ordination on species data to see the length of the first DCA axis and hence seek the recommendation whether the linear or unimodal method should be used (the length is 3.8 S.D., lying in a grey zone between 3 and 4 S.D.); we will use RDA here.

First, let's see how would the variation partitioning looks like if we do it step-by-step manually (without using `varpart`

function). We will partition variation explained by variables `dose`

and `cover`

from `fertil.env`

. For this, we need to define the global model (to see how much variation is explained by both of them), and also models with simple (marginal) effect of each explanatory variable:

# fractions [a+b+c]: rda.all <- rda (fertil.spe ~ dose + cover, data = fertil.env) # fractions [a+b]: rda.dose <- rda (fertil.spe ~ dose, data = fertil.env) # fractions [b+c]: rda.cover <- rda (fertil.spe ~ cover, data = fertil.env)

Now we can calculate variations explained by individual fractions (using `RsquareAdj`

function):

## fraction [a+b+c] RsquareAdj (rda.all) # $`r.squared` # [1] 0.1430381 # # $adj.r.squared # [1] 0.1286354 abc <- RsquareAdj (rda.all)$adj.r.squared ## fraction [a+b]: RsquareAdj (rda.dose) # $`r.squared` # [1] 0.06889364 # # $adj.r.squared # [1] 0.06113442 ab <- RsquareAdj (rda.dose)$adj.r.squared ## fraction [b+c]: RsquareAdj (rda.cover) # $`r.squared` # [1] 0.09628412 # # $adj.r.squared # [1] 0.08875315 bc <- RsquareAdj (rda.cover)$adj.r.squared

In this way, we calculated adjusted R^{2} of fractions [a+b+c], [a+b] and [b+c]. To calculate values for individual fractions [a], [b] and [c], we need to do a simple subtractions:

b <- ab + bc - abc # 0.02125216 a <- ab - b # 0.03988225 c <- bc - b # 0.06750099

After recalculating into percentage, we get:

- [a] = 4.0%
- [b] = 2.1%
- [c] = 6.8%

Conditional effect of dose (fraction [a] = 4.0%) is slightly lower than conditional effect of cover ([c] = 6.8%), and shared variance is not too high ([b] = 2.1%). Seems that weed species composition is more affected by light conditions modified by increasing cover of barley, than by fertilizer itself.

Now, let's see the same, using function `varpart`

(note that by default, `varpart`

function calculates RDA on the matrices; if you want to use CCA, you need to use argument `chisquare = TRUE`

).

varp <- varpart (fertil.spe, ~ dose, ~ cover, data = fertil.env) varp

Partition of variance in RDA Call: varpart(Y = fertil.spe, X = ~dose, ~cover, data = fertil.env) Explanatory tables: X1: ~dose X2: ~cover No. of explanatory tables: 2 Total variation (SS): 2581.2 Variance: 21.332 No. of observations: 122 Partition table: Df R.squared Adj.R.squared Testable [a+b] = X1 1 0.06889 0.06113 TRUE [b+c] = X2 1 0.09628 0.08875 TRUE [a+b+c] = X1+X2 2 0.14304 0.12864 TRUE Individual fractions [a] = X1|X2 1 0.03988 TRUE [b] 0 0.02125 FALSE [c] = X2|X1 1 0.06750 TRUE [d] = Residuals 0.87136 FALSE --- Use function ‘rda’ to test significance of fractions of interest

The results for individual fractions are identical with our results above. `varpart`

reports also simple (marginal) effect (both not-adjusted and adjusted variation) of individual predictors, i.e. fractions [a+b] and [b+c].

We may plot the results into Venn's diagram (argument `digits`

influences number of decimal digits shown in the diagram, `Xnames`

the displayed names of variables, and `bg`

background color of the diagram fractions; see `?varpart`

for details):

plot (varp, digits = 2, Xnames = c('dose', 'cover'), bg = c('navy', 'tomato'))

Now, when we know both simple and conditional effect of each variables, we may want to know whether these variances are significant, and hence worth of interpreting. Results from `varpart`

contain the column `testable`

with logical values indicating whether given fraction is testable or not. To test each of them, we will need the models defined above, and the function `anova`

, which (if applied on single object resulting from `rda`

or `cca`

method, returns Monte Carlo permutation test of the predictor effect). For this, we need to first define also partial ordination models with one variable as exlanatory and the other as covariable (`Condition`

):

# fraction [a]: rda.dose.cover <- rda (fertil.spe ~ dose + Condition (cover), data = fertil.env) # fraction [c]: rda.cover.dose <- rda (fertil.spe ~ cover + Condition (dose), data = fertil.env)

Then, to test individual fractions, we use:

## The global model (fractions [a+b+c]): anova (rda.all) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(formula = fertil.spe ~ dose + cover, data = fertil.env) # Df Variance F Pr(>F) # Model 2 3.0513 9.9313 0.001 *** # Residual 119 18.2808 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## Simple (marginal) effect of dose (fraction [a+b]): anova (rda.dose) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(formula = fertil.spe ~ dose, data = fertil.env) # Df Variance F Pr(>F) # Model 1 1.4696 8.8789 0.001 *** # Residual 120 19.8624 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## ## Simple (marginal) effect of cover (fraction [b+c]): anova (rda.cover) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(formula = fertil.spe ~ cover, data = fertil.env) # Df Variance F Pr(>F) # Model 1 2.0539 12.785 0.001 *** # Residual 120 19.2781 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## Conditional (partial) effect of dose (fraction [a]): anova (rda.dose.cover) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(formula = fertil.spe ~ dose + Condition(cover), data = fertil.env) # Df Variance F Pr(>F) # Model 1 0.9974 6.4924 0.001 *** # Residual 119 18.2808 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## Conditional (partial) effect of cover (fraction [c]): anova (rda.cover.dose) # Permutation test for rda under reduced model # Permutation: free # Number of permutations: 999 # # Model: rda(formula = fertil.spe ~ cover + Condition(dose), data = fertil.env) # Df Variance F Pr(>F) # Model 1 1.5817 10.296 0.001 *** # Residual 119 18.2808 # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

From these results, you may see that all simple (marginal) and conditional (partial) effects of both predictors are significant at P < 0.001.