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