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