Schefferville Transitions NRI without Ferns and Gymnosperms
Tammy Elliott
Thursday, January 22, 2015
Schefferville Transitions NRI/NTI without Ferns and Gymnosperms
The following analyzes calculate abundance-weighted and non-abundance weighted NRI and NTI values for the 176 plots of the large sampling grid on Mont Irony without ferns and gymnosperms. The dataset has 125 species in total, including 5 gymnosperm species and 10 fern/lycophye species.
Calculate presence/absence NRI without ferns
abd.sp<-cbind(sm.abd[,c(2:58, 60:75, 78:114)],lg.abd[,c(1:4)])
#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)
#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 topologyFind difference in species in phylogeny and plots
#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)Species included in analysis without ferns.
abd.sp.no.ferns<-abd.sp.allsp[,c(1:34,36:38, 41:44,49:53, 55:61, 63:92, 94:125)]
names(abd.sp.no.ferns)## [1] "AchMil" "AgrMer" "AlnVir" "AmeBar" "AndPol" "AnePar" "AntAlp"
## [8] "AntMon" "ArcAlp" "ArnAng" "BarAlp" "BetGla" "BisViv" "CalCan"
## [15] "CarAqu" "CarBig" "CarBru" "CarCap" "CarCat" "CarDef" "CarDis"
## [22] "CarGyn" "CarLep" "CarLim" "CarMag" "CarSci" "CarTri" "CarUtr"
## [29] "CarVag" "CasSep" "CerAlp" "ChaAng" "CopTri" "CorCan" "DasFru"
## [36] "AveFle" "DiaLap" "DryInt" "ElyTra" "EmpNig" "EpiHor" "EriVir"
## [43] "EurRad" "FraVir" "GauHis" "GeoLiv" "JunCom" "JunTri" "KalPol"
## [50] "LinBor" "LisCor" "LonVil" "LuzPar" "MaiTri" "MinBif" "MitNud"
## [57] "MoeMac" "MonUni" "MyrGal" "OrtSec" "PacAur" "ParKot" "PetFri"
## [64] "PhyCae" "PoaArc" "PyrAsa" "PyrGra" "RhiMin" "RhoGro" "RhoLap"
## [71] "RibGla" "RubArc" "RubCha" "RubIda" "SalArc" "SalArg" "SalGla"
## [78] "SalHum" "SalPed" "SalPla" "SalUva" "SalVes" "SchPur" "SolMac"
## [85] "SolMul" "SteBor" "SteLon" "TofPus" "TriAlp" "TriBor" "TriCes"
## [92] "TriSpi" "VacCes" "VacMyr" "VacOxy" "VacUli" "VacVit" "VibEdu"
## [99] "VioAdu" "VioRen" "AbiBal" "LarLar" "PicGla" "PicMar" "AreHum"
## [106] "BetMin" "BetPum" "CarBel" "CarGla" "CopLap" "CorTri" "LuzSpi"
## [113] "PoaAlp" "SalHer" "SalMyr"Richness of plots without ferns
#Change to presence/absence matrix for richness calculations
abd.sp.no.ferns.pres<-abd.sp.no.ferns
abd.sp.no.ferns.pres[abd.sp.no.ferns.pres>0.01]<-1
abd.sp.no.ferns.rich<-rowSums(abd.sp.no.ferns.pres)
abd.sp.no.ferns.rich## EA1 EB1 EC1 ED1 EE1 EF1 EG1 EH1 EI1
## 13 19 17 18 20 7 19 11 9
## EJ1 EK1 EA2 EB2 EC2 ED2 EE2 EF2 EG2
## 17 14 12 8 13 25 16 26 15
## EH2 EI2 EJ2 EK2 EA3 EB3 EC3 ED3 EE3
## 17 15 19 19 16 8 12 18 13
## EF3 EG3 EH3 EI3 EJ3 EK3 EA4 EB4 EC4
## 9 14 6 10 13 8 10 10 12
## ED4 EE4 EF4 EG4 EH4 EI4 EJ4 EK4 EA5
## 11 6 17 18 9 8 14 10 7
## EB5 EC5 ED5 EE5 EF5 EG5 EH5 EI5 EJ5
## 9 12 11 17 13 11 8 14 6
## EK5 EA6 EB6 EC6 ED6 EE6 EF6 EG6 EH6
## 11 13 6 12 12 7 9 5 8
## EI6 EJ6 EK6 EA7 EB7 EC7 ED7 EE7 EF7
## 10 15 6 12 7 13 6 11 10
## EG7 EH7 EI7 EJ7 EK7 EA8 EB8 EC8 ED8
## 8 9 7 6 7 10 5 11 7
## EE8 EF8 EG8 EH8 EI8 EJ8 EK8 E600A1 E600B2
## 12 17 7 11 11 8 10 8 12
## E600C3 E600D4 E600E5 E600F6 E600G7 E600H8 E600I9 E600J10 E600K11
## 7 9 10 8 6 10 5 7 10
## E615A1 E615B2 E615C3 E615D4 E615E5 E615F6 E615G7 E615H8 E615I9
## 10 3 7 8 6 8 8 6 8
## E615J10 E615K11 E645A1 E645B2 E645C3 E645D4 E645E5 E645F6 E645G7
## 8 7 12 14 12 3 7 11 7
## E645H8 E645I9 E645J10 E645K11 E675A1 E675B2 E675C3 E675D4 E675E5
## 3 6 10 7 8 5 6 9 4
## E675F6 E675G7 E675H8 E675I9 E675J10 E675K11 E710A1 E710B2 E710C3
## 6 5 7 8 6 6 10 10 5
## E710D4 E710E5 E710F6 E710G7 E710H8 E710I9 E710J10 E710K11 E745A1
## 7 7 11 8 6 6 5 7 11
## E745B2 E745C3 E745D4 E745E5 E745F6 E745G7 E745H8 E745I9 E745J10
## 11 9 11 2 9 6 4 4 7
## E745K11 E775A1 E775B2 E775C3 E775D4 E775E5 E775F6 E775G7 E775H8
## 3 4 6 9 8 12 11 3 5
## E775I9 E775J10 E775K11 E800A1 E800B2 E800C3 E800D4 E800E5 E800F6
## 6 3 14 5 6 9 9 6 11
## E800G7 E800H8 E800I9 E800J10 E800K11
## 7 8 10 9 10Visualization of Species Richness/plot without ferns
#Plot species richness of plots without ferns
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"), abd.sp.no.ferns.rich)
colfunc <- colorRampPalette(c("grey99","black"))
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",
col=cols, pch=16, cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Species Richness of \nplots without ferns", cex.main=1.75)
xl <- 1
yb <- 2
xr <- 1.5
yt <- 28
#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,30), 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)
)
#locator()Histogram of frequency of species richness/plot without ferns with Shapiro test
#A histogram of species richness per plot without ferns; slightly left skewed
hist(abd.sp.no.ferns.rich, col="grey50")
shapiro.test(abd.sp.no.ferns.rich) #richness is right skewed; there are more low values##
## Shapiro-Wilk normality test
##
## data: abd.sp.no.ferns.rich
## W = 0.9337, p-value = 3.114e-07Mean richness and standard deviation of NRI plot data without ferns.
sp.rich.no.ferns.mean<-mean(abd.sp.no.ferns.rich)#mean of 9.607955 / plot
sp.rich.no.ferns.mean## [1] 9.607955sp.rich.no.ferns.sd<-sd(abd.sp.no.ferns.rich) #standard deviation of 4.206083 species/ plot
sp.rich.no.ferns.sd## [1] 4.206083Calculate NRI values/plot and get values per sampling band
#do presence-absence mean pairwise difference values for vegetation matrix without ferns.
#for (i in 1:length(phy.dist.hund)){
#abd.sp.no.ferns.ses.mpd<-ses.mpd(abd.sp.no.ferns, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
#}
#Save and load R object presence-absence mean pairwise difference values for matrix without ferns analysis.
#saveRDS(abd.sp.no.ferns.ses.mpd, file="abd.sp.no.ferns.ses.mpd.rds")
abd.sp.no.ferns.ses.mpd<-readRDS("abd.sp.no.ferns.ses.mpd.rds")
#calculate presence-absence nri values
abd.sp.no.ferns.nri.names<-abd.sp.no.ferns.ses.mpd[,6, drop=FALSE]*-1
abd.sp.no.ferns.nri.names.dm<-data.matrix(abd.sp.no.ferns.nri.names, rownames.force = NA)
#get small dataframes for each of eight elevations in original transects
nri.600.nf<-abd.sp.no.ferns.nri.names.dm[grepl("600*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.615.nf<-abd.sp.no.ferns.nri.names.dm[grepl("615*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.645.nf<-abd.sp.no.ferns.nri.names.dm[grepl("645*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.675.nf<-abd.sp.no.ferns.nri.names.dm[grepl("675*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.710.nf<-abd.sp.no.ferns.nri.names.dm[grepl("710*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.745.nf<-abd.sp.no.ferns.nri.names.dm[grepl("745*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.775.nf<-abd.sp.no.ferns.nri.names.dm[grepl("775*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.800.nf<-abd.sp.no.ferns.nri.names.dm[grepl("800*", rownames(abd.sp.no.ferns.nri.names.dm)), ]
#get small dataframes for each of eight additional lines
nri.1.nf<-abd.sp.no.ferns.nri.names.dm[grepl("1", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.1.TAM.nf<-nri.1.nf[c(1:11)]
nri.2.nf<-abd.sp.no.ferns.nri.names.dm[grepl("2", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.2.TAM.nf<-nri.2.nf[c(1:11)]
nri.3.nf<-abd.sp.no.ferns.nri.names.dm[grepl("3", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.3.TAM.nf<-nri.3.nf[c(1:11)]
nri.4.nf<-abd.sp.no.ferns.nri.names.dm[grepl("4", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.4.TAM.nf<-nri.4.nf[c(1:11)]
nri.5.nf<-abd.sp.no.ferns.nri.names.dm[grepl("5", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.5.TAM.nf<-nri.5.nf[c(1:11)]
nri.6.nf<-abd.sp.no.ferns.nri.names.dm[grepl("6", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.6.TAM.nf<-nri.6.nf[c(1:11)]
nri.7.nf<-abd.sp.no.ferns.nri.names.dm[grepl("7", rownames(abd.sp.no.ferns.nri.names.dm)), ]
nri.7.TAM.nf<-nri.7.nf[c(1:11)]
nri.8.nf<-abd.sp.no.ferns.nri.names.dm[grepl("8", rownames(abd.sp.no.ferns.nri.names.dm)),]
nri.8.TAM.nf<-nri.8.nf[c(1:11)]
#find average richness per 11 row elevation bands and similar 'elevation' bands in TAM plots
nri.600.mean.nf<-mean(nri.600.nf)
nri.615.mean.nf<-mean(nri.615.nf)
nri.645.mean.nf<-mean(nri.645.nf)
nri.675.mean.nf<-mean(nri.675.nf)
nri.710.mean.nf<-mean(nri.710.nf)
nri.745.mean.nf<-mean(nri.745.nf)
nri.775.mean.nf<-mean(nri.775.nf)
nri.800.mean.nf<-mean(nri.800.nf)
nri.1.mean.nf<-mean(nri.1.TAM.nf)
nri.2.mean.nf<-mean(nri.2.TAM.nf)
nri.3.mean.nf<-mean(nri.3.TAM.nf)
nri.4.mean.nf<-mean(nri.4.TAM.nf)
nri.5.mean.nf<-mean(nri.5.TAM.nf)
nri.6.mean.nf<-mean(nri.6.TAM.nf)
nri.7.mean.nf<-mean(nri.7.TAM.nf)
nri.8.mean.nf<-mean(nri.8.TAM.nf)
nri.band.nf<-c(nri.800.mean.nf, nri.775.mean.nf,nri.745.mean.nf,nri.710.mean.nf,nri.675.mean.nf,nri.645.mean.nf,nri.615.mean.nf,nri.600.mean.nf,
nri.8.mean.nf, nri.7.mean.nf,nri.6.mean.nf,nri.5.mean.nf,nri.4.mean.nf,nri.3.mean.nf,nri.2.mean.nf,nri.1.mean.nf)
#Bring in data for NRI with all species plots for comparison
abd.sp.ses.mpd<-readRDS("abd.sp.ses.mpd.rds")
dim(abd.sp.ses.mpd)## [1] 176 8#calculate presence-absence nri values
abd.sp.nri.names<-abd.sp.ses.mpd[,6, drop=FALSE]*-1
abd.sp.nri.names.dm<-data.matrix(abd.sp.nri.names, rownames.force = NA)
#get small dataframes for each of eight elevations in original transects
nri.600<-abd.sp.nri.names.dm[grepl("600*", rownames(abd.sp.nri.names.dm)), ]
nri.615<-abd.sp.nri.names.dm[grepl("615*", rownames(abd.sp.nri.names.dm)), ]
nri.645<-abd.sp.nri.names.dm[grepl("645*", rownames(abd.sp.nri.names.dm)), ]
nri.675<-abd.sp.nri.names.dm[grepl("675*", rownames(abd.sp.nri.names.dm)), ]
nri.710<-abd.sp.nri.names.dm[grepl("710*", rownames(abd.sp.nri.names.dm)), ]
nri.745<-abd.sp.nri.names.dm[grepl("745*", rownames(abd.sp.nri.names.dm)), ]
nri.775<-abd.sp.nri.names.dm[grepl("775*", rownames(abd.sp.nri.names.dm)), ]
nri.800<-abd.sp.nri.names.dm[grepl("800*", rownames(abd.sp.nri.names.dm)), ]
#get small dataframes for each of eight additional lines
nri.1<-abd.sp.nri.names.dm[grepl("1", rownames(abd.sp.nri.names.dm)), ]
nri.1.TAM<-nri.1[c(1:11)]
nri.2<-abd.sp.nri.names.dm[grepl("2", rownames(abd.sp.nri.names.dm)), ]
nri.2.TAM<-nri.2[c(1:11)]
nri.3<-abd.sp.nri.names.dm[grepl("3", rownames(abd.sp.nri.names.dm)), ]
nri.3.TAM<-nri.3[c(1:11)]
nri.4<-abd.sp.nri.names.dm[grepl("4", rownames(abd.sp.nri.names.dm)), ]
nri.4.TAM<-nri.4[c(1:11)]
nri.5<-abd.sp.nri.names.dm[grepl("5", rownames(abd.sp.nri.names.dm)), ]
nri.5.TAM<-nri.5[c(1:11)]
nri.6<-abd.sp.nri.names.dm[grepl("6", rownames(abd.sp.nri.names.dm)), ]
nri.6.TAM<-nri.6[c(1:11)]
nri.7<-abd.sp.nri.names.dm[grepl("7", rownames(abd.sp.nri.names.dm)), ]
nri.7.TAM<-nri.7[c(1:11)]
nri.8<-abd.sp.nri.names.dm[grepl("8", rownames(abd.sp.nri.names.dm)),]
nri.8.TAM<-nri.8[c(1:11)]
#find average richness per 11 row elevation bands and similar 'elevation' bands in TAM plots
nri.600.mean<-mean(nri.600)
nri.615.mean<-mean(nri.615)
nri.645.mean<-mean(nri.645)
nri.675.mean<-mean(nri.675)
nri.710.mean<-mean(nri.710)
nri.745.mean<-mean(nri.745)
nri.775.mean<-mean(nri.775)
nri.800.mean<-mean(nri.800)
nri.1.mean<-mean(nri.1.TAM)
nri.2.mean<-mean(nri.2.TAM)
nri.3.mean<-mean(nri.3.TAM)
nri.4.mean<-mean(nri.4.TAM)
nri.5.mean<-mean(nri.5.TAM)
nri.6.mean<-mean(nri.6.TAM)
nri.7.mean<-mean(nri.7.TAM)
nri.8.mean<-mean(nri.8.TAM)
nri.band<-c(nri.800.mean, nri.775.mean,nri.745.mean,nri.710.mean,nri.675.mean,nri.645.mean,nri.615.mean,nri.600.mean,
nri.8.mean, nri.7.mean,nri.6.mean,nri.5.mean,nri.4.mean,nri.3.mean,nri.2.mean,nri.1.mean)
head(nri.band)## [1] 1.03909559 1.04213483 0.92596821 -0.12178929 0.07426231 1.05950713NRI per sampling band
#Plot average species richness at each sampling band
plot(nri.band, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Net Relatedness Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
main="Average Net Relatedness Index \nat Each Sampling Band Without Ferns", cex.main=1.75, ylim=c(-1.5,1.5), bty="c")
points(nri.band.nf, type="l",col="grey70", lwd=2)
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"), lwd=2)
axis(2, lwd=2)
box(bty="l", lwd=4)
Frequency distribution of NRI values per sampling band as well as Shapiro Test
hist(nri.band.nf, col="grey60") #right skewed
shapiro.test(nri.band.nf) #Shapiro test states that this is normal data##
## Shapiro-Wilk normality test
##
## data: nri.band.nf
## W = 0.9584, p-value = 0.6322NRI analysis without ferns and gymnosperms
Most of the code for these analyzes are not shown, as they are similar to those above.
Richness per plot without ferns and gymnosperms.
#calculate richness per plot based on presence data without ferns and gymnosperms
abd.sp.no.ferns.gym.rich<-rowSums(abd.sp.no.ferns.gym.pres)
abd.sp.no.ferns.gym.rich## EA1 EB1 EC1 ED1 EE1 EF1 EG1 EH1 EI1
## 12 17 17 16 20 6 17 10 7
## EJ1 EK1 EA2 EB2 EC2 ED2 EE2 EF2 EG2
## 15 13 11 7 12 24 15 24 14
## EH2 EI2 EJ2 EK2 EA3 EB3 EC3 ED3 EE3
## 16 15 17 18 15 7 11 16 11
## EF3 EG3 EH3 EI3 EJ3 EK3 EA4 EB4 EC4
## 9 13 6 8 13 8 10 9 10
## ED4 EE4 EF4 EG4 EH4 EI4 EJ4 EK4 EA5
## 9 5 15 16 7 6 13 8 6
## EB5 EC5 ED5 EE5 EF5 EG5 EH5 EI5 EJ5
## 8 10 9 17 12 10 7 12 5
## EK5 EA6 EB6 EC6 ED6 EE6 EF6 EG6 EH6
## 9 13 5 11 11 6 9 4 6
## EI6 EJ6 EK6 EA7 EB7 EC7 ED7 EE7 EF7
## 9 13 6 10 6 11 5 10 10
## EG7 EH7 EI7 EJ7 EK7 EA8 EB8 EC8 ED8
## 8 8 6 5 6 8 4 10 6
## EE8 EF8 EG8 EH8 EI8 EJ8 EK8 E600A1 E600B2
## 11 15 6 10 10 7 8 8 11
## E600C3 E600D4 E600E5 E600F6 E600G7 E600H8 E600I9 E600J10 E600K11
## 6 9 9 8 5 9 4 6 10
## E615A1 E615B2 E615C3 E615D4 E615E5 E615F6 E615G7 E615H8 E615I9
## 9 3 7 8 6 8 8 5 7
## E615J10 E615K11 E645A1 E645B2 E645C3 E645D4 E645E5 E645F6 E645G7
## 8 7 10 12 12 3 7 10 7
## E645H8 E645I9 E645J10 E645K11 E675A1 E675B2 E675C3 E675D4 E675E5
## 3 6 10 7 7 5 6 8 4
## E675F6 E675G7 E675H8 E675I9 E675J10 E675K11 E710A1 E710B2 E710C3
## 6 5 6 8 5 5 10 10 5
## E710D4 E710E5 E710F6 E710G7 E710H8 E710I9 E710J10 E710K11 E745A1
## 7 6 10 8 5 6 5 7 11
## E745B2 E745C3 E745D4 E745E5 E745F6 E745G7 E745H8 E745I9 E745J10
## 11 9 11 2 9 6 3 4 5
## E745K11 E775A1 E775B2 E775C3 E775D4 E775E5 E775F6 E775G7 E775H8
## 3 4 6 9 8 12 11 3 5
## E775I9 E775J10 E775K11 E800A1 E800B2 E800C3 E800D4 E800E5 E800F6
## 6 3 13 5 6 9 9 6 11
## E800G7 E800H8 E800I9 E800J10 E800K11
## 7 8 10 9 10layout(matrix(1:2,ncol=2), width = c(3,1),height = c(1,1))
cols <- myColorRamp(c("grey99", "black"), abd.sp.no.ferns.gym.rich)
colfunc <- colorRampPalette(c("grey99","black"))
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",
col=cols, pch=16, cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Species Richness without \nFerns and Gymnosperms", cex.main=1.75)
xl <- 1
yb <- 2
xr <- 1.5
yt <- 28
#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,30), 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()Histogram and Shipiro test to show structure of richness data without ferns and gymnosperms
#Histogram and Shapiro test to show structure of data
hist(abd.sp.no.ferns.gym.rich, col="grey50")
shapiro.test(abd.sp.no.ferns.gym.rich) #richness is right skewed; there are more low values #these values are not normal##
## Shapiro-Wilk normality test
##
## data: abd.sp.no.ferns.gym.rich
## W = 0.9285, p-value = 1.257e-07Calculate mean and standard deviation of richness data without ferns and gymnosperms
sp.rich.no.ferns.gym.mean<-mean(abd.sp.no.ferns.gym.rich) #mean of 8.875 / plot
sp.rich.no.ferns.gym.mean## [1] 8.875sp.rich.no.ferns.gym.sd<-sd(abd.sp.no.ferns.gym.rich) #standard deviation of 3.923009 species/ plot
sp.rich.no.ferns.gym.sd## [1] 3.923009
The following code is not shown as it is similar to that above. This code calculates NRI values without ferns and gymnosperms per plot and per sampling band.
NRI per sampling band with ferns and gymnosperms removed
#Plot NRI at each sampling band; this graph will show comparison among plots with all species and those with ferns and gymnosperms removed
plot(nri.band, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Net Relatedness Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
main="Average Net Relatedness Index \nat Each Sampling Band", cex.main=1.75, ylim=c(-3,3))
points(nri.band.nf, type="l",col="grey40", lwd=2)
points(nri.band.nf.gym, type="l",col="grey70", lwd=2)
box(bty="l", lwd=2)
abline(h=0, lty=4)
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"), lwd=2)
axis(2, lwd=2)
legend("bottomleft", c("All vascular plants", "- ferns", "- ferns and gymnosperms"), col = c("black","grey40", "grey70" ), cex=1.1,
lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99")
box(bty="l", lwd=4)
Histogram and Shapiro test of NRI values per sampling band without ferns and gymnosperms
hist(nri.band.nf.gym, col="grey60") #left skewed
shapiro.test(nri.band.nf.gym) #Shapiro test states that this is normal data##
## Shapiro-Wilk normality test
##
## data: nri.band.nf.gym
## W = 0.9183, p-value = 0.1582# Test assumption of homogeneity of variance - this is well-behaved dataComparison of NRI values without ferns and gymnosperms
#Compare NRI, NRI (without ferns), and NRI (without ferns and gymnosperms)
nri.comp<-cbind( abd.sp.nri.names.dm, abd.sp.no.ferns.nri.names.dm, abd.sp.no.ferns.gym.nri.names.dm)
colnames(nri.comp)<-c("vasculars", "no.ferns", "no.ferns.gym")
boxplot(x=as.list(as.data.frame(nri.comp)), col="grey60", main="Net Relatedness Index", ylab="Net Relatedness Index", xlab="Taxa Considered",
cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)
###t-test comparison of NRI values per plot with and without ferns and gymnosperms
#do t test
nri.comp.melt<-melt(nri.comp)
colnames(nri.comp.melt)<-c("Plot", "Taxa.included", "NRI")
# Run an ANOVA using aov()
aov.nri.taxa <- aov(NRI~Taxa.included, data=nri.comp.melt)
summary(aov.nri.taxa) ## Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 228.0 114.00 62.16 <2e-16 ***
## Residuals 525 962.9 1.83
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Or, run a linear model (MaxAbund as a function of Diet)
anov.nri.taxa.lm <- lm(NRI~Taxa.included, data=nri.comp.melt)
anova(anov.nri.taxa.lm) # same value as top function## Analysis of Variance Table
##
## Response: NRI
## Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 227.99 113.997 62.155 < 2.2e-16 ***
## Residuals 525 962.90 1.834
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnostic plots for comparisons of NRI values with and without ferns and gymnosperms
#diagnostic plots
# Plot for diagnostics - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.taxa.lm, pch=19, col="black")
par(opar)
# Test assumption of normality of residuals - not normal
shapiro.test(resid(anov.nri.taxa.lm))##
## Shapiro-Wilk normality test
##
## data: resid(anov.nri.taxa.lm)
## W = 0.9848, p-value = 2.564e-05# data is not notmal
# Test assumption of homogeneity of variance - this is not well-behaved data
bartlett.test(NRI~Taxa.included, data=nri.comp.melt)##
## Bartlett test of homogeneity of variances
##
## data: NRI by Taxa.included
## Bartlett's K-squared = 55.3737, df = 2, p-value = 9.457e-13#assumptions of variance and normality not met; cannot reliable interpret model results# But where does this difference lie? We can do a post-hoc test:
nri.taxa.tuk<-TukeyHSD(aov(anov.nri.taxa.lm),ordered=T)
# or equivalently
nri.taxa.aov<-TukeyHSD(aov.nri.taxa,ordered=T) # The
# to make list easier to read
Tukey.ordered.nri.taxa <- TukeyHSD(aov.nri.taxa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nri.taxa$Taxa.included[which(Tukey.ordered.nri.taxa$Taxa.included[,4] < 0.05),]## diff lwr upr p adj
## no.ferns.gym-no.ferns 1.409698 1.070380 1.749017 2.710689e-10
## no.ferns.gym-vasculars 1.377680 1.038361 1.716998 2.710693e-10#All NRI values combined for the three different taxa compared do not make assumptions according to Barlett and Shapiro tests, even though
#they visually look all right; signficant differnect between no.fern.gym and other two groupingsAbundance-weighted NRI analyzes
The code for the calculation of abundance-weighted NRI values and NRI calculations per sampling bands are not shown because it resembles that above.This code includes calculations for abundance-weighted NRI values without ferns and without ferns and gymnosperms
Abundance-weighted NRI at each sampling band
par(mar=c(5.1,5.7,5.5,2.1))
plot(nri.band.abd, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Abundance-weighted \nNet Relatedness Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
main="Average Abundance-weighted \nNet Relatedness Index \nat Each Sampling Band", cex.main=1.75, ylim=c(-3,3))
points(nri.band.abd.nf, type="l",col="grey40", lwd=2)
points(nri.band.abd.nf.gym, type="l",col="grey70", lwd=2)
box(bty="l", lwd=2)
abline(h=0, lty=4)
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"), lwd=2)
axis(2, lwd=2)
legend("bottomleft", c("All vascular plants", "- ferns", "- ferns and gymnosperms"), col = c("black","grey40", "grey70" ), cex=1.1,
lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99")
box(bty="l", lwd=4)
Histogram and Shapiro test of NRI-abundance data without ferns and gymnosperms
#Show histogram of data distribution of NRI-abundance data for data without ferns and gymnosperms
hist(nri.band.abd.nf.gym, col="grey60") #right skewed
shapiro.test(nri.band.abd.nf.gym) #Shapiro test states that this is normal data##
## Shapiro-Wilk normality test
##
## data: nri.band.abd.nf.gym
## W = 0.9429, p-value = 0.3866Visualization of Abundance-weighted NRI compared among plots without gymnosperms and ferns.
#Compare NRI, NRI (without ferns), and NRI (without ferns and gymnosperms)
nri.comp.abd<-cbind( abd.sp.nri.names.dm.abd, abd.sp.no.ferns.abd.nri.names.dm, abd.sp.no.ferns.gym.abd.nri.names.dm)
colnames(nri.comp.abd)<-c("vasculars", "no.ferns", "no.ferns.gym")
nri.comp.abd.melt<-melt(nri.comp.abd)
colnames(nri.comp.abd.melt)<-c("Plot", "Taxa.included", "NRI")
boxplot(x=as.list(as.data.frame(nri.comp.abd)), col="grey60", main="Abundance-weighted Net Relatedness Index ", ylab="Abundance-weigted Net Relatedness Index", xlab="Taxa Considered",
cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)
ANOVA comparisons of abundance-weighted NRI per plot including different taxa
# Run an ANOVA using aov()
aov.nri.abd.taxa <- aov(NRI~Taxa.included, data=nri.comp.abd.melt)
summary(aov.nri.abd.taxa) ### Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 222.3 111.17 54.73 <2e-16 ***
## Residuals 525 1066.4 2.03
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Or, run a linear model (MaxAbund as a function of Diet)
anov.nri.abd.taxa.lm <- lm(NRI~Taxa.included, data=nri.comp.abd.melt)
anova(anov.nri.abd.taxa.lm) # same value as top function## Analysis of Variance Table
##
## Response: NRI
## Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 222.34 111.170 54.73 < 2.2e-16 ***
## Residuals 525 1066.41 2.031
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Diagnostic tests for abundance-weighted NRI per plots including different taxa
#diagnostic plots
# Plot for diagnostics - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.abd.taxa.lm, pch=19, col="black")
par(opar)
# Test assumption of normality of residuals - not normal
shapiro.test(resid(anov.nri.abd.taxa.lm))##
## Shapiro-Wilk normality test
##
## data: resid(anov.nri.abd.taxa.lm)
## W = 0.9897, p-value = 0.0009528# data is not notmal
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(NRI~Taxa.included, data=nri.comp.abd.melt)##
## Bartlett test of homogeneity of variances
##
## data: NRI by Taxa.included
## Bartlett's K-squared = 27.4089, df = 2, p-value = 1.117e-06 #assumptions of variance and normality not metTukey tests to show between which taxa groupings differ
# But where does this difference lie? We can do a post-hoc test:
nri.abd.taxa.tuk<-TukeyHSD(aov(anov.nri.abd.taxa.lm),ordered=T)
# or equivalently
nri.abd.taxa.aov<-TukeyHSD(aov.nri.abd.taxa,ordered=T) # The
# to make list easier to read
Tukey.ordered.nri.abd.taxa <- TukeyHSD(aov.nri.abd.taxa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nri.abd.taxa$Taxa.included[which(Tukey.ordered.nri.abd.taxa$Taxa.included[,4] < 0.05),]## diff lwr upr p adj
## no.ferns.gym-no.ferns 1.518258 1.1611664 1.875349 2.710697e-10
## no.ferns.gym-vasculars 1.166705 0.8096142 1.523797 2.712999e-10Presence/absence NTI values compared among different taxon groupings
The code to calculate NTI values to compare different taxon groupings is not shown because it resembles that above.
Plot NTI compared among sampling bands for different taxon groupings
#Plot average species richness at each sampling band
par(mar=c(5.1,5.1,5.5,2.1))
plot(nti.band, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Nearest Taxon Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
main="Average Nearest Taxon Index \nat Each Sampling Band", cex.main=1.75, ylim=c(-3,3))
points(nti.band.nf, type="l",col="grey40", lwd=2)
points(nti.band.nf.gym, type="l",col="grey70", lwd=2)
box(bty="l", lwd=2)
abline(h=0, lty=4)
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"), lwd=2)
axis(2, lwd=2)
legend("bottomleft", c("All vascular plants", "- ferns", "- ferns and gymnosperms"), col = c("black","grey40", "grey70" ), cex=1.1,
lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99")
box(bty="l", lwd=4)
Plot histogram and show Shapiro test of NTI per sampling band.
hist(nti.band.nf.gym, col="grey60") #right skewed
shapiro.test(nti.band.nf.gym) #Shapiro test states that this is almost normal##
## Shapiro-Wilk normality test
##
## data: nti.band.nf.gym
## W = 0.8683, p-value = 0.02561# Test assumption of homogeneity of variance - this is well-behaved dataCompare NTI among different taxon groupings
#Compare NRI, NRI (without ferns), and NRI (without ferns and gymnosperms)
nti.comp<-cbind( abd.sp.nti.names.dm, abd.sp.no.ferns.nti.names.dm, abd.sp.no.ferns.gym.nti.names.dm)
colnames(nti.comp)<-c("vasculars", "no.ferns", "no.ferns.gym")
boxplot(x=as.list(as.data.frame(nti.comp)), col="grey60", main="Nearest Taxon Index", ylab="Nearest Taxon Index", xlab="Taxa Considered",
cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)
t-test comparison among different taxon groupings
#do t test
nti.comp.melt<-melt(nti.comp)
colnames(nti.comp.melt)<-c("Plot", "Taxa.included", "NRI")
# Run an ANOVA using aov()
aov.nti.taxa <- aov(NRI~Taxa.included, data=nti.comp.melt)
summary(aov.nti.taxa) ### Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 7.5 3.746 3.745 0.0243 *
## Residuals 525 525.2 1.000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Or, run a linear model (MaxAbund as a function of Diet)
anov.nti.taxa.lm <- lm(NRI~Taxa.included, data=nti.comp.melt)
anova(anov.nti.taxa.lm) # same value as top function## Analysis of Variance Table
##
## Response: NRI
## Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 7.49 3.7465 3.7448 0.02427 *
## Residuals 525 525.23 1.0004
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics for NTI comparisons among different taxon groupings
#diagnostic plots
# Plot for diagnostics - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nti.taxa.lm, pch=19, col="black")
par(opar)
# Test assumption of normality of residuals - good
shapiro.test(resid(anov.nti.taxa.lm))##
## Shapiro-Wilk normality test
##
## data: resid(anov.nti.taxa.lm)
## W = 0.9971, p-value = 0.4627# data is not notmal
# Test assumption of homogeneity of variance - this is not that good
bartlett.test(NRI~Taxa.included, data=nti.comp.melt)##
## Bartlett test of homogeneity of variances
##
## data: NRI by Taxa.included
## Bartlett's K-squared = 9.9798, df = 2, p-value = 0.006806 #assumptions of variance and normality not metTukey test comparison for NTI differences amongst taxon groupings
# But where does this difference lie? We can do a post-hoc test:
nti.taxa.tuk<-TukeyHSD(aov(anov.nti.taxa.lm),ordered=T)
# or equivalently
nti.taxa.aov<-TukeyHSD(aov.nti.taxa,ordered=T) # The
# to make list easier to read
Tukey.ordered.nti.taxa <- TukeyHSD(aov.nti.taxa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nti.taxa$Taxa.included[which(Tukey.ordered.nti.taxa$Taxa.included[,4] < 0.05),]## diff lwr upr p adj
## 0.28362346 0.03301649 0.53423043 0.02191381#assumption of normality good; variance no; however, there is no difference between taxa included levelsAbundance-weighted NTI compared among taxon groupings and sampling bands
The code for the following analyzes is not shown as it is similar to that above.
Visualization of abundance-weighted NTI per sampling band with different taxon groupings
#Plot average species richness at each sampling band
par(mar=c(5.1,5.7,5.5,2.1))
plot(nti.band.abd, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Abundance-weighted \nNearest Taxon Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
main="Average Abundance-weighted \nNearest Taxon Index \nat Each Sampling Band", cex.main=1.75, ylim=c(-3,3))
points(nti.band.abd.nf, type="l",col="grey40", lwd=2)
points(nti.band.abd.nf.gym, type="l",col="grey70", lwd=2)
box(bty="l", lwd=2)
abline(h=0, lty=4)
axis(1,1:16,labels=c("800m", "775m","745m", "710m", "675m", "645m", "615m", "600m", "8", "7", "6", "5", "4", "3","2", "1"), lwd=2)
axis(2, lwd=2)
legend("bottomleft", c("All vascular plants", "- ferns", "- ferns and gymnosperms"), col = c("black","grey40", "grey70" ), cex=1.1,
lty = c(1, 1, 1), pch = c(NA, NA, NA), lwd=3 , bg = "gray99")
box(bty="l", lwd=4)
Histogram and Shapiro test of abundance-weighted NTI per sampling band with different taxon groupings
hist(nti.band.abd.nf.gym, col="grey60") #right skewed
shapiro.test(nti.band.abd.nf.gym) #Shapiro test states that this is normal data##
## Shapiro-Wilk normality test
##
## data: nti.band.abd.nf.gym
## W = 0.9267, p-value = 0.2159Comparisons of abundance-weighted NTI for different taxon groupings
#Compare NRI, NRI (without ferns), and NRI (without ferns and gymnosperms)
nti.comp.abd<-cbind( abd.sp.nti.names.dm.abd, abd.sp.no.ferns.abd.nti.names.dm, abd.sp.no.ferns.gym.abd.nti.names.dm)
colnames(nti.comp.abd)<-c("Vasculars", "No.ferns", "No.ferns.gymnosperms")
nti.comp.abd.melt<-melt(nti.comp.abd)
colnames(nti.comp.abd.melt)<-c("Plot", "Taxa.included", "NTI.abd")
boxplot(NTI.abd~Taxa.included, data=nti.comp.abd.melt, col="grey60", main="Abundance-weighted Net Relatedness Index ", ylab="Abundance- weigted Net Relatedness Index", xlab="Taxa Considered",
cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)
ANOVA comparisons of abundance-weighted NTI for different taxon groupings
# Run an ANOVA using aov()
aov.nti.abd.taxa <- aov(NTI.abd~Taxa.included, data=nti.comp.abd.melt)
summary(aov.nti.abd.taxa) ### Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 22.0 11.019 11.14 1.82e-05 ***
## Residuals 525 519.2 0.989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Or, run a linear model (MaxAbund as a function of Diet)
anov.nti.abd.taxa.lm <- lm(NTI.abd~Taxa.included, data=nti.comp.abd.melt)
anova(anov.nti.abd.taxa.lm) # same value as top function## Analysis of Variance Table
##
## Response: NTI.abd
## Df Sum Sq Mean Sq F value Pr(>F)
## Taxa.included 2 22.04 11.019 11.142 1.823e-05 ***
## Residuals 525 519.20 0.989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics for abundance-weighted NTI for different taxon groupings
#diagnostic plots
# Plot for diagnostics - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nti.abd.taxa.lm, pch=19, col="black")
par(opar)
# Test assumption of normality of residuals - not normal
shapiro.test(resid(anov.nti.abd.taxa.lm))##
## Shapiro-Wilk normality test
##
## data: resid(anov.nti.abd.taxa.lm)
## W = 0.9852, p-value = 3.367e-05# data is not notmal
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(NTI.abd~Taxa.included, data=nti.comp.abd.melt)##
## Bartlett test of homogeneity of variances
##
## data: NTI.abd by Taxa.included
## Bartlett's K-squared = 1.2663, df = 2, p-value = 0.5309 #assumptions of variance and normality not metTukey test comparisons for abundance-weighted NTI for different taxon groupings
# But where does this difference lie? We can do a post-hoc test:
nti.abd.taxa.tuk<-TukeyHSD(aov(anov.nti.abd.taxa.lm),ordered=T)
# or equivalently
nti.abd.taxa.aov<-TukeyHSD(aov.nti.abd.taxa,ordered=T) # The
# to make list easier to read
Tukey.ordered.nti.abd.taxa <- TukeyHSD(aov.nti.abd.taxa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nti.abd.taxa$Taxa.included[which(Tukey.ordered.nti.abd.taxa$Taxa.included[,4] < 0.05),]## diff lwr upr p adj
## No.ferns.gymnosperms-No.ferns 0.4560609 0.2068979 0.7052239 5.988489e-05
## No.ferns.gymnosperms-Vasculars 0.4064454 0.1572824 0.6556084 4.148471e-04