Title: Basic clustering visualizations using Rao’s quadratic entropy

Name: Tammy L. Elliott

Date: April 8, 2016

R version 3.1

Gower’s distance calculation - four variables

# Drop soil moisture based on previous collinearity analyses
site.dist<-cbind(site.elev,site.vis, site.slope,site.soil)

#Calculate Gower's distance; the results dissimilarity distances
site.env.dist<-gowdis(site.dist)

Justification for using Rao’s quadratic entropy and how I applied it

# 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)

---

Rao’s quadratic entropy calculations

#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
raoD.pa<-raoD(pa.sp.allsp, phy=NULL)
raoD.pa.dist<-dist(raoD.pa$H)
#Calculate Rao's quadratic entropy on abundance data, with no phylogeny
raoD.abd<-raoD(abd.sp.allsp, phy=NULL)
raoD.abd.dist<-dist(raoD.abd$H)

#Phylogenetic analyses
raoD.pa.phy<-raoD(pa.sp.allsp, phy=trans.one.tree)

## [1] "Dropping tips from the tree because they are not present in the community data:"
## [1] "CorAlt" "CorSto" "HydArb" "ComUmb" "TheHum" "XimAme"

raoD.pa.phy.dist<-dist(raoD.pa.phy$H)
raoD.abd.phy<-raoD(abd.sp.allsp, phy=trans.one.tree)

## [1] "Dropping tips from the tree because they are not present in the community data:"
## [1] "CorAlt" "CorSto" "HydArb" "ComUmb" "TheHum" "XimAme"

raoD.abd.phy.dist<-dist(raoD.abd.phy$H)

---

Two clusters, presence-absence

#Plot this relationship; site.env distance and 2 cluster presence-absence and phylogenetic presence-absence
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.2.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.2.cmdscale, site.env.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.2.cmdscale, site.env.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(site.env.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

raoD.pa.dist.2.cmdscale<-cmdscale(raoD.pa.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.dist.2.cmdscale, raoD.pa.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.dist.2.cmdscale, raoD.pa.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.raoD.pa.dist<-cascadeKM((raoD.pa.dist),2,15)
plot(ccas.raoD.pa.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(raoD.pa.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

raoD.pa.phy.dist.2.cmdscale<-cmdscale(raoD.pa.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.phy.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.phy.dist.2.cmdscale, raoD.pa.phy.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.phy.dist.2.cmdscale, raoD.pa.phy.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.raoD.pa.phy.dist<-cascadeKM((raoD.pa.phy.dist),2,15)
plot(ccas.raoD.pa.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(raoD.pa.phy.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

Two clusters, abundance

#Plot relationship with abundance-weighted Rao's D and phylogenetic version
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.2.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.2.cmdscale, site.env.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.2.cmdscale, site.env.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(site.env.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

raoD.abd.dist.2.cmdscale<-cmdscale(raoD.abd.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.dist.2.cmdscale, raoD.abd.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.dist.2.cmdscale, raoD.abd.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.dist<-cascadeKM((raoD.abd.dist),2,15)
plot(ccas.raoD.abd.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(raoD.abd.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

raoD.abd.phy.dist.2.cmdscale<-cmdscale(raoD.abd.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.phy.dist.2.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.phy.dist.2.cmdscale, raoD.abd.phy.dist.cfuz.2$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.phy.dist.2.cmdscale, raoD.abd.phy.dist.cfuz.2$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.phy.dist<-cascadeKM((raoD.abd.phy.dist),2,15)
plot(ccas.raoD.abd.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black",  "gray85"), as.numeric(raoD.abd.phy.dist.cfuz.2$cluster))
colfunc <- colorRampPalette(c("black", "gray85"))
#par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomright", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black",  "gray85"),col="black",  cex=0.75,  bty="n")

Six clusters, presence-absence

#Plot this relationship; site.env distance and 6 cluster presence-absence and phylogenetic presence-absence
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.6.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.6.cmdscale, site.env.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.6.cmdscale, site.env.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(site.env.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
    bg=cols,pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.pa.dist.6.cmdscale<-cmdscale(raoD.pa.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.dist.6.cmdscale, raoD.pa.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.dist.6.cmdscale, raoD.pa.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.raoD.pa.dist<-cascadeKM((raoD.pa.dist),2,15)
plot(ccas.raoD.pa.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(raoD.pa.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.pa.phy.dist.6.cmdscale<-cmdscale(raoD.pa.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.phy.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.phy.dist.6.cmdscale, raoD.pa.phy.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.phy.dist.6.cmdscale, raoD.pa.phy.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.raoD.pa.phy.dist<-cascadeKM((raoD.pa.phy.dist),2,15)
plot(ccas.raoD.pa.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(raoD.pa.phy.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

Six clusters, abundance

#Plot relationship with abundance-weighted Rao's D and phylogenetic version
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.6.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.6.cmdscale, site.env.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.6.cmdscale, site.env.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(site.env.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.abd.dist.6.cmdscale<-cmdscale(raoD.abd.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.dist.6.cmdscale, raoD.abd.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.dist.6.cmdscale, raoD.abd.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.dist<-cascadeKM((raoD.abd.dist),2,15)
plot(ccas.raoD.abd.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(raoD.abd.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.abd.phy.dist.6.cmdscale<-cmdscale(raoD.abd.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.phy.dist.6.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.phy.dist.6.cmdscale, raoD.abd.phy.dist.cfuz.6$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.phy.dist.6.cmdscale, raoD.abd.phy.dist.cfuz.6$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.phy.dist<-cascadeKM((raoD.abd.phy.dist),2,15)
plot(ccas.raoD.abd.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray20", "gray40", "gray60", "gray80", "white"), as.numeric(raoD.abd.phy.dist.cfuz.6$cluster))
colfunc <- colorRampPalette(c("black", "gray20", "gray40", "gray60", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray20"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6"), pch=21,
    pt.bg=c("gray40","gray60", "gray80", "white"), col="black",  cex=0.75,  bty="n")

Seven clusters, presence-absence

#Plot this relationship; site.env distance and 6 cluster presence-absence and phylogenetic presence-absence
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.7.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.7.cmdscale, site.env.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.7.cmdscale, site.env.dist.cfuz.7$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(site.env.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
    bg=cols,pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.pa.dist.7.cmdscale<-cmdscale(raoD.pa.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.dist.7.cmdscale, raoD.pa.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.dist.7.cmdscale, raoD.pa.dist.cfuz.7$cluster, col="grey70")

## Warning in chol.default(cov, pivot = TRUE): the matrix is either rank-
## deficient or indefinite

box(lwd=1)

ccas.raoD.pa.dist<-cascadeKM((raoD.pa.dist),2,15)
plot(ccas.raoD.pa.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(raoD.pa.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.pa.phy.dist.7.cmdscale<-cmdscale(raoD.pa.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.pa.phy.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.pa.phy.dist.7.cmdscale, raoD.pa.phy.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.pa.phy.dist.7.cmdscale, raoD.pa.phy.dist.cfuz.7$cluster, col="grey70")
box(lwd=1)

ccas.raoD.pa.phy.dist<-cascadeKM((raoD.pa.phy.dist),2,15)
plot(ccas.raoD.pa.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(raoD.pa.phy.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")

Seven clusters, abundance

#Plot relationship with abundance-weighted Rao's D and phylogenetic version
#dev.new(width=11.8, height=11.8)
par(mfrow=c(3,3))
par(mar=c(5.1,4.1,4.1,2.1))
site.env.dist.7.cmdscale<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(site.env.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(site.env.dist.7.cmdscale, site.env.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(site.env.dist.7.cmdscale, site.env.dist.cfuz.7$cluster, col="grey70")
box(lwd=1)

ccas.site.env.dist<-cascadeKM((site.env.dist),2,15)
plot(ccas.site.env.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(site.env.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.abd.dist.7.cmdscale<-cmdscale(raoD.abd.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.dist.7.cmdscale, raoD.abd.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.dist.7.cmdscale, raoD.abd.dist.cfuz.7$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.dist<-cascadeKM((raoD.abd.dist),2,15)
plot(ccas.raoD.abd.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(raoD.abd.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")

raoD.abd.phy.dist.7.cmdscale<-cmdscale(raoD.abd.phy.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
plot(raoD.abd.phy.dist.7.cmdscale, type="n", xlab="", ylab="", main="", cex=1, pch=16, col="black", xaxt='n', yaxt="n")
ordispider(raoD.abd.phy.dist.7.cmdscale, raoD.abd.phy.dist.cfuz.7$cluster, col="grey40", label=TRUE)
ordiellipse(raoD.abd.phy.dist.7.cmdscale, raoD.abd.phy.dist.cfuz.7$cluster, col="grey70")
box(lwd=1)

ccas.raoD.abd.phy.dist<-cascadeKM((raoD.abd.phy.dist),2,15)
plot(ccas.raoD.abd.phy.dist,sortq=TRUE)

cols <- myColorRamp(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"), as.numeric(raoD.abd.phy.dist.cfuz.7$cluster))
colfunc <- colorRampPalette(c("black", "gray18", "gray35", "gray48","grey63", "gray80", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.coord$Latitude..WGS.84.datum., site.coord$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",bg=cols,pch=21, 
    col="black",cex=1.75, cex.lab=1.5, cex.axis=1, main="", cex.main=1.25)
#xl <- 1
#yb <- -2.3
#xr <- 1.5
#yt <- 2.7
legend("bottomleft", c("Group 1","Group 2"), pch=21,
    pt.bg=c("black", "gray18"),col="black",  cex=0.75,  bty="n")
legend("bottomright", c("Group 3", "Group 4", "Group5", "Group6", "Group7"), pch=21,
    pt.bg=c( "gray35", "gray48","grey63", "gray80", "white"), col="black",  cex=0.75,  bty="n")