Schefferville basic net relatedness index (nri) and nearest taxon index (nti) plot information



The following analyzes show basic nri and nti values for each plot on the large sampling grid. This does not include the plots in the transition areas. The grid has 11 transects with 16 plots per transect. The plots towards the top right-hand part of the grid are closer to the top of Mont Irony, whereas the plots closer to the lake are in the lower left-hand portion of the grid.

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



The following shows the basic structure of the site data for the large sampling grid.

#make dataframe for large transects; with plot number, Habitat.description, latitude, longitude
site.trunc<-site[,c(3:5)]
head(site.trunc)
##     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
dim(site.trunc)
## [1] 176   3



The following dataframe shows the structure of the abundance data for the large sampling grid.

abd.sp<-cbind(sm.abd[,c(2:58, 60:75, 78:114)],lg.abd[,c(1:4)])
colnames(abd.sp)
##   [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" "CysMon"
##  [36] "DasFru" "AveFle" "DiaLap" "DipCom" "DryExp" "DryInt" "ElyTra"
##  [43] "EmpNig" "EpiHor" "EquArv" "EquSci" "EquSyl" "EquVar" "EriVir"
##  [50] "EurRad" "FraVir" "GauHis" "GeoLiv" "HupApr" "JunCom" "JunTri"
##  [57] "KalPol" "LinBor" "LisCor" "LonVil" "LuzPar" "LycAno" "MaiTri"
##  [64] "MinBif" "MitNud" "MoeMac" "MonUni" "MyrGal" "OrtSec" "PacAur"
##  [71] "ParKot" "PetFri" "PhyCae" "PoaArc" "PyrAsa" "PyrGra" "RhiMin"
##  [78] "RhoGro" "RhoLap" "RibGla" "RubArc" "RubCha" "RubIda" "SalArc"
##  [85] "SalArg" "SalGla" "SalHum" "SalPed" "SalPla" "SalUva" "SalVes"
##  [92] "SchPur" "SelSel" "SolMac" "SolMul" "SteBor" "SteLon" "TofPus"
##  [99] "TriAlp" "TriBor" "TriCes" "TriSpi" "VacCes" "VacMyr" "VacOxy"
## [106] "VacUli" "VacVit" "VibEdu" "VioAdu" "VioRen" "AbiBal" "LarLar"
## [113] "PicGla" "PicMar"
dim(abd.sp)
## [1] 176 114



Read in data for one hundred random trees from 30006 trees from the Bayesian posterior distribution. Calculate a phylogenetic distance matrix for each one of these trees.

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


  Eleven species were not in abundance data, but were in tree. These species are important to include in species pool because they are found in either one of the transition area plots.

#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)
colnames(abd.sp.allsp)
##   [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" "CysMon"
##  [36] "DasFru" "AveFle" "DiaLap" "DipCom" "DryExp" "DryInt" "ElyTra"
##  [43] "EmpNig" "EpiHor" "EquArv" "EquSci" "EquSyl" "EquVar" "EriVir"
##  [50] "EurRad" "FraVir" "GauHis" "GeoLiv" "HupApr" "JunCom" "JunTri"
##  [57] "KalPol" "LinBor" "LisCor" "LonVil" "LuzPar" "LycAno" "MaiTri"
##  [64] "MinBif" "MitNud" "MoeMac" "MonUni" "MyrGal" "OrtSec" "PacAur"
##  [71] "ParKot" "PetFri" "PhyCae" "PoaArc" "PyrAsa" "PyrGra" "RhiMin"
##  [78] "RhoGro" "RhoLap" "RibGla" "RubArc" "RubCha" "RubIda" "SalArc"
##  [85] "SalArg" "SalGla" "SalHum" "SalPed" "SalPla" "SalUva" "SalVes"
##  [92] "SchPur" "SelSel" "SolMac" "SolMul" "SteBor" "SteLon" "TofPus"
##  [99] "TriAlp" "TriBor" "TriCes" "TriSpi" "VacCes" "VacMyr" "VacOxy"
## [106] "VacUli" "VacVit" "VibEdu" "VioAdu" "VioRen" "AbiBal" "LarLar"
## [113] "PicGla" "PicMar" "AreHum" "BetMin" "BetPum" "CarBel" "CarGla"
## [120] "CopLap" "CorTri" "LuzSpi" "PoaAlp" "SalHer" "SalMyr"
dim(abd.sp.allsp)
## [1] 176 125
#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)

Presence-absence NRI per plot in large grid

#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)){
abd.sp.ses.mpd<-ses.mpd(abd.sp.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
 }
saveRDS(abd.sp.ses.mpd, file="abd.sp.ses.mpd.rds")
abd.sp.ses.mpd<-readRDS("abd.sp.ses.mpd.rds")
dim(abd.sp.ses.mpd)
## [1] 176   8
#do 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)

Plot NRI/plot

#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("grey99", "black"), as.numeric(abd.sp.nri.names.dm))
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="Net Relatedness Index", cex.main=1.75)
xl <- 1
yb <- -2.3
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.5,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)
    )


  Calculate NRI averaged per per sampling band

#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.05950713



The following plot shows mean NRI 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", cex.main=1.75)
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=4)

hist(nri.band, col="grey60") #right skewed

shapiro.test(nri.band)   #Shapiro test states that this is normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  nri.band
## W = 0.947, p-value = 0.4442

Abundance-weighted NRI



The following analysis show abundance-weighted NRI values for each plot on the large sampling grid.

#abundance-weighted nri
#do abundance-weighted mean pairwise difference values; comment out for loop so doesn't run again
#for (i in 1:length(phy.dist.hund)){
abd.sp.ses.mpd.abd<-ses.mpd(abd.sp.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=TRUE, runs=999, iterations=1000)
 }



Import abundance-weighted data

#saveRDS(abd.sp.ses.mpd.abd, file="abd.sp.ses.mpd.abd.rds")
abd.sp.ses.mpd.abd<-readRDS("abd.sp.ses.mpd.abd.rds")
#do abundance nri values
abd.sp.nri.names.abd<-abd.sp.ses.mpd.abd[,6, drop=FALSE]*-1
abd.sp.nri.names.dm.abd<-data.matrix(abd.sp.nri.names.abd, rownames.force = NA)
colnames(abd.sp.nri.names.dm.abd)
## [1] "mpd.obs.z"

Plots of abundance-weighted NRI data

#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("grey99", "black"), as.numeric(abd.sp.nri.names.dm.abd))
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="Abundance-weighted \nNet Relatedness Index", cex.main=1.75)
xl <- 1
yb <- -2.3
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(-3,3.0), 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)
    )

 
Calculate mean abundance-weighted NRI values per sampling band (code not shown)

Plot average abundance-weighted NRI values at each sampling band

plot(nri.band.abd, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Weighted \nNet Relatedness Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
 main="Average Weighted \nNet Relatedness Index \nat Each Sampling Band", cex.main=1.75)
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=4)



Evaluate data to see if it is normal.

hist(nri.band.abd, col="grey60") #a little left skewed

shapiro.test(nri.band.abd)   #Shapiro test states that this is normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  nri.band.abd
## W = 0.9273, p-value = 0.221

Nearest taxon index per plot



Import NTI data and prepare data for further analyzes. Code (not shown) resembles code for similar functions above.

#do presence-absence mean nearest taxon nidex
#for (i in 1:length(phy.dist.hund)){
#abd.sp.ses.mntd<-ses.mntd(abd.sp.allsp, phy.dist.hund[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
}
#saveRDS(abd.sp.ses.mntd, file="abd.sp.ses.mntd.rds")
abd.sp.ses.mntd<-readRDS("abd.sp.ses.mntd.rds")
#do presence-absence nti values
abd.sp.nti.names<-abd.sp.ses.mntd[,6, drop=FALSE]*-1
abd.sp.nti.names.dm<-data.matrix(abd.sp.nti.names, rownames.force = NA)

NTI per plot

#plot out nti 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("grey99", "black"), as.numeric(abd.sp.nti.names.dm))
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="Nearest Taxon Index", cex.main=1.75)
xl <- 1
yb <- -2.3
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.5,2.5), las=2, cex.axis=1.25)
rect(
     xl,
     head(seq(yb,yt,(yt-yb)/25),-1),
     xr,
     tail(seq(yb,yt,(yt-yb)/25),-1),
     col=colfunc(25)
    )



Prepare NTI data for sampling band analyzes (code not shown). Code similar to that above.

Plot NTI at each sampling band

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(-1.5,1.5))
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=4)



Show normality of NTI per band data

hist(nti.band, col="grey60") #normal

shapiro.test(nti.band)   #Shapiro test states that this is normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  nti.band
## W = 0.9774, p-value = 0.9395

Abundance-weighted net taxon index


  Import data. Code is similar to that above.

Plot abundance-weighted NTI

#plot out nti 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("grey99", "black"), as.numeric(abd.sp.nti.names.dm.abd))
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="Abundance-weighted \nNearest Taxon Index Index", cex.main=1.75)
xl <- 1
yb <- -3.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(-4,2.0), las=2, cex.axis=1.25)
rect(
     xl,
     head(seq(yb,yt,(yt-yb)/27),-1),
     xr,
     tail(seq(yb,yt,(yt-yb)/27),-1),
     col=colfunc(27)
    )



Get small dataframes for all of the sampling bands. Code same as above and not include.

Plot abundance-weighted NTI at each sampling band

plot(nti.band.abd, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Weighted \nNearest Taxon Index", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
 main="Average Abundance-weighted \nNearest Taxon Index at \nEach Sampling Band", cex.main=1.75, ylim=c(-1.5,1.5))
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=4)



Show histogram for abundance-weighted NTI values per sampling band.

hist(nti.band.abd, col="grey60") #normal

shapiro.test(nti.band.abd)   #Shapiro test states that this is normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  nti.band.abd
## W = 0.9355, p-value = 0.2974

Evaluate if significant differences between sampling bands



First, set up the data for NRI per sampling band analyzes

#Tukey test per band for NRI with presence-absence
nri.pa<-c(nri.800, nri.775,nri.745,nri.710,nri.675,nri.645,nri.615,nri.600,
nri.8.TAM, nri.7.TAM,nri.6.TAM,nri.5.TAM,nri.4.TAM,nri.3.TAM,nri.2.TAM,nri.1.TAM)



Specify sampling bands

band<-as.factor(c("800m", "800m", "800m","800m","800m","800m","800m","800m","800m","800m","800m",
"775m", "775m","775m","775m","775m","775m","775m","775m","775m","775m","775m",
"745m","745m","745m","745m","745m","745m","745m","745m","745m","745m","745m",
"710m", "710m","710m","710m","710m","710m","710m","710m","710m","710m","710m",
"675m","675m","675m","675m","675m","675m","675m","675m","675m","675m","675m",
"645m","645m","645m","645m","645m","645m","645m","645m","645m","645m","645m",
"615m", "615m","615m","615m","615m","615m","615m","615m","615m","615m","615m",
"600m","600m","600m","600m","600m","600m","600m","600m","600m","600m","600m",
"8","8","8","8","8","8","8","8","8","8","8",
"7", "7","7","7","7","7","7","7","7","7","7",
"6","6","6","6","6","6","6","6","6","6","6",
"5","5","5","5","5","5","5","5","5","5","5",
"4","4","4","4","4","4","4","4","4","4","4",
"3","3","3","3","3","3","3","3","3","3","3",
"2","2","2","2","2","2","2","2","2","2","2",
"1","1","1","1","1","1","1","1","1","1","1"))



Next, visualize the data

# First visualize the data
boxplot(nri.pa~ band, col="grey60")
 abline(h=0, lty=4)

#nri.pa.boxplot <- ggplot(nri.pa, aes(x=band, y=nri.pa)) + geom_boxplot()



Reorder according to median of each band NRI values. Plots with lowest NRI values are towards the lake.

nri.pa.med <- sort(tapply(nri.pa, band, median))
boxplot(nri.pa ~ factor(band, levels=names(nri.pa.med)),col="grey60")
abline(h=0, lty=4)



Show effect sizes of NRI with Plot design. Net relatedness index is generally higher towards the top of the grid.

# The best way to view the effect sizes graphically is to use plot.design()
plot.design(nri.pa ~ band, ylab = "Net Relatedness Index", cex=0.75)

ANOVA of NRI per sampling band



ANOVA shows that there are significant differences in NRI values among sampling bands.

# Run an ANOVA using aov()
aov.nri.pa <- aov(nri.pa~band)
summary(aov.nri.pa) 
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## band         15  48.19   3.213   4.029 3.47e-06 ***
## Residuals   160 127.59   0.797                     
## ---
## 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.pa.lm <- lm(nri.pa~band)
anova(anov.nri.pa.lm) # same value as top function
## Analysis of Variance Table
## 
## Response: nri.pa
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## band       15  48.192  3.2128  4.0289 3.468e-06 ***
## Residuals 160 127.591  0.7974                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Based on the following diagnostic plots, residual variances is homogenous and residuals are normally distributed.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.pa.lm, pch=19, col="black")

par(opar)

# Test assumption of normality of residuals
shapiro.test(resid(anov.nri.pa.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nri.pa.lm)
## W = 0.9914, p-value = 0.3719
# data is normal

# Test assumption of homogeneity of variance
bartlett.test(nri.pa, band)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nri.pa and band
## Bartlett's K-squared = 18.6622, df = 15, p-value = 0.2295
# variance is homogenous

Compare NRI among sampling bands with Tukey tests

# But where does this difference lie? We can do a post-hoc test:
nri.pa.tuk<-TukeyHSD(aov(anov.nri.pa.lm),ordered=T)
# or equivalently
nri.pa.tuk.aov<-TukeyHSD(aov.nri.pa,ordered=T) # The

# to make list easier to read
Tukey.ordered.nri.pa <- TukeyHSD(aov.nri.pa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nri.pa$band[which(Tukey.ordered.nri.pa$band[,4] < 0.05),]
##            diff        lwr      upr        p adj
## 615m-8 1.610890 0.28523019 2.936549 0.0039058675
## 745m-8 1.744991 0.41933174 3.070651 0.0009775686
## 800m-8 1.858119 0.53245912 3.183778 0.0002810220
## 775m-8 1.861158 0.53549836 3.186818 0.0002715173
## 645m-8 1.878530 0.55287067 3.204190 0.0002228436
## 800m-4 1.417642 0.09198186 2.743301 0.0234659251
## 775m-4 1.420681 0.09502110 2.746340 0.0228605889
## 645m-4 1.438053 0.11239341 2.763713 0.0196619143

Compare NRI between old grid verses extended grid using ANOVAS

#Prepare data for ANOVA between top of grid and grid extension
top<-rep("top", 88)
bot<-rep("bot", 88)
top.bot<-c(top, bot)
nri.pa.top.bot<-as.data.frame(cbind(nri.pa, top.bot))



Visualize the top and bottom of grid data.

#create boxplots
boxplot(nri.pa~top.bot, col="grey60", main="Net Relatedness Index", ylab="Net Relatedness Index", xlab="Top or Bottom of Grid", cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)

ANOVAS of NRI comparison for old grid verses grid extension plots



ANOVA results show that there are signficant differences in NRI values between 88 plots at top of grid and the 88 new plots of the extended grid.

# Run an ANOVA using aov()
aov.nri.pa.top.bot <- aov(nri.pa~top.bot)
summary(aov.nri.pa.top.bot) #
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## top.bot       1  17.28  17.276   18.96 2.27e-05 ***
## Residuals   174 158.51   0.911                     
## ---
## 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.pa.top.botlm <- lm(nri.pa~top.bot)
anova(anov.nri.pa.top.botlm) # same value as top function
## Analysis of Variance Table
## 
## Response: nri.pa
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## top.bot     1  17.276  17.276  18.965 2.268e-05 ***
## Residuals 174 158.507   0.911                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Diagnostic plots for this ANOVA show that data meet the assumptions of homogenous variance and normal distribution of the residuals. This is ‘well-behaved’ data.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.pa.top.botlm, pch=19, col="black")
## hat values (leverages) are all = 0.01136364
##  and there are no factor predictors; no plot no. 5

par(opar)

# Test assumption of normality of residuals     
shapiro.test(resid(anov.nri.pa.top.botlm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nri.pa.top.botlm)
## W = 0.9887, p-value = 0.1737
# data is normal

# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(nri.pa, top.bot)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nri.pa and top.bot
## Bartlett's K-squared = 0.175, df = 1, p-value = 0.6757

Abundance-weighted NRI ANOVA comparisons among sampling bands



Set-up data. Code not shown.


The following boxplots show that abundance-weighted NRI values are generally lower in the extended portion of the grid.

# First visualize the data
boxplot(nri.abd~ band, col="grey60")
 abline(h=0, lty=4)

#nri.pa.boxplot <- ggplot(nri.pa, aes(x=band, y=nri.pa)) + geom_boxplot()

# We can also reorder according to the median of each band nri.abd level
nri.abd.med <- sort(tapply(nri.abd, band, median))
boxplot(nri.abd ~ factor(band, levels=names(nri.abd.med)),col="grey60")
abline(h=0, lty=4)

# The best way to view the effect sizes graphically is to use plot.design()
plot.design(nri.abd ~ band, ylab = "Abundance-Weighted Net Relatedness Index", cex=0.75)

Compare NRI.abd per sampling band with ANOVAS and Tukey tests



ANOVA results show that there is a significant difference in abundance-weighted NRI values among sampling bands.

# Run an ANOVA using aov()
aov.nri.abd <- aov(nri.abd~band)
summary(aov.nri.abd) #
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## band         15  62.62   4.175   4.558 3.53e-07 ***
## Residuals   160 146.55   0.916                     
## ---
## 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.lm <- lm(nri.abd~band)
anova(anov.nri.abd.lm) # same value as top function
## Analysis of Variance Table
## 
## Response: nri.abd
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## band       15  62.621  4.1748  4.5579 3.529e-07 ***
## Residuals 160 146.551  0.9159                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.abd.lm, pch=19, col="black")

par(opar)

# Test assumption of normality of residuals;data is not as well behaved as before; normality is not normal
shapiro.test(resid(anov.nri.abd.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nri.abd.lm)
## W = 0.9789, p-value = 0.008895
# data is not notmal

# Test assumption of homogeneity of variance
bartlett.test(nri.abd, band)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nri.abd and band
## Bartlett's K-squared = 37.9772, df = 15, p-value = 0.0009093
# variance is not homogenous

Post-hoc Tukey tests to show where there are signficant differences in abundance-weighted NRI amongst sampling bands.

# But where does this difference lie? We can do a post-hoc test:
nri.abd.tuk<-TukeyHSD(aov(anov.nri.abd.lm),ordered=T)
# or equivalently
nri.abd.tuk.aov<-TukeyHSD(aov.nri.abd,ordered=T) # The

# to make list easier to read
Tukey.ordered.nri.abd <- TukeyHSD(aov.nri.abd,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nri.abd$band[which(Tukey.ordered.nri.abd$band[,4] < 0.05),]
##            diff        lwr      upr        p adj
## 600m-4 1.564370 0.14362168 2.985118 0.0162482738
## 745m-4 1.757708 0.33696015 3.178456 0.0029160335
## 645m-4 1.797667 0.37691909 3.218415 0.0019921698
## 615m-4 1.814924 0.39417610 3.235672 0.0016856182
## 800m-4 1.905215 0.48446689 3.325963 0.0006864901
## 775m-4 1.918787 0.49803885 3.339535 0.0005977607
## 745m-8 1.444647 0.02389931 2.865395 0.0418847670
## 645m-8 1.484606 0.06385825 2.905354 0.0308676296
## 615m-8 1.501863 0.08111525 2.922611 0.0269630565
## 800m-8 1.592154 0.17140604 3.012902 0.0128684594
## 775m-8 1.605726 0.18497800 3.026474 0.0114629175

Abundance-weighted NRI comparison between old grid and extension grid



Visualize the abundance-weighted NRI data divided into top and bottom of the grid plots.

# First visualize the data
boxplot(nri.abd~top.bot, col="grey60", main="Abundance-weighted \nNet Relatedness Index", ylab="Abundance-weighted \nNet Relatedness Index", xlab="Top or Bottom of Grid", cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)

ANOVA between top of grid and extension for abundance-weighted NRI data


  ANOVA results show that there is a significant difference in abundance-weighted NRI values between old grid plots and those plots from extended portion of grid.

# Run an ANOVA using aov()
aov.nri.abd.top.bot <- aov(nri.abd~top.bot)
summary(aov.nri.abd.top.bot) #
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## top.bot       1  49.66   49.66   54.16 7.01e-12 ***
## Residuals   174 159.52    0.92                     
## ---
## 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.top.botlm <- lm(nri.abd~top.bot)
anova(anov.nri.abd.top.botlm) # same value as top function
## Analysis of Variance Table
## 
## Response: nri.abd
##            Df  Sum Sq Mean Sq F value   Pr(>F)    
## top.bot     1  49.657  49.657  54.166 7.01e-12 ***
## Residuals 174 159.516   0.917                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Diagnostic plots suggest that data meets the assumptions of homogeneous variance of the residuals and normality of residuals.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nri.abd.top.botlm, pch=19, col="black")
## hat values (leverages) are all = 0.01136364
##  and there are no factor predictors; no plot no. 5

par(opar)

# Test assumption of normality of residuals - marginal, but alright.
shapiro.test(resid(anov.nri.abd.top.botlm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nri.abd.top.botlm)
## W = 0.9848, p-value = 0.05277
# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(nri.abd, top.bot)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nri.abd and top.bot
## Bartlett's K-squared = 1.6886, df = 1, p-value = 0.1938
# variance is homogenous

#there is also a significant difference in nri.abd between the top and bottom of grid ; assumptions of variance and normality are met

Nearest taxon index per sampling band (presence/absence data only)

#Prepare the data
nti.pa<-c(nti.800, nti.775,nti.745,nti.710,nti.675,nti.645,nti.615,nti.600,
nti.8.TAM, nti.7.TAM,nti.6.TAM,nti.5.TAM,nti.4.TAM,nti.3.TAM,nti.2.TAM,nti.1.TAM)

Visualize the data

#Visualize the data
boxplot(nti.pa~ band, col="grey60")
 abline(h=0, lty=4)

#nri.pa.boxplot <- ggplot(nri.pa, aes(x=band, y=nri.pa)) + geom_boxplot()

# We can also reorder according to the median of each band richness level
nti.pa.med <- sort(tapply(nti.pa, band, median))
boxplot(nti.pa ~ factor(band, levels=names(nti.pa.med)),col="grey60")
abline(h=0, lty=4)

# The best way to view the effect sizes graphically is to use plot.design()
plot.design(nti.pa ~ band, ylab = "Nearest Taxon Index", cex=0.75)

ANOVA comparisons for nearest taxon index values among sampling bands

# Run an ANOVA using aov()
aov.nti.pa <- aov(nti.pa~band)
summary(aov.nti.pa) #
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## band         15  30.14  2.0094    3.04 0.000249 ***
## Residuals   160 105.75  0.6609                     
## ---
## 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.pa.lm <- lm(nti.pa~band)
anova(anov.nti.pa.lm) # same value as top function
## Analysis of Variance Table
## 
## Response: nti.pa
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## band       15  30.142 2.00944  3.0403 0.0002493 ***
## Residuals 160 105.749 0.66093                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Diagnostic plots indicated that nearest taxon index values per sampling band data meet the assumptions of homogeneous variance of the residuals and normality of the residuals.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nti.pa.lm, pch=19, col="black")

par(opar)

# Test assumption of normality of residuals
shapiro.test(resid(anov.nti.pa.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nti.pa.lm)
## W = 0.9939, p-value = 0.6775
# data is normal

# Test assumption of homogeneity of variance
bartlett.test(nti.pa, band)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nti.pa and band
## Bartlett's K-squared = 15.7649, df = 15, p-value = 0.3978
# variance is homogeneous

Post-hoc Tukey tests to show where signficant differences in NTI amongst sampling bands.

# But where does this difference lie? We can do a post-hoc test:
nti.pa.tuk<-TukeyHSD(aov(anov.nti.pa.lm),ordered=T)
# or equivalently
nti.pa.tuk.aov<-TukeyHSD(aov.nti.pa,ordered=T) # The

# to make list easier to read
Tukey.ordered.nti.pa <- TukeyHSD(aov.nti.pa,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nti.pa$band[which(Tukey.ordered.nti.pa$band[,4] < 0.05),]
##               diff        lwr      upr       p adj
## 745m-4    1.208897 0.00203011 2.415765 0.049131499
## 800m-4    1.387816 0.18094839 2.594683 0.008994598
## 775m-4    1.485489 0.27862192 2.692357 0.003172471
## 800m-675m 1.270518 0.06365110 2.477386 0.028308547
## 775m-675m 1.368192 0.16132462 2.575059 0.010987602
## 775m-600m 1.246087 0.03921981 2.452954 0.035381735

Comparison of NTI values between old grid and 88 plots of grid extention

#prepare data
top<-rep("top", 88)
bot<-rep("bot", 88)
top.bot<-c(top, bot)
nti.pa.top.bot<-as.data.frame(cbind(nti.pa, top.bot))

Visualize data for comparison between plots between top and bottom of grid

boxplot(nti.pa~ top.bot, col="grey60", main="Nearest Taxon Index", ylab="Nearest Taxon Index", xlab="Top or Bottom of Grid", cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)

ANOVA comparisons between 88 plots at top of grid and the 88 plots of extended grid



ANOVA results indicate that there is not a significant difference in NTI values in plots in bottom compared to plots at top of grid.

# Run an ANOVA using aov()
aov.nti.pa.top.bot <- aov(nti.pa~top.bot)
summary(aov.nti.pa.top.bot) #
##              Df Sum Sq Mean Sq F value Pr(>F)
## top.bot       1   0.75  0.7465   0.961  0.328
## Residuals   174 135.14  0.7767
# Or, run a linear model 
anov.nti.pa.top.botlm <- lm(nti.pa~top.bot)
anova(anov.nti.pa.top.botlm) # same value as top function
## Analysis of Variance Table
## 
## Response: nti.pa
##            Df  Sum Sq Mean Sq F value Pr(>F)
## top.bot     1   0.747 0.74653  0.9612 0.3283
## Residuals 174 135.144 0.77669



Diagnostic plots indicate that NTI data meet the assumptions of homogeneous variance of the residuals and normality when comparing plots at the top and bottom of the sampling grid.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nti.pa.top.botlm, pch=19, col="black")
## hat values (leverages) are all = 0.01136364
##  and there are no factor predictors; no plot no. 5

par(opar)

# Test assumption of normality of residuals  
shapiro.test(resid(anov.nti.pa.top.botlm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nti.pa.top.botlm)
## W = 0.9871, p-value = 0.1063
# data is normal

# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(nti.pa, top.bot)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nti.pa and top.bot
## Bartlett's K-squared = 0.225, df = 1, p-value = 0.6352
#There is a significant difference in NTI values between top and bottom of grid

Abundance-weighted nearest taxon index comparisons between top of grid and 88 extension plots

#Abundance weighted NTI
nti.abd<-c(nti.800.abd, nti.775.abd,nti.745.abd,nti.710.abd,nti.675.abd,nti.645.abd,nti.615.abd,nti.600.abd,
nti.8.TAM.abd, nti.7.TAM.abd,nti.6.TAM.abd,nti.5.TAM.abd,nti.4.TAM.abd,nti.3.TAM.abd,nti.2.TAM.abd,nti.1.TAM.abd)

Visualization of nearest taxon index per sampling band data

# First visualize the data
boxplot(nti.abd ~ band, col="grey60")

# We can also reorder according to the median of each band richness level
nti.abd.med <- sort(tapply(nti.abd, band, median))
boxplot(nti.abd ~ factor(band, levels=names(nti.abd.med)),col="grey60")

# The best way to view the effect sizes graphically is to use plot.design()
plot.design(nti.abd ~ band, ylab = "Abundance-weighted nearest taxon index", cex=0.75)

ANOVA comparison of abundance-weighted nearest taxon index values between different sampling bands



ANOVA indicates that there is not a significant difference in nearest taxon index values between sampling bands.

# Run an ANOVA using aov()
aov.nti.abd <- aov(nti.abd~band)
summary(aov.nti.abd) 
##              Df Sum Sq Mean Sq F value Pr(>F)
## band         15  20.23   1.349   1.278  0.222
## Residuals   160 168.90   1.056
# Or, run a linear model (MaxAbund as a function of Diet)
anov.nti.abd.lm <- lm(nti.abd~band)
anova(anov.nti.abd.lm) # same value as top function
## Analysis of Variance Table
## 
## Response: nti.abd
##            Df  Sum Sq Mean Sq F value Pr(>F)
## band       15  20.229  1.3486  1.2775 0.2219
## Residuals 160 168.902  1.0556



Diagnostic plots show that residuals have heterogenous variance and are not normally distributed.

#diagnostic plots
# Plot for diagnostics
opar <- par(mfrow=c(2,2))
plot(anov.nti.abd.lm, pch=19, col="black")

par(opar)

# Test assumption of normality of residuals
shapiro.test(resid(anov.nti.abd.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nti.abd.lm)
## W = 0.9434, p-value = 1.891e-06
# data is not notmal

# Test assumption of homogeneity of variance
bartlett.test(nti.abd, band)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nti.abd and band
## Bartlett's K-squared = 35.8325, df = 15, p-value = 0.001869
# variance is not homogenous

Post-hoc Tukey comparison amongst sampling bands for abundance-weighted nearest taxon index

nti.abd.tuk<-TukeyHSD(aov(anov.nti.abd.lm),ordered=T)
# or equivalently
nti.abd.tuk.aov<-TukeyHSD(aov.nti.abd,ordered=T) # The

# to make list easier to read
Tukey.ordered.nti.abd <- TukeyHSD(aov.nti.abd,ordered=T)
#give only those results with p<0.05
Tukey.ordered.nti.abd$band[which(Tukey.ordered.nti.abd$band[,4] < 0.05),]
##      diff lwr upr p adj

Abundance-weighted nearest taxon index comparison between 88 plots at top of grid and 88 plots of grid extension.



First, prepare the data and visualize it.

nti.abd.top.bot<-as.data.frame(cbind(nti.abd, top.bot))
# First visualize the data
boxplot(nti.abd~ top.bot, col="grey60", main="Abundance-weighted Nearest Taxon Index", ylab="Abundance-Weighted Nearest Taxon Index", xlab="Top or Bottom of Grid", cex.lab=1.25, cex.main=1.5)
abline(h=0, lty=4)

ANOVA comparison between plots at top of grid and those 88 plots at bottom of grid.



ANOVA indicates that there is a significant difference in abundance-weighted nearest taxon index values in 88 plots at top of grid compared to plots at bottom of grid.

# Run an ANOVA using aov()
aov.nti.abd.top.bot <- aov(nti.abd~top.bot)
summary(aov.nti.abd.top.bot) #
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## top.bot       1   9.69   9.690   9.396 0.00252 **
## Residuals   174 179.44   1.031                   
## ---
## 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.top.botlm <- lm(nti.abd~top.bot)
anova(anov.nti.abd.top.botlm) # same value as top function
## Analysis of Variance Table
## 
## Response: nti.abd
##            Df Sum Sq Mean Sq F value   Pr(>F)   
## top.bot     1   9.69  9.6896  9.3958 0.002522 **
## Residuals 174 179.44  1.0313                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Diagnostic plots indicate that there is homogeneous variance of the residuals and that the residuals are not normally distributed.

#diagnostic plots
# Plot for diagnostics  - really nicely behaved data
opar <- par(mfrow=c(2,2))
plot(anov.nti.abd.top.botlm, pch=19, col="black")
## hat values (leverages) are all = 0.01136364
##  and there are no factor predictors; no plot no. 5

par(opar)

# Test assumption of normality of residuals - not normal
shapiro.test(resid(anov.nti.abd.top.botlm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anov.nti.abd.top.botlm)
## W = 0.939, p-value = 8.212e-07
# data is not notmal

# Test assumption of homogeneity of variance - this is well-behaved data
bartlett.test(nti.abd, top.bot)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  nti.abd and top.bot
## Bartlett's K-squared = 1.3827, df = 1, p-value = 0.2396