Constrained ordination is a method to explore effects of environmental fators to the assemblage ordination by an extension of multivariate regression. Here, we use the same data sets of the Unconstrained Ordination; dat and edat.

Formatting environmental variables

We use environmental factors as explanation variables. We have to check assumptions as in ordinal multivariage analyses. Three points are examined in the example below.

• Normality
• Multiple-colinearlity
• Standardization

Although multivariage regression is somowhat robust to violation of normality, it should be better to transform variables of excessively skewed distribution. Check distribution of variables by a function hist().

hist(edat$MdG)

A variable MdG is the median grain size of the sediment measured by mm scale. Because it looks skewed, it is transformed into phi scale.

edat$MdPh <- -log(edat$MdG, 2)
hist(edat$MdPh)

It's still skewed, but looks better. Other tree variables are transformed into their square-rooted values.

edat$sqPM <- sqrt(edat$PM)
edat$sqTN <- sqrt(edat$TN)
edat$sqTOC <- sqrt(edat$TOC)

Still skewed, but move to the next.

Multiple-colinearlity between variables are examined by variance inflation factor (VIF). When VIF > 10, problems of multiple-colinearlity become severe. So, remove variables of large VIF, one by one. First, prepare a new data set and check pair-wise correlations by functions cor() and pairs(). Following example uses pairs.panels function in psych package.

edat2 <- edat[,c(1:7,12:15)]
psych::pairs.panels(edat2)

vif <- diag(solve(cor(edat2)))
vif

##       PSU        PH      ORP3      ORP5      ORP8        IL      Chla 
##  4.638988  4.147434 12.680515 29.958645 19.160637  6.747144 22.484479 
##      MdPh      sqPM      sqTN     sqTOC 
##  6.649999  7.595803 17.796963  9.618832

VIF exceeds 10. Three variables of sediment redox potentials (ORP3, ORP5, and ORP8) seem to have problems. Now, we remove ORP5.

vif <- diag(solve(cor(edat2[,-4])))
vif

##       PSU        PH      ORP3      ORP8        IL      Chla      MdPh 
##  3.792806  3.433859  9.190892 12.508433  6.046272 20.499752  5.710017 
##      sqPM      sqTN     sqTOC 
##  7.401397 16.304317  9.043974

Check the VIF and remove variables, consecutively.

vif <- diag(solve(cor(edat2[,c(-4, -7, -11)])))
vif

##      PSU       PH     ORP3     ORP8       IL     MdPh     sqPM     sqTN 
## 2.545664 1.474089 2.594832 2.617503 2.188238 1.948444 2.347553 1.865053

Now, VIF is 2.6175032: well below to the reference value of 10. Name the new variable as edat2 and standardize it.

edat2 <- edat2[,c(-4, -7, -11)]
edat2 <- data.frame(scale(edat2))
colMeans(edat2)

##           PSU            PH          ORP3          ORP8            IL 
##  2.016905e-16 -2.035409e-15  9.066821e-17  5.551115e-17 -5.458597e-17 
##          MdPh          sqPM          sqTN 
##  4.070818e-17  1.850372e-17 -6.476301e-17

apply(edat2, 2, sd)

##  PSU   PH ORP3 ORP8   IL MdPh sqPM sqTN 
##    1    1    1    1    1    1    1    1

Now, each variables has its mean = 0 and sd = 1.

Constrained ordination by dbRDA

Distance-based Redunduncy Analysis (dbRDA) performs multiple regression of assemblage by environmental factors as explanation variables. In the analysis, ordination by PCoA is dependent variables. Because it is an extension of multiple regression, it is possible to forcast a new data set, but linear relationships are assumed.

To perform dbRDA, we can use functions capscale() or dbrda() in vegan package. Here, we use an example of dbrda().

res1.rda <- dbrda(dat~., data=edat2, distance="horn")
res1.rda

## Call: dbrda(formula = dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh +
## sqPM + sqTN, data = edat2, distance = "horn")
## 
##               Inertia Proportion Rank RealDims
## Total          4.6471     1.0000              
## Constrained    3.1141     0.6701    8        7
## Unconstrained  1.5329     0.3299    6        5
## Inertia is squared Horn distance 
## 
## Eigenvalues for constrained axes:
##  dbRDA1  dbRDA2  dbRDA3  dbRDA4  dbRDA5  dbRDA6  dbRDA7 idbRDA1 
##  1.1001  0.7293  0.5035  0.4615  0.2383  0.1006  0.0233 -0.0425 
## 
## Eigenvalues for unconstrained axes:
##    MDS1    MDS2    MDS3    MDS4    MDS5   iMDS1 
##  0.8778  0.3836  0.1772  0.1471  0.0103 -0.0631

plot(res1.rda)

Check valuse of "Total Inertia" and "Constrained Inertia". Also check eigenvalues > 1.

anova.cca(res1.rda)

## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh + sqPM + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)  
## Model     8   3.1141 1.5236  0.066 .
## Residual  6   1.5329                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Significance test of dbRDA shows that this is not significant, unfortunately.

Variable selection is performed by a function ordistep(). There are two options: directions of increasing variables or decreasing variables. Here, we show an example of decreasing variables from the full model. A function anova.cca() examines significance of the dbRDA.

res0.rda <- dbrda(dat~1, data=edat2, distance="horn")
res2.rda <- ordistep(res1.rda, scope=formula(res0.rda), direction="both")

## 
## Start: dat ~ PSU + PH + ORP3 + ORP8 + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)  
## - ORP3  1 22.216 0.3468  0.905  
## - ORP8  1 22.365 0.4105  0.865  
## - IL    1 24.001 1.1492  0.370  
## - sqPM  1 24.180 1.2346  0.285  
## - sqTN  1 25.044 1.6637  0.140  
## - MdPh  1 25.209 1.7486  0.080 .
## - PSU   1 25.427 1.8619  0.070 .
## - PH    1 28.043 3.3598  0.020 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + ORP8 + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - ORP8  1 20.833 0.2943  0.930   
## - IL    1 22.617 1.2153  0.305   
## - sqPM  1 22.752 1.2898  0.250   
## - MdPh  1 23.915 1.9581  0.085 . 
## - sqTN  1 23.591 1.7665  0.075 . 
## - PSU   1 23.934 1.9695  0.060 . 
## - PH    1 27.119 4.0906  0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + IL + MdPh + sqPM + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - sqPM  1 21.302 1.4309  0.175   
## - IL    1 21.874 1.7980  0.115   
## - sqTN  1 22.046 1.9108  0.090 . 
## - PSU   1 22.390 2.1408  0.065 . 
## - MdPh  1 22.555 2.2531  0.050 * 
## - PH    1 25.491 4.4692  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PSU + PH + IL + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - PSU   1 21.790 1.6242  0.150   
## - IL    1 21.885 1.6916  0.150   
## - sqTN  1 21.640 1.5180  0.145   
## - MdPh  1 22.804 2.3672  0.045 * 
## - PH    1 25.445 4.5555  0.010 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PH + IL + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - IL    1 22.052 1.6276   0.14   
## - sqTN  1 22.478 1.9623   0.07 . 
## - MdPh  1 22.825 2.2420   0.03 * 
## - PH    1 24.876 4.0363   0.01 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: dat ~ PH + MdPh + sqTN 
## 
##        Df    AIC      F Pr(>F)   
## - MdPh  1 22.653 2.0825  0.075 . 
## - sqTN  1 23.014 2.4014  0.015 * 
## - PH    1 24.547 3.8428  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

res2.rda

## Call: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2,
## distance = "horn")
## 
##               Inertia Proportion Rank RealDims
## Total          4.6471     1.0000              
## Constrained    1.9129     0.4116    3        3
## Unconstrained  2.7341     0.5884   11        8
## Inertia is squared Horn distance 
## 
## Eigenvalues for constrained axes:
## dbRDA1 dbRDA2 dbRDA3 
## 0.9743 0.6191 0.3195 
## 
## Eigenvalues for unconstrained axes:
##    MDS1    MDS2    MDS3    MDS4    MDS5    MDS6    MDS7    MDS8   iMDS1 
##  1.0463  0.6487  0.4820  0.3541  0.2177  0.1143  0.0509  0.0016 -0.0159 
##   iMDS2   iMDS3 
## -0.0579 -0.1077

plot(res2.rda)

anova.cca(res2.rda)

## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)    
## Model     3   1.9129 2.5654  0.001 ***
## Residual 11   2.7341                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova.cca(res2.rda, by="axis")

## Permutation test for dbrda under reduced model
## Forward tests for axes
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)  
## dbRDA1    1  0.97431 3.9199  0.014 *
## dbRDA2    1  0.61912 2.4909  0.106  
## dbRDA3    1  0.31950 1.2854  0.283  
## Residual 11  2.73412                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova.cca(res2.rda, by="terms")

## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = dat ~ PH + MdPh + sqTN, data = edat2, distance = "horn")
##          Df SumOfSqs      F Pr(>F)   
## PH        1  0.85141 3.4254  0.010 **
## MdPh      1  0.46465 1.8694  0.080 . 
## sqTN      1  0.59687 2.4014  0.018 * 
## Residual 11  2.73412                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

As a result of reducing the number of explanation variables, the value of Constrained Proportion decreases, but dbRDA is now significant. Coefficients of variables can be seen by a function coef().

coef(res2.rda)

##           dbRDA1      dbRDA2     dbRDA3
## PH    0.27375601 -0.02673265  0.0433661
## MdPh -0.03684029  0.22352641  0.1695031
## sqTN  0.10816218  0.22639967 -0.1415342

---

References

• Borcard D, Gille F, Legendre P (2011) Numerical Ecology with R. Springer
• Legendre P, Legendre L (2012) Numerical Ecology. Elsevier
• Magurran AE, McGill BJ (2011) Biological Diversity. Oxford Univ.

Back to "Analyses on Ecological Assemblages by R"