Rao’s quadratic entropy - explanation
# Rao's quadratic entropy
# 1. calculate presence/absence based metrics by first creating a 0-1 presence-absence species-site dataframe out of abundance dataframe. I then use this as input dataframe for Rao's quadratic entropy formula in picante
# 2. Phylogeny is ultrametric - I made my BEAST phylogeny ultrametric by using the chrons function in the ape package; lambda=0.1, correlated model - this produced a tree with similar
# branch lengths to the BEAST output phylogeny; ultrametric phylogeny is required for raoD function (see Pavoine et al. 2005)
# 3. For Rao's quadratic entropy - Rao’s quadratic entropy summarizes the expected
# dissimilarity between two sample plots randomly drawn from the set of N plots (Ricotto and Marignani 2007) or another way of saying it -
# If no phylogeny is supplied, Dkk is equivalent to Simpson’s diversity (probability that two individuals
# drawn from a community are from different taxa)
# 4. Phylogenetic Rao's quadratic entropy - MPD based metric of phylogenetic beta diversity; that is - Mean phylogenetic distance
# between a species from assemblage 1 and a species from assemblage 2
# 5. I am extracting the H values, which are the Among-community diversities excluding within-community diversity, this gives a diagonal of zero as it is a dissimilarity index (just like the Gower's index)
#Calculate Rao's quadratic entropy on presence-absence data, with no phylogeny
# $H gives a dissimilarity matrix with 0's as the diagonals; convert to distance object for subsequent analyses
Plot four cluster presence-absence indices
par(mfrow=c(3,2))
par(mar=c(3.5,3.5,4.1,4.1))
plot(site.env.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
site.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = site.env.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
legend("topright", legend = paste( c('1', '2', '3', '4')), pch = 16, col=c("grey80", "gray55", "gray25", "black"),
cex=1.4,bty='n',inset=c(0.01,0.05))
par(mar=c(3.5,3.5,4.1,4.1))
plot(raoD.pa.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.dist.4.cmdscale, raoD.pa.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.dist.4.cmdscale, raoD.pa.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
vasc.pa.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = raoD.pa.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
par(mar=c(3.5,3.5,4.1,4.1))
plot(raoD.pa.phy.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.phy.dist.4.cmdscale, raoD.pa.phy.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.phy.dist.4.cmdscale, raoD.pa.phy.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
vasc.pa.phy.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = raoD.pa.phy.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)", main = "", ticktype = "simple", theta = 35, d = 5,colkey = FALSE))
#textClick.bold("(a)",cex=c(size=1.5, font=2) )
#textClick.bold("(b)", cex=c(size=1.5, font=2))
#textClick.bold("(c)", cex=c(size=1.5, font=2))
#textClick.bold("(d)", cex=c(size=1.5, font=2))
#textClick.bold("(e)", cex=c(size=1.5, font=2))
#textClick.bold("(f)", cex=c(size=1.5, font=2))
Plot four cluster abundance indices
par(mfrow=c(3,2))
par(mar=c(3.5,3.5,4.1,4.1))
plot(site.env.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
site.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = site.env.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
legend("topright", legend = paste( c('1', '2', '3', '4')), pch = 16, col=c("grey80", "gray55", "gray25", "black"),
cex=1.4,bty='n',inset=c(0.01,0.05))
par(mar=c(3.5,3.5,4.1,4.1))
plot(raoD.abd.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.dist.4.cmdscale, raoD.abd.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.dist.4.cmdscale, raoD.abd.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
vasc.abd.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = raoD.abd.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
par(mar=c(3.5,3.5,4.1,4.1))
plot(raoD.abd.phy.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.phy.dist.4.cmdscale, raoD.abd.phy.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.phy.dist.4.cmdscale, raoD.abd.phy.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
vasc.abd.phy.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = raoD.abd.phy.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)", main = "", ticktype = "simple", theta = 35, d = 5,colkey = FALSE))
#textClick.bold("(a)",cex=c(size=1.5, font=2) )
#textClick.bold("(b)", cex=c(size=1.5, font=2))
#textClick.bold("(c)", cex=c(size=1.5, font=2))
#textClick.bold("(d)", cex=c(size=1.5, font=2))
#textClick.bold("(e)", cex=c(size=1.5, font=2))
#textClick.bold("(f)", cex=c(size=1.5, font=2))
Angiosperm only analyses
Plot four cluster presence-absence indices for angiosperms
par(mfrow=c(3,2))
par(mar=c(3.5,3.5,4.1,4.1))
plot(site.env.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
site.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = site.env.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
legend("topright", legend = paste( c('1', '2', '3', '4')), pch = 16, col=c("grey80", "gray55", "gray25", "black"),
cex=1.4,bty='n',inset=c(0.01,0.05))
par(mar=c(3.5,3.5,4.1,4.1))
plot(angio.raoD.pa.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(angio.raoD.pa.dist.4.cmdscale, angio.raoD.pa.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(angio.raoD.pa.dist.4.cmdscale, angio.raoD.pa.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
angio.pa.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = angio.raoD.pa.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
par(mar=c(3.5,3.5,4.1,4.1))
plot(angio.raoD.pa.phy.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(angio.raoD.pa.phy.dist.4.cmdscale, angio.raoD.pa.phy.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(angio.raoD.pa.phy.dist.4.cmdscale, angio.raoD.pa.phy.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
angio.pa.phy.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = angio.raoD.pa.phy.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)", main = "", ticktype = "simple", theta = 35, d = 5,colkey = FALSE))
#textClick.bold("(a)",cex=c(size=1.5, font=2) )
#textClick.bold("(b)", cex=c(size=1.5, font=2))
#textClick.bold("(c)", cex=c(size=1.5, font=2))
#textClick.bold("(d)", cex=c(size=1.5, font=2))
#textClick.bold("(e)", cex=c(size=1.5, font=2))
#textClick.bold("(f)", cex=c(size=1.5, font=2))
Plot four cluster abundance indices for angiosperms
par(mfrow=c(3,2))
par(mar=c(3.5,3.5,4.1,4.1))
plot(site.env.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.4.cmdscale, site.env.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
site.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = site.env.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
legend("topright", legend = paste( c('1', '2', '3', '4')), pch = 16, col=c("grey80", "gray55", "gray25", "black"),
cex=1.4,bty='n',inset=c(0.01,0.05))
par(mar=c(3.5,3.5,4.1,4.1))
plot(angio.raoD.abd.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(angio.raoD.abd.dist.4.cmdscale, angio.raoD.abd.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(angio.raoD.abd.dist.4.cmdscale, angio.raoD.abd.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
angio.abd.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = angio.raoD.abd.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)",
main = "", ticktype = "simple", theta = 35, d = 5,
colkey = FALSE))
par(mar=c(3.5,3.5,4.1,4.1))
plot(angio.raoD.abd.phy.dist.4.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(angio.raoD.abd.phy.dist.4.cmdscale, angio.raoD.abd.phy.dist.cfuz.4$cluster, col="grey40", label=TRUE)
ordiellipse(angio.raoD.abd.phy.dist.4.cmdscale, angio.raoD.abd.phy.dist.cfuz.4$cluster, col="grey70")
box(lwd=1)
par(mar=c(1,0.25,1,1))
angio.abd.phy.dist.4.3d<-with(site.coord, scatter3D(x = site.lon[,1], y = site.lat[,1], z = site.elev[,1], colvar = angio.raoD.abd.phy.dist.cfuz.4$clustering,
col=c("grey80", "gray55", "gray25", "black"),pch = 16, cex = 1.25, cex.lab=1.35, cex.axis=1,xlab = "Longitude", ylab = "Latitude",
zlab = "Elevation (m)", main = "", ticktype = "simple", theta = 35, d = 5,colkey = FALSE))
#textClick.bold("(a)",cex=c(size=1.5, font=2) )
#textClick.bold("(b)", cex=c(size=1.5, font=2))
#textClick.bold("(c)", cex=c(size=1.5, font=2))
#textClick.bold("(d)", cex=c(size=1.5, font=2))
#textClick.bold("(e)", cex=c(size=1.5, font=2))
#textClick.bold("(f)", cex=c(size=1.5, font=2))