Tundra-shrub transition basic phylogenetic community metric visualizations

Format data-frames so that they are compatitble with phylogeny and include all 125 species

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

names(TS.lg.ab)
## [1] "PicGla"
names(TS.sm.ab)
##  [1] "AlnVir"   "AnePar"   "AntAlp"   "AntMon"   "ArcAlp"   "AreHum"  
##  [7] "ArnAng"   "BarAlp"   "BetGla"   "BisViv"   "CalCan"   "CarBel"  
## [13] "CarBig"   "CarCap"   "CarCapit" "CarDef"   "CarGla"   "CarSci"  
## [19] "CarVag"   "CerArc"   "ChaAng"   "CopTri"   "CysMon"   "DesFle"  
## [25] "DiaLap"   "DryInt"   "ElyTra"   "EmpNig"   "HupApp"   "JunCom"  
## [31] "JunTri"   "LinBor"   "LuzSpi"   "LycAnn"   "MinBif"   "MitNud"  
## [37] "MoeMac"   "OrtSec"   "PhyCae"   "PicGla"   "PoaAlp"   "PoaArc"  
## [43] "PyrGra"   "RhiMin"   "RhoGro"   "RhoLap"   "RubArc"   "RubIda"  
## [49] "SalGla"   "SalHer"   "SalUva"   "SalVes"   "SelSel"   "SolMac"  
## [55] "SolMul"   "SteLon"   "TofPus"   "TriBor"   "TriCes"   "TriSpi"  
## [61] "VacUli"   "VacVit"   "VioAdu"   "VioRen"   "MOSS"     "LICHEN"  
## [67] "ROCK"
#Remove PicGla from small plots; use PicGla values from large plots instead
TS.sm.ab<-TS.sm.ab[,c(1:39, 41:64)]
TS.ab<-cbind(TS.sm.ab, TS.lg.ab)

#upload site data for each plot
TS.site<-read.csv("Schefferville_2013_TS_sites.csv", head=TRUE, row.names=1)

#make dataframe for large transects; with plot number, Habitat.description, latitude, longitude
TS.site.trunc<-TS.site[,c(3:5)]
head(TS.site.trunc)
##         Habitat.Description Latitude..WGS.84.datum.
## TSA1.1a        Tundra-shrub                54.90393
## TSA1.1b        Tundra-shrub                54.90392
## TSA1.2a        Tundra-shrub                54.90408
## TSA1.2b        Tundra-shrub                54.90414
## TSA1.3a        Tundra-shrub                54.90430
## TSA1.3b               Shrub                54.90435
##         Longitude..WGS.84.datum.
## TSA1.1a                -67.17024
## TSA1.1b                -67.17020
## TSA1.2a                -67.17003
## TSA1.2b                -67.16995
## TSA1.3a                -67.16973
## TSA1.3b                -67.16963
#read in 100 random Transition trees
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)

 #Find species difference in trees and TS.ab
 x<-names(TS.ab)
 y<-trans.hundred.trees[[1]]$tip.label
 setdiff(x,y)
## [1] "CarCapit" "CerArc"   "DesFle"   "HupApp"   "LycAnn"
 setdiff(y,x)
##  [1] "AbiBal" "AchMil" "AgrMer" "AmeBar" "AndPol" "AveFle" "BetMin"
##  [8] "BetPum" "CarAqu" "CarBru" "CarCat" "CarDis" "CarGyn" "CarLep"
## [15] "CarLim" "CarMag" "CarTri" "CarUtr" "CasSep" "CerAlp" "CopLap"
## [22] "CorCan" "CorTri" "DasFru" "DipCom" "DryExp" "EpiHor" "EquArv"
## [29] "EquSci" "EquSyl" "EquVar" "EriVir" "EurRad" "FraVir" "GauHis"
## [36] "GeoLiv" "HupApr" "KalPol" "LarLar" "LisCor" "LonVil" "LuzPar"
## [43] "LycAno" "MaiTri" "MonUni" "MyrGal" "PacAur" "ParKot" "PetFri"
## [50] "PicMar" "PyrAsa" "RibGla" "RubCha" "SalArc" "SalArg" "SalHum"
## [57] "SalMyr" "SalPed" "SalPla" "SchPur" "SteBor" "TriAlp" "VacCes"
## [64] "VacMyr" "VacOxy" "VibEdu"
#Change spelling of following names in TS.ab so that they match with phylogeny
#"CarCapit" "CerArc"   "DesFle"   "HupApp"   "LycAnn"
colnames(TS.ab)[15]<-"CarCat"
colnames(TS.ab)[20]<-"CerAlp"
colnames(TS.ab)[24]<-"AveFle"
colnames(TS.ab)[29]<-"HupApr"
colnames(TS.ab)[34]<-"LycAno"

#create community matrix with 0 values for "AbiBal" "AchMil" "AgrMer" "AmeBar" "AndPol" "BetMin" "BetPum" "CarAqu" "CarBru" "CarDis" "CarGyn" "CarLep" "CarLim" "CarMag" "CarTri" "CarUtr" "CasSep" "CopLap" "CorCan" "CorTri" "DasFru"
#"DipCom" "DryExp" "EpiHor" "EquArv" "EquSci" "EquSyl" "EquVar" "EriVir" "EurRad" "FraVir" "GauHis" "GeoLiv" "KalPol" "LarLar" "LisCor" "LonVil" "LuzPar" "MaiTri" "MonUni" "MyrGal" "PacAur"
#"ParKot" "PetFri" "PicMar" "PyrAsa" "RibGla" "RubCha" "SalArc" "SalArg" "SalHum" "SalMyr" "SalPed" "SalPla" "SchPur" "SteBor" "TriAlp" "VacCes" "VacMyr" "VacOxy" "VibEdu"

sp.zero.com.matrix.names<-c("AbiBal", "AchMil", "AgrMer", "AmeBar", "AndPol", "BetMin", "BetPum", "CarAqu", "CarBru", "CarDis", "CarGyn", "CarLep", "CarLim", "CarMag",
"CarTri", "CarUtr", "CasSep", "CopLap", "CorCan", "CorTri", "DasFru","DipCom", "DryExp", "EpiHor", "EquArv", "EquSci", "EquSyl", "EquVar", "EriVir", "EurRad", "FraVir",
"GauHis", "GeoLiv", "KalPol", "LarLar", "LisCor", "LonVil", "LuzPar", "MaiTri", "MonUni", "MyrGal", "PacAur","ParKot", "PetFri", "PicMar", "PyrAsa", "RibGla",
"RubCha", "SalArc", "SalArg", "SalHum", "SalMyr", "SalPed", "SalPla", "SchPur", "SteBor", "TriAlp", "VacCes", "VacMyr", "VacOxy", "VibEdu")

sp.zero.com.matrix.zeros<-matrix(0, 110, 61)

#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(TS.ab)

#add 0 values for 61 species to large grid community matrix
TS.ab.allsp<-cbind(TS.ab, sp.zero.com.matrix.zeros)
dim(TS.ab.allsp)
## [1] 110 125
#again, find difference in what is in phylogeny and plots
xx<-names(TS.ab.allsp)
yy<-trans.hundred.trees[[1]]$tip.label
xxyy<-setdiff(xx,yy)
yyxx<-setdiff(yy,xx)
xxyy
## character(0)
yyxx
## character(0)

Calculate and visualize species richness of tundra-shrub transition area

#First examine species richness
TS.sp<-TS.ab.allsp
TS.sp[TS.sp>0]<-1
TS.sp.rich<-rowSums(TS.sp)


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("white", "black"), TS.sp.rich)
colfunc <- colorRampPalette(c("white","black"))

layout(matrix(1:2,ncol=2), width = c(3,1),height = c(1,1))
par(mar=c(5.1,4.1,4.1,2.1))
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, col="black", pch=21, cex=1,  cex.lab=1.5, cex.axis=1.25, main="Species Richness", cex.main=1.75)
xl <- 1
yb <- 0
xr <- 1.5
yt <- 20

#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(0,16), las=2, cex.axis=1.25)
rect(
     xl,
     head(seq(yb,yt,(yt-yb)/30),-1),
     xr,
     tail(seq(yb,yt,(yt-yb)/30),-1),
     col=colfunc(30)
    )

#locator()


hist(TS.sp.rich, col="grey50")
shapiro.test(TS.sp.rich)  
## 
##  Shapiro-Wilk normality test
## 
## data:  TS.sp.rich
## W = 0.9404, p-value = 9.478e-05
sp.rich.mean<-mean(TS.sp.rich)#mean 
sp.rich.mean
## [1] 6.945455
sp.rich.sd<-sd(TS.sp.rich)#standard deviation 
sp.rich.sd
## [1] 4.422288

```

Presence-absence NRI visualization

#do presence-absence mean pairwise difference values
#commented out because output saved as R object - abd.sp.ses.mpd
#for (i in 1:length(phy.dist.hund)){
#TS.pa.ses.mpd<-ses.mpd(TS.ab.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
# }

#saveRDS(TS.pa.ses.mpd, file="TS.pa.ses.mpd.rds")
TS.pa.ses.mpd<-readRDS("TS.pa.ses.mpd.rds")

#do presence-absence nri values
TS.pa.nri.names<-TS.pa.ses.mpd[,6, drop=FALSE]*-1
TS.pa.nri.names.dm<-data.matrix(TS.pa.nri.names, rownames.force = NA)

TS.pa.nri.names.dm[is.na(TS.pa.nri.names.dm)] <- -1

#There are lots of NA values because of the rock plots - ho
#plot out nri per plot with color ramp############################################################
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))


bg <- myColorRamp(c("white", "black"), as.numeric(TS.pa.nri.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, col="black", pch=21, cex=1, cex.lab=1.5, cex.axis=1.25, main="Net Relatedness Index", cex.main=1.75)
xl <- 1
yb <- -1
xr <- 1.5
yt <- 2.7

#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(-1,3.0), 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)
    )

Abundance-weighted NRI visualization

#do abundance mean pairwise difference values
#commented out because output saved as R object - abd.sp.ses.mpd
#for (i in 1:length(phy.dist.hund)){
#TS.abd.ses.mpd<-ses.mpd(TS.ab.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=TRUE, runs=999, iterations=1000)
# }

#saveRDS(TS.abd.ses.mpd, file="TS.abd.ses.mpd.rds")
TS.abd.ses.mpd<-readRDS("TS.abd.ses.mpd.rds")

#calculate abundance-weighted nri values
TS.abd.nri.names<-TS.abd.ses.mpd[,6, drop=FALSE]*-1
TS.abd.nri.names.dm<-data.matrix(TS.abd.nri.names, rownames.force = NA)

#Give NA values the lower than the lowest value so they appear whitle
TS.abd.nri.names.dm[is.na(TS.abd.nri.names.dm)] <- -2.25

#plot out nri per plot with color ramp############################################################
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("white", "black"), as.numeric(TS.abd.nri.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, col="black", pch=21, cex=1, cex.lab=1.5, cex.axis=1.25, main="Abundance-weighted \nNet Relatedness Index", cex.main=1.75)
xl <- 1
yb <- -2.25
xr <- 1.5
yt <- 2.7

#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(-2.25,3.0), 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)
    )

Presence-absence NTI visualizations

#Calculate presence-absence mean nearest taxon distance values
#commented out because output saved as R object
#for (i in 1:length(phy.dist.hund)){
#TS.pa.ses.mntd<-ses.mntd(TS.ab.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
# }

#saveRDS(TS.pa.ses.mntd, file="TS.pa.ses.mntd.rds")
TS.pa.ses.mntd<-readRDS("TS.pa.ses.mntd.rds")

#calculate abundance-weighted nri values
TS.pa.nti.names<-TS.pa.ses.mntd[,6, drop=FALSE]*-1
TS.pa.nti.names.dm<-data.matrix(TS.pa.nti.names, rownames.force = NA)

#Give NA values the lower than the lowest value so they appear whitle
TS.pa.nti.names.dm[is.na(TS.pa.nti.names.dm)] <- -1.5

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("white", "black"), as.numeric(TS.pa.nti.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, col="black", pch=21, cex=1, cex.lab=1.5, cex.axis=1.25, main="Nearest Taxon Index", cex.main=1.75)
xl <- 1
yb <- -1.5
xr <- 1.5
yt <- 2.5

#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(-1.5,2.5), 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)
    )

Abundance-weighted NTI visualizations

    #Calculate presence-absence mean nearest taxon distance values
#commented out because output saved as R object
#for (i in 1:length(phy.dist.hund)){
#TS.abd.ses.mntd<-ses.mntd(TS.ab.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=TRUE, runs=999, iterations=1000)
# }

#saveRDS(TS.abd.ses.mntd, file="TS.abd.ses.mntd.rds")
TS.abd.ses.mntd<-readRDS("TS.abd.ses.mntd.rds")

#calculate abundance-weighted nri values
TS.abd.nti.names<-TS.abd.ses.mntd[,6, drop=FALSE]*-1
TS.abd.nti.names.dm<-data.matrix(TS.abd.nti.names, rownames.force = NA)

#Give NA values the lower than the lowest value so they appear whitle
TS.abd.nti.names.dm[is.na(TS.abd.nti.names.dm)] <- -2.75

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("white", "black"), as.numeric(TS.abd.nti.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
bg=cols, col="black", pch=21, cex=1, cex.lab=1.5, cex.axis=1.25, main="Abundance-weighted \nNearest Taxon Index", cex.main=1.75)
xl <- 1
yb <- -2.75
xr <- 1.5
yt <- 2

#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(-2.75,2), 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)
    )