Schefferville transitions - Local contribution to beta diversity on clusters

Tammy L. Elliott

March 4, 2015

R version 2.12 and 3.01 (vers. 2012 for constrained.clust and mvpart)


Site characteristics data preparation

#Truncate site data with different environmental variables
site.lat<-site[,4,drop=FALSE]
site.lon<-site[,5,drop=FALSE]
site.elev<-site[,6,drop=FALSE]
site.vis<-site[,10,drop=FALSE]
site.slope<-site[,8,drop=FALSE]
site.soil.dep.ind<-site[,c(14:17,drop=FALSE)]
site.soil<-apply(site.soil.dep.ind,1,mean)

#Calculate soil moist
soil.moist<-site$Soil.moisture
soil.moist.hydric<-gsub("hydric", "1", soil.moist)
soil.moist.mesic.hydric<-gsub("mesic-1", "2", soil.moist.hydric)
soil.moist.mesic<-gsub("mesic", "3", soil.moist.mesic.hydric)
soil.moist.mesic.xeric<-gsub("xeric-3", "4", soil.moist.mesic)
soil.moist.xeric<-gsub("xeric", "5", soil.moist.mesic.xeric)

site.moist<-as.numeric(soil.moist.xeric)

#Create matrix of variables 
site.env<-cbind(site.elev,site.vis, site.slope,site.soil, site.moist)
head(site.env)
##     Elevation..m. Sky.visible.... Slope..Degrees. site.soil site.moist
## EA1           523              95               0    22.250          2
## EB1           522             100               0    22.500          2
## EC1           524             100               0    14.500          2
## ED1           532              80               0    18.125          2
## EE1           522             100               0    26.375          1
## EF1           531              80               0    30.000          2
#Calculate Gower's distance; the results are a dissimilarity matrix
site.env.dist<-gowdis(site.env)

K-means clustering of environmental distances

#Now do a K-means clustering analysis of these distances

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

#Try K-means clustering with 4-6 clusters

#vasc.env.clust.4<-kmeans(site.env.dist,4)
#saveRDS(vasc.env.clust.4, file="vasc.env.clust.4.rds")
vasc.env.clust.4<-readRDS("vasc.env.clust.4.rds")

#vasc.env.clust.5<-kmeans(site.env.dist,5)
#saveRDS(vasc.env.clust.5, file="vasc.env.clust.5.rds")
vasc.env.clust.5<-readRDS("vasc.env.clust.5.rds")

#vasc.env.clust.6<-kmeans(site.env.dist,6)
#saveRDS(vasc.env.clust.6, file="vasc.env.clust.6.rds")
vasc.env.clust.6<-readRDS("vasc.env.clust.6.rds")

#Plot PCoA with 4 groupings based on k-means clustering for environmental distances
#dev.new(width=11.8, height=8)
par(mfrow=c(1,2))
vasc.env.4<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.4[,1]
y<-vasc.env.4[,2]
plot(vasc.env.4, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with four \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.4, vasc.env.clust.4$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.4, vasc.env.clust.4$cluster, col="grey70")
box(lwd=2)
text(x, y, rownames(vasc.env.4), cex=0.4)

vasc.env.4<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.4[,1]
y<-vasc.env.4[,2]
plot(vasc.env.4, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with four \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.4, vasc.env.clust.4$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.4, vasc.env.clust.4$cluster, col="grey70")
box(lwd=2)

#Plot PCoA with 5 groupings based on k-means clustering for environmental distances
#dev.new(width=11.8, height=8)
par(mfrow=c(1,2))
vasc.env.5<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.5[,1]
y<-vasc.env.5[,2]
plot(vasc.env.5, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with five \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.5, vasc.env.clust.5$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.5, vasc.env.clust.5$cluster, col="grey70")
box(lwd=2)
text(x, y, rownames(vasc.env.5), cex=0.4)

vasc.env.5<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.5[,1]
y<-vasc.env.5[,2]
plot(vasc.env.5, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with five \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.5, vasc.env.clust.5$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.5, vasc.env.clust.5$cluster, col="grey70")
box(lwd=2)

#Plot PCoA with 6 groupings based on k-means clustering for environmental distances
#dev.new(width=11.8, height=8)
par(mfrow=c(1,2))
vasc.env.6<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.6[,1]
y<-vasc.env.6[,2]
plot(vasc.env.6, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with six \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.6, vasc.env.clust.6$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.6, vasc.env.clust.6$cluster, col="grey70")
box(lwd=2)
text(x, y, rownames(vasc.env.6), cex=0.4)

vasc.env.6<-cmdscale(site.env.dist, k = 54, eig = FALSE, add = FALSE, x.ret = FALSE)
x<-vasc.env.6[,1]
y<-vasc.env.6[,2]
plot(vasc.env.6, type="n", xlab="", ylab="", main="Principal Coordinate Analysis (PCoA) with six \nclusters based on K-means clustering on \npairwise environmental distances", cex=1, pch=16, col="black")
ordispider(vasc.env.6, vasc.env.clust.6$cluster, col="grey40", label=TRUE)
ordiellipse(vasc.env.6, vasc.env.clust.6$cluster, col="grey70")
box(lwd=2)

### K-means cluster attributes
# Show structure of 5 and 6 clusters
vasc.env.clust.5$cluster
##     EA1     EB1     EC1     ED1     EE1     EF1     EG1     EH1     EI1 
##       1       1       4       1       1       1       1       4       1 
##     EJ1     EK1     EA2     EB2     EC2     ED2     EE2     EF2     EG2 
##       4       1       4       4       1       1       1       4       1 
##     EH2     EI2     EJ2     EK2     EA3     EB3     EC3     ED3     EE3 
##       4       4       4       4       1       4       1       4       1 
##     EF3     EG3     EH3     EI3     EJ3     EK3     EA4     EB4     EC4 
##       1       1       4       1       1       4       1       4       1 
##     ED4     EE4     EF4     EG4     EH4     EI4     EJ4     EK4     EA5 
##       1       4       1       4       4       4       1       4       4 
##     EB5     EC5     ED5     EE5     EF5     EG5     EH5     EI5     EJ5 
##       4       4       1       4       1       4       4       4       4 
##     EK5     EA6     EB6     EC6     ED6     EE6     EF6     EG6     EH6 
##       4       4       4       4       4       4       4       5       1 
##     EI6     EJ6     EK6     EA7     EB7     EC7     ED7     EE7     EF7 
##       4       4       4       1       4       1       1       4       4 
##     EG7     EH7     EI7     EJ7     EK7     EA8     EB8     EC8     ED8 
##       4       4       5       4       4       4       4       1       1 
##     EE8     EF8     EG8     EH8     EI8     EJ8     EK8  E600A1  E600B2 
##       4       1       1       4       4       5       5       4       4 
##  E600C3  E600D4  E600E5  E600F6  E600G7  E600H8  E600I9 E600J10 E600K11 
##       5       5       5       5       5       5       5       5       5 
##  E615A1  E615B2  E615C3  E615D4  E615E5  E615F6  E615G7  E615H8  E615I9 
##       4       5       5       5       5       2       5       5       5 
## E615J10 E615K11  E645A1  E645B2  E645C3  E645D4  E645E5  E645F6  E645G7 
##       5       2       5       5       5       2       2       2       5 
##  E645H8  E645I9 E645J10 E645K11  E675A1  E675B2  E675C3  E675D4  E675E5 
##       2       5       5       5       5       2       2       2       2 
##  E675F6  E675G7  E675H8  E675I9 E675J10 E675K11  E710A1  E710B2  E710C3 
##       2       3       2       2       2       2       2       2       2 
##  E710D4  E710E5  E710F6  E710G7  E710H8  E710I9 E710J10 E710K11  E745A1 
##       5       2       2       5       2       2       2       3       2 
##  E745B2  E745C3  E745D4  E745E5  E745F6  E745G7  E745H8  E745I9 E745J10 
##       3       2       3       3       3       3       2       3       3 
## E745K11  E775A1  E775B2  E775C3  E775D4  E775E5  E775F6  E775G7  E775H8 
##       3       3       2       2       3       3       3       3       3 
##  E775I9 E775J10 E775K11  E800A1  E800B2  E800C3  E800D4  E800E5  E800F6 
##       3       3       3       3       3       2       3       3       3 
##  E800G7  E800H8  E800I9 E800J10 E800K11 
##       3       3       3       3       3
vasc.env.clust.5$totss
## [1] 586.5044
vasc.env.clust.5$withinss
## [1]  9.450085 18.481272 17.924102 16.956836 12.062646
vasc.env.clust.5$tot.withinss
## [1] 74.87494
vasc.env.clust.5$betweenss
## [1] 511.6295
vasc.env.clust.5$size
## [1] 34 29 29 53 31
vasc.env.clust.6$cluster
##     EA1     EB1     EC1     ED1     EE1     EF1     EG1     EH1     EI1 
##       1       1       1       1       3       3       1       1       3 
##     EJ1     EK1     EA2     EB2     EC2     ED2     EE2     EF2     EG2 
##       1       1       1       1       3       3       1       1       3 
##     EH2     EI2     EJ2     EK2     EA3     EB3     EC3     ED3     EE3 
##       1       1       1       1       1       4       3       1       1 
##     EF3     EG3     EH3     EI3     EJ3     EK3     EA4     EB4     EC4 
##       3       3       4       3       3       4       3       4       1 
##     ED4     EE4     EF4     EG4     EH4     EI4     EJ4     EK4     EA5 
##       3       4       3       4       4       1       3       4       4 
##     EB5     EC5     ED5     EE5     EF5     EG5     EH5     EI5     EJ5 
##       4       4       3       4       3       4       4       4       4 
##     EK5     EA6     EB6     EC6     ED6     EE6     EF6     EG6     EH6 
##       4       1       1       1       4       4       4       2       3 
##     EI6     EJ6     EK6     EA7     EB7     EC7     ED7     EE7     EF7 
##       4       4       4       3       4       3       3       4       4 
##     EG7     EH7     EI7     EJ7     EK7     EA8     EB8     EC8     ED8 
##       4       4       2       4       4       4       4       1       3 
##     EE8     EF8     EG8     EH8     EI8     EJ8     EK8  E600A1  E600B2 
##       4       3       1       4       4       2       2       4       4 
##  E600C3  E600D4  E600E5  E600F6  E600G7  E600H8  E600I9 E600J10 E600K11 
##       2       2       2       2       2       2       2       2       2 
##  E615A1  E615B2  E615C3  E615D4  E615E5  E615F6  E615G7  E615H8  E615I9 
##       4       2       2       2       2       5       2       2       2 
## E615J10 E615K11  E645A1  E645B2  E645C3  E645D4  E645E5  E645F6  E645G7 
##       2       5       2       2       2       5       5       5       2 
##  E645H8  E645I9 E645J10 E645K11  E675A1  E675B2  E675C3  E675D4  E675E5 
##       5       2       2       2       2       5       5       5       5 
##  E675F6  E675G7  E675H8  E675I9 E675J10 E675K11  E710A1  E710B2  E710C3 
##       5       6       5       5       5       5       5       5       5 
##  E710D4  E710E5  E710F6  E710G7  E710H8  E710I9 E710J10 E710K11  E745A1 
##       2       5       5       2       5       5       5       6       5 
##  E745B2  E745C3  E745D4  E745E5  E745F6  E745G7  E745H8  E745I9 E745J10 
##       6       5       6       6       6       6       5       6       6 
## E745K11  E775A1  E775B2  E775C3  E775D4  E775E5  E775F6  E775G7  E775H8 
##       6       6       5       5       6       6       6       6       6 
##  E775I9 E775J10 E775K11  E800A1  E800B2  E800C3  E800D4  E800E5  E800F6 
##       6       6       6       6       6       5       6       6       6 
##  E800G7  E800H8  E800I9 E800J10 E800K11 
##       6       6       6       6       6
vasc.env.clust.6$totss
## [1] 586.5044
vasc.env.clust.6$withinss
## [1]  5.379836 12.062646  5.258928 10.120233 18.481272 17.924102
vasc.env.clust.6$tot.withinss
## [1] 69.22702
vasc.env.clust.6$betweenss
## [1] 517.2774
vasc.env.clust.6$size
## [1] 26 31 23 38 29 29

Subset different clusters for the above two K-means clusters

vasc.5.grp.one<-vasc.env.clust.5$clust[vasc.env.clust.5$clust==1]
vasc.5.grp.one
## EA1 EB1 ED1 EE1 EF1 EG1 EI1 EK1 EC2 ED2 EE2 EG2 EA3 EC3 EE3 EF3 EG3 EI3 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## EJ3 EA4 EC4 ED4 EF4 EJ4 ED5 EF5 EH6 EA7 EC7 ED7 EC8 ED8 EF8 EG8 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1
# Corresponds to partially sunny, moist, flat sites at low elevations

vasc.5.grp.two<-vasc.env.clust.5$clust[vasc.env.clust.5$clust==2]
vasc.5.grp.two
##  E615F6 E615K11  E645D4  E645E5  E645F6  E645H8  E675B2  E675C3  E675D4 
##       2       2       2       2       2       2       2       2       2 
##  E675E5  E675F6  E675H8  E675I9 E675J10 E675K11  E710A1  E710B2  E710C3 
##       2       2       2       2       2       2       2       2       2 
##  E710E5  E710F6  E710H8  E710I9 E710J10  E745A1  E745C3  E745H8  E775B2 
##       2       2       2       2       2       2       2       2       2 
##  E775C3  E800C3 
##       2       2
# Corresponds to shaded (from shrubs), medium-dry, sloped sites with shallow rooting distance to impenetrable layer 

vasc.5.grp.three<-vasc.env.clust.5$clust[vasc.env.clust.5$clust==3]
vasc.5.grp.three
##  E675G7 E710K11  E745B2  E745D4  E745E5  E745F6  E745G7  E745I9 E745J10 
##       3       3       3       3       3       3       3       3       3 
## E745K11  E775A1  E775D4  E775E5  E775F6  E775G7  E775H8  E775I9 E775J10 
##       3       3       3       3       3       3       3       3       3 
## E775K11  E800A1  E800B2  E800D4  E800E5  E800F6  E800G7  E800H8  E800I9 
##       3       3       3       3       3       3       3       3       3 
## E800J10 E800K11 
##       3       3
#Corresponds to open, high elevation sites with shallow rooting distance to impenetrable layer 

vasc.5.grp.four<-vasc.env.clust.5$clust[vasc.env.clust.5$clust==4]
vasc.5.grp.four
##    EC1    EH1    EJ1    EA2    EB2    EF2    EH2    EI2    EJ2    EK2 
##      4      4      4      4      4      4      4      4      4      4 
##    EB3    ED3    EH3    EK3    EB4    EE4    EG4    EH4    EI4    EK4 
##      4      4      4      4      4      4      4      4      4      4 
##    EA5    EB5    EC5    EE5    EG5    EH5    EI5    EJ5    EK5    EA6 
##      4      4      4      4      4      4      4      4      4      4 
##    EB6    EC6    ED6    EE6    EF6    EI6    EJ6    EK6    EB7    EE7 
##      4      4      4      4      4      4      4      4      4      4 
##    EF7    EG7    EH7    EJ7    EK7    EA8    EB8    EE8    EH8    EI8 
##      4      4      4      4      4      4      4      4      4      4 
## E600A1 E600B2 E615A1 
##      4      4      4
# Corresponds to partially sunny, moist, flat sites at low elevations
# Sites low on grid similar to group 2, but K-means is clearly distinguishing the two groups

vasc.5.grp.five<-vasc.env.clust.5$clust[vasc.env.clust.5$clust==5]
vasc.5.grp.five
##     EG6     EI7     EJ8     EK8  E600C3  E600D4  E600E5  E600F6  E600G7 
##       5       5       5       5       5       5       5       5       5 
##  E600H8  E600I9 E600J10 E600K11  E615B2  E615C3  E615D4  E615E5  E615G7 
##       5       5       5       5       5       5       5       5       5 
##  E615H8  E615I9 E615J10  E645A1  E645B2  E645C3  E645G7  E645I9 E645J10 
##       5       5       5       5       5       5       5       5       5 
## E645K11  E675A1  E710D4  E710G7 
##       5       5       5       5
# Corresponds to medium elevation, sloped, medium soil moisture sites 

# Plot the following clusters on the grid showing the general locations

Position of five K-means clusters on large grid

#make dataframe for large transects; with plot number, Habitat.description, latitude, longitude
site.trunc<-site[,c(3:5)]

#dev.new(width=11.8, height=8)
myColorRamp <- function(colors, values) {
    v <- (values - min(values))/diff(range(values))
    x <- colorRamp(colors)(v)
    rgb(x[,1], x[,2], x[,3], maxColorValue = 255)
}

cols <- myColorRamp(c("black", "gray25", "gray50", "gray85", "white"), as.numeric(vasc.env.clust.5$cluster))
colfunc <- colorRampPalette(c("black", "gray25", "gray50", "gray85", "white"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, pch=21, col="black",cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Locations of different groups \nbased on K-means clustering", cex.main=1.75)
xl <- 1
yb <- -2.3
xr <- 1.5
yt <- 2.7

legend("bottomright", c("Group 1","Group 2","Group 3","Group 4","Group 5"), 
    col=c("black", "gray25", "gray50", "gray85", "white"),  cex=1.1, pch = 16, bty="n")

Sorensen’s local contribution to beta diversity for vascular plants

#First make dataframe with longitude and latitude coordinates
site.xy<-site[,c(4:5)]
res <- beta.div(abd.sp.allsp, "hellinger", nperm=999)

#Get location LCBD scores
loc.LCBD <-res$LCBD


#subset based on values from different clusters
loc.LCBD.five.groups<- cbind(loc.LCBD, vasc.env.clust.5$cluster)
colnames(loc.LCBD.five.groups)<-c("LCBD", "Group")
loc.LCBD.five.groups.df<-as.data.frame(loc.LCBD.five.groups)

beta.div.LCBD.five.groups.mean<-aggregate(LCBD~Group,loc.LCBD.five.groups.df, mean)

boxplot(LCBD~Group, data=loc.LCBD.five.groups.df, col="grey60", main="Local contribution to beta diversity by group",
 ylab="Local contribution to beta diversity", xlab="Group", cex.lab=1.25, cex.main=1.5)

# ANOVA with a linear model for Sorensen's beta diversity
loc.LCBD.five.groups.lm <- lm(LCBD~Group, data=loc.LCBD.five.groups.df)
summary(loc.LCBD.five.groups.lm)
## 
## Call:
## lm(formula = LCBD ~ Group, data = loc.LCBD.five.groups.df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0022685 -0.0007427 -0.0001980  0.0006524  0.0028861 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.760e-03  1.976e-04  29.158   <2e-16 ***
## Group       -2.529e-05  5.811e-05  -0.435    0.664    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.001072 on 174 degrees of freedom
## Multiple R-squared:  0.001087,   Adjusted R-squared:  -0.004654 
## F-statistic: 0.1893 on 1 and 174 DF,  p-value: 0.664
# check diagnostic plots
opar <- par(mfrow=c(2,2))
plot(loc.LCBD.five.groups.lm, pch=19, col="black")

par(opar)

# Shapiro test
shapiro.test(resid(loc.LCBD.five.groups.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(loc.LCBD.five.groups.lm)
## W = 0.9665, p-value = 0.0003057
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(loc.LCBD.five.groups.df$LCBD, loc.LCBD.five.groups.df$Group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  loc.LCBD.five.groups.df$LCBD and loc.LCBD.five.groups.df$Group
## Bartlett's K-squared = 1.7332, df = 4, p-value = 0.7847

Sorensen’s local contribution to beta diversity for angiosperms

#Analyze angiosperm only data
#subset transdata to remove angiosperms
angio.sp.abd<-abd.sp.allsp[,-match(c("LycAno", "HupApr", "DipCom", "SelSel", "CysMon", "DryExp", "EquSyl", "EquArv", "EquSci", "EquVar", "JunCom", "PicMar", "PicGla", "AbiBal", "LarLar"), names(abd.sp.allsp))]

#First make dataframe with longitude and latitude coordinates
site.xy<-site[,c(4:5)]
res.bd.angio <- beta.div(angio.sp.abd, "hellinger", nperm=999)

#Get location LCBD scores
loc.LCBD.angio <-res.bd.angio$LCBD

#subset based on values from different clusters
loc.LCBD.five.groups.angio.bd<- cbind(loc.LCBD.five.groups, loc.LCBD.angio)
colnames(loc.LCBD.five.groups.angio.bd)<-c("LCBD", "Group", "LCBD.angios")
loc.LCBD.five.groups.angio.bd.df<-as.data.frame(loc.LCBD.five.groups.angio.bd)

loc.LCBD.five.groups.angio.bd.mean<-aggregate(LCBD.angios~Group,loc.LCBD.five.groups.angio.bd.df, mean)

boxplot(LCBD.angios~Group, data=loc.LCBD.five.groups.angio.bd.df, col="grey60", main="Local contribution to beta diversity by group \nfor angiosperms only",
 ylab="Local contribution to beta diversity", xlab="Group", cex.lab=1.25, cex.main=1.5)

# ANOVA with a linear model for Sorensen's beta diversity
loc.LCBD.five.groups.angios.lm <- lm(LCBD.angios~Group, data=loc.LCBD.five.groups.angio.bd.df)
summary(loc.LCBD.five.groups.angios.lm)
## 
## Call:
## lm(formula = LCBD.angios ~ Group, data = loc.LCBD.five.groups.angio.bd.df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0027416 -0.0008203 -0.0001980  0.0006086  0.0030826 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.660e-03  2.169e-04  26.102   <2e-16 ***
## Group       6.945e-06  6.379e-05   0.109    0.913    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.001176 on 174 degrees of freedom
## Multiple R-squared:  6.812e-05,  Adjusted R-squared:  -0.005679 
## F-statistic: 0.01185 on 1 and 174 DF,  p-value: 0.9134
# check diagnostic plots
opar <- par(mfrow=c(2,2))
plot(loc.LCBD.five.groups.angios.lm, pch=19, col="black")

par(opar)

# Shapiro test
shapiro.test(resid(loc.LCBD.five.groups.angios.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(loc.LCBD.five.groups.angios.lm)
## W = 0.9695, p-value = 0.0006607
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(loc.LCBD.five.groups.angio.bd.df$LCBD, loc.LCBD.five.groups.angio.bd.df$Group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  loc.LCBD.five.groups.angio.bd.df$LCBD and loc.LCBD.five.groups.angio.bd.df$Group
## Bartlett's K-squared = 1.7332, df = 4, p-value = 0.7847

Sorensen’s phylogenetic beta diversity calculations for vascular plants

# Calculate Sorensen's phylogenetic beta diversity LCBD for vasculars and angiosperms separately
vasc.ps<-phylosor(abd.sp.allsp, trans.one.tree)
angio.ps<-phylosor(angio.sp.abd, trans.one.tree)

LCBD.vasc.ps <- LCBD.comp(vasc.ps, sqrt.x=TRUE)
LCBD.vasc.ps.vals <-LCBD.vasc.ps$LCBD

LCBD.angio.ps <-LCBD.comp(angio.ps, sqrt.x=TRUE)
LCBD.angio.ps.vals <-LCBD.angio.ps$LCBD

#subset based on values from different clusters
loc.LCBD.five.groups.pbd<- cbind(loc.LCBD.five.groups.angio.bd.df, LCBD.vasc.ps.vals,LCBD.angio.ps.vals )
colnames(loc.LCBD.five.groups.pbd)<-c("LCBD", "Group", "LCBD.angios", "LCBD.vasc.ps", "LCBD.angios.ps")
loc.LCBD.five.groups.pbd.df<-as.data.frame(loc.LCBD.five.groups.pbd)

loc.LCBD.five.groups.pbd.df.mean<-aggregate(LCBD.vasc.ps~Group,loc.LCBD.five.groups.pbd.df, mean)
loc.LCBD.five.groups.pbd.df.mean
##   Group LCBD.vasc.ps
## 1     1  0.006170717
## 2     2  0.005290057
## 3     3  0.005195965
## 4     4  0.005950608
## 5     5  0.005507058
loc.LCBD.five.groups.pbd.df.angio.pbd.mean<-aggregate(LCBD.angios.ps~Group,loc.LCBD.five.groups.pbd.df, mean)
loc.LCBD.five.groups.pbd.df.mean
##   Group LCBD.vasc.ps
## 1     1  0.006170717
## 2     2  0.005290057
## 3     3  0.005195965
## 4     4  0.005950608
## 5     5  0.005507058
# Local contribution of Sorensen's phylogenetic beta diversity to groups for vascular plants
boxplot(LCBD.vasc.ps~Group, data=loc.LCBD.five.groups.pbd.df, col="grey60", main="Local contribution to phylogenetic beta diversity by group \nfor vascular plants",
 ylab="Local contribution to phylogenetic beta diversity", xlab="Group", cex.lab=1.25, cex.main=1.5)

# ANOVA with a linear model for Sorensen's beta diversity
loc.LCBD.five.groups.vasc.ps.lm <- lm(LCBD.vasc.ps~Group, data=loc.LCBD.five.groups.pbd.df)
summary(loc.LCBD.five.groups.vasc.ps.lm)
## 
## Call:
## lm(formula = LCBD.vasc.ps ~ Group, data = loc.LCBD.five.groups.pbd.df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0031434 -0.0006226  0.0001237  0.0006099  0.0017013 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.850e-03  1.580e-04  37.035   <2e-16 ***
## Group       -5.431e-05  4.647e-05  -1.169    0.244    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.000857 on 174 degrees of freedom
## Multiple R-squared:  0.007789,   Adjusted R-squared:  0.002087 
## F-statistic: 1.366 on 1 and 174 DF,  p-value: 0.2441
# check diagnostic plots
opar <- par(mfrow=c(2,2))
plot(loc.LCBD.five.groups.vasc.ps.lm, pch=19, col="black")

par(opar)

# Shapiro test
shapiro.test(resid(loc.LCBD.five.groups.vasc.ps.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(loc.LCBD.five.groups.vasc.ps.lm)
## W = 0.9797, p-value = 0.01139
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(loc.LCBD.five.groups.pbd.df$LCBD.vasc.ps, loc.LCBD.five.groups.pbd.df$Group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  loc.LCBD.five.groups.pbd.df$LCBD.vasc.ps and loc.LCBD.five.groups.pbd.df$Group
## Bartlett's K-squared = 12.6517, df = 4, p-value = 0.01311

Local contribution of Sorensen’s phylogenetic beta diversity to groups for angiosperms

boxplot(LCBD.angios.ps~Group, data=loc.LCBD.five.groups.pbd.df, col="grey60", main="Local contribution to phylogenetic beta diversity \nby group for angiosperms only",
 ylab="Local contribution to phylogenetic beta diversity", xlab="Group", cex.lab=1.25, cex.main=1.5)

# ANOVA with a linear model for Sorensen's beta diversity
loc.LCBD.five.groups.angios.ps.lm <- lm(LCBD.angios.ps~Group, data=loc.LCBD.five.groups.pbd.df)
summary(loc.LCBD.five.groups.angios.ps.lm)
## 
## Call:
## lm(formula = LCBD.angios.ps ~ Group, data = loc.LCBD.five.groups.pbd.df)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.760e-03 -4.279e-04  3.937e-05  4.153e-04  1.580e-03 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.815e-03  1.199e-04  48.493   <2e-16 ***
## Group       -4.299e-05  3.527e-05  -1.219    0.225    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0006506 on 174 degrees of freedom
## Multiple R-squared:  0.008465,   Adjusted R-squared:  0.002767 
## F-statistic: 1.485 on 1 and 174 DF,  p-value: 0.2246
# check diagnostic plots
opar <- par(mfrow=c(2,2))
plot(loc.LCBD.five.groups.angios.ps.lm, pch=19, col="black")

par(opar)

# Shapiro test
shapiro.test(resid(loc.LCBD.five.groups.angios.ps.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(loc.LCBD.five.groups.angios.ps.lm)
## W = 0.9914, p-value = 0.3762
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(loc.LCBD.five.groups.pbd.df$LCBD.angios.ps, loc.LCBD.five.groups.pbd.df$Group)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  loc.LCBD.five.groups.pbd.df$LCBD.angios.ps and loc.LCBD.five.groups.pbd.df$Group
## Bartlett's K-squared = 2.4598, df = 4, p-value = 0.6519