Basic NRI and NTI visualizations for spruce-moss forest transition

Spruce-moss forest data structure

#Load data files
#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
SMF.site<-read.csv("Schefferville_2013_SMF_sites.csv", head=TRUE, row.names=1)

names(SMF.lg.ab)
## [1] "AbeBal" "LarLar" "PicGla" "PicMar"
names(SMF.sm.ab)
##  [1] "AbiBal" "AlnVir" "AmeBar" "BetGla" "BetPum" "CalCan" "CarAqu"
##  [8] "CarCap" "CarDis" "CarGyn" "CarLep" "CarLim" "CarSci" "CarVag"
## [15] "CasSep" "ChaAng" "CopLap" "CopTri" "CorCan" "CorTri" "CysMon"
## [22] "DesFle" "EmpNig" "EquArv" "EquSyl" "EquVar" "GauHis" "GeoLiv"
## [29] "JunCom" "LarLar" "LinBor" "LisCor" "LuzPar" "LycAnn" "MaiTri"
## [36] "MitNud" "MoeMac" "MonUni" "OrtSec" "PetFri" "PicGla" "PicMar"
## [43] "PyrAsa" "RhoGro" "RibGla" "RubArc" "RubCha" "RubIda" "SalArc"
## [50] "SalGla" "SalPla" "SalVes" "SchPur" "SolMac" "SteBor" "TriBor"
## [57] "VacCae" "VacMyr" "VacOxy" "VacUli" "VacVit" "VibEdu" "VioAdu"
## [64] "VioRen" "MOSS"   "LICHEN" "ROCKS"  "WATER"
#Remove AbiBal, LarLar, PicMar, PicGla from small plots; use  values of these species from large plots instead
SMF.sm.ab<-SMF.sm.ab[,c(2:29,31:40, 43:64)]
SMF.ab<-cbind(SMF.sm.ab, SMF.lg.ab)

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

#make dataframe for large transects; with plot number, Habitat.description, latitude, longitude
SMF.site.trunc<-SMF.site[,c(3:5)]
head(SMF.site.trunc)
##          Habitat.Description Latitude..WGS.84.datum.
## SMFA1.1a         Spruce moss                54.90061
## SMFA1.1b         Spruce moss                54.90066
## SMFA1.2a         Spruce moss                54.90076
## SMFA1.2b         Spruce moss                54.90075
## SMFA1.3a         Spruce moss                54.90090
## SMFA1.3b         Spruce moss                54.90090
##          Longitude..WGS.84.datum.
## SMFA1.1a                -67.17426
## SMFA1.1b                -67.17425
## SMFA1.2a                -67.17412
## SMFA1.2b                -67.17393
## SMFA1.3a                -67.17387
## SMFA1.3b                -67.17380
#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.hund<- lapply(trans.hundred.trees, cophenetic)

#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(SMF.ab)
y<-trans.hundred.trees[[1]]$tip.label
xy<-setdiff(x,y)
yx<-setdiff(y,x)
xy
## [1] "DesFle" "LycAnn" "VacCae" "AbeBal"
yx
##  [1] "AbiBal" "AchMil" "AgrMer" "AndPol" "AnePar" "AntAlp" "AntMon"
##  [8] "ArcAlp" "AreHum" "ArnAng" "AveFle" "BarAlp" "BetMin" "BisViv"
## [15] "CarBel" "CarBig" "CarBru" "CarCat" "CarDef" "CarGla" "CarMag"
## [22] "CarTri" "CarUtr" "CerAlp" "DasFru" "DiaLap" "DipCom" "DryExp"
## [29] "DryInt" "ElyTra" "EpiHor" "EquSci" "EriVir" "EurRad" "FraVir"
## [36] "HupApr" "JunTri" "KalPol" "LonVil" "LuzSpi" "LycAno" "MinBif"
## [43] "MyrGal" "PacAur" "ParKot" "PhyCae" "PoaAlp" "PoaArc" "PyrGra"
## [50] "RhiMin" "RhoLap" "SalArg" "SalHer" "SalHum" "SalMyr" "SalPed"
## [57] "SalUva" "SelSel" "SolMul" "SteLon" "TofPus" "TriAlp" "TriCes"
## [64] "TriSpi" "VacCes"
#Change spelling of following names in SMF.ab so that they match with phylogeny
#"DesFle"  "LycAnn" "VacCae" "AbeBal"
colnames(SMF.ab)[21]<-"AveFle"
colnames(SMF.ab)[32]<-"LycAno"
colnames(SMF.ab)[53]<-"VacCes"
colnames(SMF.ab)[61]<-"AbiBal"

#create community matrix with 0 values for "AchMil" "AgrMer" "AndPol" "AnePar" "AntAlp" "AntMon" "ArcAlp" "AreHum" "ArnAng" "BarAlp" "BetMin" "BisViv" "CarBel"
#"CarBig" "CarBru" "CarCat" "CarDef" "CarGla" "CarMag" "CarTri" "CarUtr"
#"CerAlp" "DasFru" "DiaLap" "DipCom" "DryExp" "DryInt" "ElyTra" "EpiHor"
#"EquSci" "EriVir" "EurRad" "FraVir" "HupApr" "JunTri" "KalPol" "LonVil" "LuzSpi" "MinBif" "MyrGal" "PacAur" "ParKot"
# "PhyCae" "PoaAlp" "PoaArc" "PyrGra" "RhiMin" "RhoLap" "SalArg" "SalHer" "SalHum" "SalMyr" "SalPed" "SalUva" "SelSel"
#"SolMul" "SteLon" "TofPus" "TriAlp" "TriCes" "TriSpi"

sp.zero.com.matrix.names<-c("AchMil", "AgrMer", "AndPol", "AnePar", "AntAlp", "AntMon", "ArcAlp", "AreHum", "ArnAng", "BarAlp", "BetMin", "BisViv", "CarBel",
"CarBig", "CarBru", "CarCat", "CarDef", "CarGla", "CarMag", "CarTri", "CarUtr","CerAlp", "DasFru", "DiaLap", "DipCom", "DryExp", "DryInt", "ElyTra", "EpiHor",
"EquSci", "EriVir", "EurRad", "FraVir", "HupApr", "JunTri", "KalPol", "LonVil", "LuzSpi", "MinBif", "MyrGal", "PacAur", "ParKot",
"PhyCae", "PoaAlp", "PoaArc", "PyrGra", "RhiMin", "RhoLap", "SalArg", "SalHer", "SalHum", "SalMyr", "SalPed", "SalUva", "SelSel",
"SolMul", "SteLon", "TofPus", "TriAlp", "TriCes", "TriSpi")
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(SMF.ab)

#add 0 values for 61 species to large grid community matrix
SMF.ab.allsp<-cbind(SMF.ab, sp.zero.com.matrix.zeros)
dim(SMF.ab.allsp)
## [1] 110 125
#again, find difference in what is in phylogeny and plots
xx<-names(SMF.ab.allsp)
yy<-trans.hundred.trees[[1]]$tip.label
xxyy<-setdiff(xx,yy)
yyxx<-setdiff(yy,xx)
xxyy
## character(0)
yyxx
## character(0)
#First examine species richness
SMF.sp<-SMF.ab.allsp
SMF.sp[SMF.sp>0]<-1
SMF.sp.rich<-rowSums(SMF.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"), SMF.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(SMF.site.trunc$Latitude..WGS.84.datum., SMF.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 <- 21

#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,21), 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(SMF.sp.rich, col="grey50")
shapiro.test(SMF.sp.rich)  #richness is right skewed; there are more low values
## 
##  Shapiro-Wilk normality test
## 
## data:  SMF.sp.rich
## W = 0.9606, p-value = 0.00251
sp.rich.mean<-mean(SMF.sp.rich) #mean
sp.rich.mean
## [1] 10.45455
sp.rich.sd<-sd(SMF.sp.rich)   #standard deviation
sp.rich.sd
## [1] 3.204459

***

Spruce-moss forest net relatedness index visualization

#calculate 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)){
#SMF.pa.ses.mpd<-ses.mpd(SMF.ab.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
# }

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

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



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


cols<- myColorRamp(c("white", "black"), as.numeric(SMF.pa.nri.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(SMF.site.trunc$Latitude..WGS.84.datum., SMF.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 <- -2
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,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)
    )

***

Spruce-moss forest abundance-weighted NRI visualization

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

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

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

#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(SMF.abd.nri.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(SMF.site.trunc$Latitude..WGS.84.datum., SMF.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.25

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

***

Spruce-moss forest nearest taxon index visualization

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

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

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


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(SMF.pa.nti.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(SMF.site.trunc$Latitude..WGS.84.datum., SMF.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.75
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.75,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)
    )

***

Spruce-moss forest transition abundance-weighted nearest taxon index

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

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

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


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(SMF.abd.nti.names.dm))
colfunc <- colorRampPalette(c("white","black"))
par(mar=c(5.1,4.1,4.1,2.1))
plot(SMF.site.trunc$Latitude..WGS.84.datum., SMF.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.5
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.5,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)
    )