Worked example 5

WE5 provides R code to perform the worked examples in Pélissier & Couteron (2007, Appendix), using standard functions of the base library or other libraries (stats, ade4, vegan) available at http://cran.r-project.org/, and dedicated routines the diversity library available at http://pelissier.free.fr/Diversity.html.

1 - Data example

spot is the table of species occurrences (Fig. 1a):

spot<-cbind(c(1,0,0,1,0,0,0),c(0,1,1,0,1,0,0),c(0,0,0,0,0,1,1))
spot

#    [,1] [,2] [,3]
#[1,]    1    0    0
#[2,]    0    1    0
#[3,]    0    1    0
#[4,]    1    0    0
#[5,]    0    1    0
#[6,]    0    0    1
#[7,]    0    0    1


hat is the vector of discrete habitat types (Fig. 1b):

hat<-as.factor(c(1,2,2,3,3,3,3))
hat

#[1] 1 2 2 3 3 3 3
#Levels: 1 2 3


ph is the vector of soil pH variable ("CONTINUOUS ENVIRONMENTAL GRADIENT"):

ph<-c(4.6,5.3,5.3,5.8,5.8,5.8,5.8)
ph

#[1] 4.6 5.3 5.3 5.8 5.8 5.8 5.8

n is the total number of observations:

n<-sum(spot)
n

#[1] 7


wrich and wshan are vectors of species weights for the richness and Shannon metrics, respectively (for Simpson metric the weight is 1 for all species):

wrich<-1/apply(spot,2,mean)
wrich

#[1] 3.500000 2.333333 3.500000

wshan<-log(1/apply(spot,2,mean))/(1-apply(spot,2,mean))
wshan

#[1] 1.753868 1.482771 1.753868

2 - Partitioning species diversity according to discrete habitat types (Tab. 1)

The standard aov{base}function performs (M)ANOVA, while the summary.div{diversity}converts sum of squares into diversity measures:

library(diversity)
summary.div(aov(spot~hat))
#class:  summary.anodiv
#call: aov(formula = spot ~ hat)
#mod:  1
#
#$Richness:
#            Df Diversity F-value
#Total          2               
#hat         2  0.875     1.56  
#Residuals   4  1.13            
#
#$Shannon:
#            Df Diversity F-value
#Total          1.08            
#hat         2  0.482     1.61  
#Residuals   4  0.597           
#
#$Simpson:
#            Df Diversity F-value
#Total          0.653           
#hat         2  0.296     1.66  
#Residuals   4  0.357           
#
#Do you want to print by-species results ? (y/n) :

y
#Response $1 (n=2):
#            Df Variance F-value
#Total          0.204          
#hat         2  0.0969   1.81  
#Residuals   4  0.107          
#
#Response $2 (n=3):
#            Df Variance F-value
#Total          0.245          
#hat         2  0.138    2.57  
#Residuals   4  0.107       
#
#Response $3 (n=2):
#            Df Variance F-value
#Total          0.204          
#hat         2  0.0612   0.857 
#Residuals   4  0.143

The anodiv{diversity}function performs MANOVA from a specific odf object more efficient for large data sets. Moreover, summary{diversity}performs randomisation tests of statistical significance:

ovec<-odf(spot)
summary(anodiv(ovec$sp~hat),nrepet=100)

#class:  summary.anodiv
#call: anodiv.formula(formula = ovec$sp ~ hat)
#mod:  1
#Monte Carlo tests based on 100 replicates
#
#$Richness:
#          Df Diversity F-value Pr(>F)
#Total        2                            
#hat       2  0.875     1.56    0.257
#Residuals 4  1.12                    
#
#$Shannon:
#          Df Diversity F-value Pr(>F)
#Total        1.08                    
#hat       2  0.482     1.61    0.406
#Residuals 4  0.597                   
#
#$Simpson:
#          Df Diversity F-value Pr(>F)
#Total        0.653                   
#hat       2  0.296     1.66    0.257
#Residuals 4  0.357                   
#
#Do you want to print by-species results ? (y/n) :

y
#Response $Sp1 (n=2):
#          Df Variance F-value Pr(>F)
#Total        0.204                  
#hat       2  0.0969   1.81    0.99  
#Residuals 4  0.107                  
#
#Response $Sp2 (n=3):
#          Df Variance F-value Pr(>F)
#Total        0.245                  
#hat       2  0.138    2.57    0.406
#Residuals 4  0.107                  
#
#Response $Sp3 (n=2):
#          Df Variance F-value Pr(>F)
#Total        0.204                  
#hat       2  0.0612   0.857   0.446
#Residuals 4  0.143


3 - Partitioning species diversity according to pH gradient (Tab. 2)

The standard aov{base}function performs (M)LM and (M)ANOVA, while the summary.div{diversity}converts sum of squares into diversity measures (randomisation tests for continuous covariates are not yet implemented in package diversity, but see for instance Anderson 2001): 

summary.div(aov(spot~ph))
#class:  summary.anodiv
#call: aov(formula = spot ~ ph)
#mod:  1
#
#$Richness:
#            Df Diversity F-value
#Total          2               
#ph          1  0.29      0.847 
#Residuals   5  1.71            
#
#$Shannon:
#            Df Diversity F-value
#Total          1.08            
#ph          1  0.145     0.778 
#Residuals   5  0.934           
#
#$Simpson:
#            Df Diversity F-value
#Total          0.653           
#ph          1  0.0829    0.727 
#Residuals   5  0.57            
#
#Do you want to print by-species results ? (y/n) :

y
#Response $1 (n=2):
#            Df Variance F-value
#Total          0.204          
#ph          1  0.0374   1.12  
#Residuals   5  0.167          
#
#Response $2 (n=3):
#            Df Variance F-value
#Total          0.245          
#ph          1  0.000374 0.00764
#Residuals   5  0.245          
#
#Response $3 (n=2):
#            Df Variance F-value
#Total          0.204          
#ph          1  0.0452   1.42  
#Residuals   5  0.159


4 - Computing MSS from inertia of an abundance data table (Tab. 1)

corresp{MASS}function performs Correspondence Analysis (CA) providing MSS in the richness metric: 

A<-apply(spot,2,function(x) tapply(x,hat,sum))
library(MASS)
sum(corresp(A,nf=2)$cor^2)

#[1] 0.875


dudi.coa{ade4}and dudi.nsc{ade4}perform CA and Non-Symmetric CA providing MSS in both richness and Simpson metrics, respectively: 

library(ade4)
sum(dudi.coa(as.data.frame(A),scannf=F)$eig)

#[1] 0.875

sum(dudi.nsc(as.data.frame(A),scannf=F)$eig/ncol(A))
#[1] 0.2959184


ca.richness{diversity}, nsca.simpson{diversity} and cwca.shannon{diversity} are wrapper functions to as.dudi{ade4}, which perform MSS in the richness, Shannon and Simpson metrics, respectively:

rich<-ca.richness(A)
#Select the number of axes:
2

summary(rich)
#class:  summary.dudiv
#metric: Richness
#call: ca.richness.default(Y = A)
#
#total diversity: 2
#explained diversity: 0.875
#ratio of explained diversity: 0.437
#Pr(>ratio): 0.371 based on 999 replicates

shan<-cwca.shannon(A)
#Select the number of axes: 
2

summary(shan)
#class:  summary.dudiv
#metric: Shannon
#call: cwca.shannon.default(Y = A)
#
#total diversity: 1.08
#explained diversity: 0.482
#ratio of explained diversity: 0.446
#Pr(>ratio): 0.371 based on 999 replicates

simp<-nsca.simpson(A)
#Select the number of axes:
2

summary(simp)
#class:  summary.dudiv
#metric: Simpson
#call: nsca.simpson.default(Y = A)
#
#total diversity: 0.653
#explained diversity: 0.296
#ratio of explained diversity: 0.453 


5 - Computing MSS from dissimilarity matrices (Fig. 3)

dist{base}function computes Euclidean distances:

ni<-apply(A,1,sum)
w<-c(ni[1]*ni[2],ni[1]*ni[3],ni[2]*ni[3])
B<-apply(A,2,function(x) x/ni)
D2<-apply(B,2,dist)^2
rich.ssq<-sum(apply(t(apply(D2,1,function(x) x*wrich)),2,function(x) x*w))/n^2
shan.ssq<-sum(apply(t(apply(D2,1,function(x) x*wshan)),2,function(x) x*w))/n^2
simp.ssq<-sum(D2*w)/n^2
cbind(rich.ssq,shan.ssq,simp.ssq)

#     rich.ssq  shan.ssq  simp.ssq
#[1,]    0.875 0.4816568 0.2959184


dissim{diversity}function computes dissimilarities:

d<-dissim(A)
D2<-d$delta
apply(d$dissim,2,function(x) sum(x*d$weight))

# richness   Shannon   Simpson
#0.8750000 0.4816568 0.2959184

 

Literature cited

Anderson, M.J. 2001. Permutation tests for univariate or multivariate analysis of variance and regression. Canadian Journal of Fisheries and Aquatic Science, 58:626-639.

Pélissier, R. and Couteron, P. 2007. An operational, addittive framework for species diversity partitioning and beta-diversity analysis. Journal of Ecology, 95:294-300.