Schefferville - transitions; phylobetadiversity
Large Grid Transition Betadiversity
Date: February 12, 2015
R version 3.1.0
Prepare data files
#load data files; TS are fortundra-shrub transition area; SMF files are for spruce-moss forest abundance data
sm.abd<-read.csv("Schefferville_2013_abundance_TLE.csv", head=TRUE, row.names=1)
lg.abd<-read.csv("Schefferville_large_plots_TLE.csv", head=TRUE, row.names=1)
#tundra-shrub transition area
TS.sm.ab<-read.csv("Schefferville_2013_TS_abundance.csv", head=TRUE, row.names=1)
TS.med.ab<-read.csv("Schefferville_2013_TS_med_abundance.csv", head=TRUE, row.names=1)
TS.lg.ab<-read.csv("Schefferville_2013_TS_large_abundance.csv", head=TRUE, row.names=1)
#Spruce-moss forest transition area
SMF.sm.ab<-read.csv("Schefferville_2013_SMF_abundance.csv", head=TRUE, row.names=1)
SMF.med.ab<-read.csv("Schefferville_2013_SMF_med_abundance.csv", head=TRUE, row.names=1)
SMF.lg.ab<-read.csv("Schefferville_2013_SMF_large_abundance.csv", head=TRUE, row.names=1)
#upload site data for each plot
site<-read.csv("Schefferville_2013_sites.csv", head=TRUE, row.names=1)
TS.site<-read.csv("Schefferville_2013_TS_sites.csv", head=TRUE, row.names=1)
SMF.site<-read.csv("Schefferville_2013_SMF_sites.csv", head=TRUE, row.names=1)
Set-up dataframes for analyses
#make dataframe for large transects; with plot number, Habitat.description, latitude, longitude
site.trunc<-site[,c(3:5)]
abd.sp<-cbind(sm.abd[,c(2:58, 60:75, 78:114)],lg.abd[,c(1:4)])
#read in 1000 and 100 random Transition trees
#trans.thousand.trees<-read.nexus("trans.thousand.rand.nex")
trans.hundred.trees<-read.nexus("trans.hundred.rand.nex")
#calculate phy.dist for all 100 random trees
#phy.dist.thous<- lapply(trans.thousand.trees, cophenetic)
phy.dist.hund<- lapply(trans.hundred.trees, cophenetic)
#further analyzes will all be past on this random sample of 100 trees from the posterior distrubution; note, this will only account for uncertainty in
#....branch lengths and not topology
#find difference in what is in phylogeny and plots
x<-names(abd.sp)
y<-trans.hundred.trees[[1]]$tip.label
xy<-setdiff(x,y)
yx<-setdiff(y,x)
xy
## character(0)
yx
## [1] "AreHum" "BetMin" "BetPum" "CarBel" "CarGla" "CopLap" "CorTri"
## [8] "LuzSpi" "PoaAlp" "SalHer" "SalMyr"
#this is large grid plots; should calculate nri from species pool from all species together in plots
#create community matrix with 0 values for "AreHum" "BetMin" "BetPum" "CarBel" "CarGla" "CopLap" "CorTri" "LuzSpi" "PoaAlp" "SalHer" "SalMyr"
sp.zero.com.matrix.names<-c("AreHum", "BetMin", "BetPum", "CarBel", "CarGla", "CopLap", "CorTri", "LuzSpi", "PoaAlp", "SalHer", "SalMyr")
sp.zero.com.matrix.zeros<-matrix(0, 176, 11)
#add missing species names as colnames and plot numbers as rownames
colnames(sp.zero.com.matrix.zeros) <-sp.zero.com.matrix.names
rownames(sp.zero.com.matrix.zeros)<-rownames(abd.sp)
#add 0 values for 11 species to large grid community matrix
abd.sp.allsp<-cbind(abd.sp, sp.zero.com.matrix.zeros)
#again, find difference in what is in phylogeny and plots
xx<-names(abd.sp.allsp)
yy<-trans.hundred.trees[[1]]$tip.label
xxyy<-setdiff(xx,yy)
yyxx<-setdiff(yy,xx)
xxyy
## character(0)
yyxx
## character(0)
Show elevations per plot
#Trancate site data with elevation included
site.elev<-site[,c(3:6)]
myColorRamp <- function(colors, values) {
v <- (values - min(values))/diff(range(values))
x <- colorRamp(colors)(v)
rgb(x[,1], x[,2], x[,3], maxColorValue = 255)
}
layout(matrix(1:2,ncol=2), width = c(3,1),height = c(1,1))
cols <- myColorRamp(c("grey99", "black"), as.numeric(site.elev$Elevation..m.))
colfunc <- colorRampPalette(c("grey99","black"))
par(mar=c(5.1,4.1,4.1,2.1))
head(site.elev)
## Habitat.Description Latitude..WGS.84.datum. Longitude..WGS.84.datum.
## EA1 Spruce moss 54.89533 -67.18067
## EB1 Fen 54.89466 -67.17943
## EC1 Fen 54.89388 -67.17813
## ED1 Spruce moss - Fen 54.89323 -67.17700
## EE1 Fen 54.89268 -67.17585
## EF1 Spruce moss 54.89206 -67.17488
## Elevation..m.
## EA1 523
## EB1 522
## EC1 524
## ED1 532
## EE1 522
## EF1 531
plot(site.elev$Latitude..WGS.84.datum., site.elev$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
col=cols, pch=16, cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Elevation (metres)", cex.main=1.75)
xl <- 1
yb <- 522
xr <- 1.5
yt <- 806
#plot(NA,type="n",ann=FALSE,xlim=c(1,2),ylim=c(1,2),xaxt="n",yaxt="n",bty="n")
plot(NA,type="n",ann=FALSE,xlim=c(1,2),xaxt="n",bty="n", ylim=c(522,806), las=2, cex.axis=1.25)
rect(
xl,
head(seq(yb,yt,(yt-yb)/28),-1),
xr,
tail(seq(yb,yt,(yt-yb)/28),-1),
col=colfunc(28)
)
Plots showing species contribution to beta diversity
#Practice with beta.div function in R; this data includes all vascular plants; make sure to do final analyzes with only angios.
#First make dataframe with longitude and latitude co-occordinates
site.xy<-site[,c(4:5)]
res <- beta.div(abd.sp.allsp, "hellinger", nperm=999)
plot(site.xy, asp=1, type="n", xlab="Latitude", ylab="Longitude", main="Map
of transition LCBD")
points(site.xy, pch=21, col="black", bg="black", cex=250*res$LCBD)
#Get species CBD scores
sp.LCBD<-res$SCBD
#Load maximum clade credibility tree
one.tree<-read.nexus("trans.one.tree.nex")
#again, find difference in what is in phylogeny and plots
xx<-names(sp.LCBD)
yy<-one.tree$tip.label
xxyy<-setdiff(xx,yy)
yyxx<-setdiff(yy,xx)
xxyy
## character(0)
yyxx
## [1] "ComUmb" "CorAlt" "CorSto" "HydArb" "TheHum" "XimAme"
#drop extra tips from tree
trans.one.tree<-drop.tip(one.tree,c("ComUmb", "TheHum", "XimAme", "CorAlt", "CorSto", "HydArb"))
#use phytools library
library(phytools)
#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))]
#drop extra tips from tree
angio.one.tree<-drop.tip(trans.one.tree,c("LycAno", "HupApr", "DipCom", "SelSel", "CysMon", "DryExp", "EquSyl", "EquArv", "EquSci", "EquVar", "JunCom", "PicMar", "PicGla", "AbiBal", "LarLar"))
#First calculate the local contribution to betadiversity for the angiosperm only data
res.angio <- beta.div(angio.sp.abd, "hellinger", nperm=999)
#get species BD scores
angio.sp.LCBD<-res.angio$SCBD
#plot all vascular and angiosperm contMaps side-by-side
#dev.new(width=5.9, height=8)
par(mfrow=c(1,2), mar=c(4, 4, 0, 1), mai=c(4, 0, 0, 0))
#Black and white ContMap
trans.map$cols[]<-grey(seq(1,0,length.out=length(trans.map$cols)))
plot(trans.map,outline=TRUE,lwd=3, res=10, ftype="b",sig=2,legend=FALSE, fsize=0.35)
#add.color.bar(420, trans.map$cols, lims=c(0,0.106),title="Species contribution to beta diversity", digits=3, prompt=TRUE, lwd=3, outline=TRUE, fsize=0.6, ftype="b")
#Black and white ContMap for angiosperms only
angio.map$cols[]<-grey(seq(1,0,length.out=length(angio.map$cols)))
plot(angio.map,outline=TRUE,lwd=3, res=10, ftype="b",sig=2,legend=FALSE, fsize=0.35)
#add.color.bar(200, angio.map$cols, lims=c(0,0.120),title="Species contribution to beta diversity (angios.)", digits=3, prompt=TRUE, lwd=3, outline=TRUE, fsize=0.6, ftype="b")
Plot local contribution to beta diversity per plot on the large grid
#plot LCBD values per plot on the large grid for both vascular plants and angiosperms
#first do the calculations for significant percentages
res.per <- beta.div(abd.sp.allsp, "percentage", nperm=999, clock=TRUE)
## Time for computation = 47.110000 sec
res.per.angio <- beta.div(angio.sp.abd, "percentage", nperm=999, clock=TRUE)
## Time for computation = 42.190000 sec
#dev.new(width=11.8, height=8)
par(mfrow=c(2,2))
plot(site.xy, asp=1, type="n", xlab="Latitude", ylab="Longitude", main="Map of transition LCBD", cex.lab=1.25)
points(site.xy, pch=16, col="black", cex=200*res$LCBD)
plot(site.xy, asp=1, type="n", xlab="Latitude", ylab="Longitude", main="Map of angiosperm \ntransition LCBD angiosperms", cex.lab=1.25)
points(site.xy, pch=16, col="black", cex=200*res.angio$LCBD)
### Example using the transition data and the percentage difference dissimilarity
# Plot a map of the LCDB indices
# First, load the file of cartesian coordinates of the 70 mite sampling sites
signif <- which(res.per$p.LCBD <= 0.05) # Which are the significant LCDB indices?
nonsignif <- which(res.per$p.LCBD > 0.05) # Which are the non-significant LCDB indices?
plot(site.xy, asp=1, type="n", xlab="Latitude", ylab="Longitude", main="Map of transition \nLCBD (black= significant plots)", cex.lab=1.25)
points(site.xy[nonsignif,], pch=21, col="grey70", bg="grey70", cex=200*res.per$LCBD[nonsignif])
points(site.xy[signif,], pch=21, col="black", bg="black", cex=200*res.per$LCBD[signif])
# Plot a map of the LCDB indices
# First, load the file of cartesian coordinates of the 70 mite sampling sites
signif <- which(res.per.angio$p.LCBD <= 0.05) # Which are the significant LCDB indices?
nonsignif <- which(res.per.angio$p.LCBD > 0.05) # Which are the non-significant LCDB indices?
plot(site.xy, asp=1, type="n", xlab="Latitude", ylab="Longitude", main="Map of angiosperm transition \nLCBD (black= significant plots)", cex.lab=1.25)
points(site.xy[nonsignif,], pch=21, col="grey70", bg="grey70", cex=200*res.per.angio$LCBD[nonsignif])
points(site.xy[signif,], pch=21, col="black", bg="black", cex=200*res.per.angio$LCBD[signif])
Local contribution of beta diversity per sampling band
#LCBD for sampling band
vasc.800<-abd.sp.allsp[grepl("800*", rownames(abd.sp.allsp)),]
m.vasc.800<-colMeans(vasc.800)
vasc.775<-abd.sp.allsp[grepl("775*", rownames(abd.sp.allsp)),]
m.vasc.775<-colMeans(vasc.775)
vasc.745<-abd.sp.allsp[grepl("745*", rownames(abd.sp.allsp)),]
m.vasc.745<-colMeans(vasc.745)
vasc.710<-abd.sp.allsp[grepl("710*", rownames(abd.sp.allsp)),]
m.vasc.710<-colMeans(vasc.710)
vasc.710<-abd.sp.allsp[grepl("710*", rownames(abd.sp.allsp)),]
m.vasc.710<-colMeans(vasc.710)
vasc.675<-abd.sp.allsp[grepl("675*", rownames(abd.sp.allsp)),]
m.vasc.675<-colMeans(vasc.675)
vasc.645<-abd.sp.allsp[grepl("645*", rownames(abd.sp.allsp)),]
m.vasc.645<-colMeans(vasc.645)
vasc.615<-abd.sp.allsp[grepl("615*", rownames(abd.sp.allsp)),]
m.vasc.615<-colMeans(vasc.615)
vasc.600<-abd.sp.allsp[grepl("600*", rownames(abd.sp.allsp)),]
m.vasc.600<-colMeans(vasc.600)
vasc.1<-abd.sp.allsp[c(1:11),]
m.vasc.1<-colMeans(vasc.1)
vasc.2<-abd.sp.allsp[c(12:22),]
m.vasc.2<-colMeans(vasc.2)
vasc.3<-abd.sp.allsp[c(23:33),]
m.vasc.3<-colMeans(vasc.3)
vasc.4<-abd.sp.allsp[c(34:44),]
m.vasc.4<-colMeans(vasc.4)
vasc.5<-abd.sp.allsp[c(45:55),]
m.vasc.5<-colMeans(vasc.5)
vasc.6<-abd.sp.allsp[c(56:66),]
m.vasc.6<-colMeans(vasc.6)
vasc.7<-abd.sp.allsp[c(67:77),]
m.vasc.7<-colMeans(vasc.7)
vasc.8<-abd.sp.allsp[c(78:88),]
m.vasc.8<-colMeans(vasc.8)
res.vasc.elev.means<-rbind(m.vasc.800, m.vasc.775, m.vasc.745, m.vasc.710, m.vasc.675, m.vasc.645, m.vasc.615, m.vasc.600, m.vasc.8, m.vasc.7, m.vasc.6, m.vasc.5, m.vasc.4, m.vasc.3, m.vasc.2, m.vasc.1)
res.elev <- beta.div(res.vasc.elev.means, "hellinger", nperm=999)
res.elev.LCBD<-res.elev$LCBD
angio.800<-angio.sp.abd[grepl("800*", rownames(angio.sp.abd)),]
m.angio.800<-colMeans(angio.800)
angio.775<-angio.sp.abd[grepl("775*", rownames(angio.sp.abd)),]
m.angio.775<-colMeans(angio.775)
angio.745<-angio.sp.abd[grepl("745*", rownames(angio.sp.abd)),]
m.angio.745<-colMeans(angio.745)
angio.710<-angio.sp.abd[grepl("710*", rownames(angio.sp.abd)),]
m.angio.710<-colMeans(angio.710)
angio.710<-angio.sp.abd[grepl("710*", rownames(angio.sp.abd)),]
m.angio.710<-colMeans(angio.710)
angio.675<-angio.sp.abd[grepl("675*", rownames(angio.sp.abd)),]
m.angio.675<-colMeans(angio.675)
angio.645<-angio.sp.abd[grepl("645*", rownames(angio.sp.abd)),]
m.angio.645<-colMeans(angio.645)
angio.615<-angio.sp.abd[grepl("615*", rownames(angio.sp.abd)),]
m.angio.615<-colMeans(angio.615)
angio.600<-angio.sp.abd[grepl("600*", rownames(angio.sp.abd)),]
m.angio.600<-colMeans(angio.600)
angio.1<-angio.sp.abd[c(1:11),]
m.angio.1<-colMeans(angio.1)
angio.2<-angio.sp.abd[c(12:22),]
m.angio.2<-colMeans(angio.2)
angio.3<-angio.sp.abd[c(23:33),]
m.angio.3<-colMeans(angio.3)
angio.4<-angio.sp.abd[c(34:44),]
m.angio.4<-colMeans(angio.4)
angio.5<-angio.sp.abd[c(45:55),]
m.angio.5<-colMeans(angio.5)
angio.6<-angio.sp.abd[c(56:66),]
m.angio.6<-colMeans(angio.6)
angio.7<-angio.sp.abd[c(67:77),]
m.angio.7<-colMeans(angio.7)
angio.8<-angio.sp.abd[c(78:88),]
m.angio.8<-colMeans(angio.8)
angio.elev.means<-rbind(m.angio.800, m.angio.775, m.angio.745, m.angio.710, m.angio.675, m.angio.645, m.angio.615, m.angio.600, m.angio.8, m.angio.7, m.angio.6, m.angio.5, m.angio.4, m.angio.3, m.angio.2, m.angio.1)
res.angio.elev <- beta.div(angio.elev.means, "hellinger", nperm=999)
res.angio.elev.LCBD<-res.angio.elev$LCBD
#dev.new(width=11.8, height=4)
par(mfrow=c(1,2))
plot(res.elev.LCBD, type="l",axes=FALSE, col="black", lwd=2, ylab="Local contribution to beta diversity", xlab="Sampling band",las=1, cex.axis=0.65, cex.lab=1,
main="Mean LCBD at each sampling band \n(vasculars)", cex.main=1, ylim=c(0.03,0.10), bty="c")
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"))
axis(2)
box(bty="l", lwd=2)
plot(res.angio.elev.LCBD, type="l",axes=FALSE, col="black", lwd=2, ylab="Local contribution to beta diversity", xlab="Sampling band",las=1, cex.axis=0.65, cex.lab=1,
main="Mean LCBD at each sampling band \n(angiosperms)", cex.main=1, ylim=c(0.03,0.11), bty="c")
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"))
axis(2)
box(bty="l", lwd=2)
