User Tools

Site Tools


en:varpart_examples

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 R2 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.

en/varpart_examples.txt · Last modified: 2019/02/26 23:37 by David Zelený