Plot of large sampling grid with latitude, longitude for x and y axes

site.trunc<-site[,c(3:5)]
plot(site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum.)

#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum.)

Plot of tundra-shrub transition area with latitude, longitude for x and y axes

#make dataframe for TS (tundra shrub) transition; with latitude, longitude for x and y axes
TS.site.trunc<-TS.site[,c(3:5)]
plot(TS.site.trunc$Latitude..WGS.84.datum., TS.site.trunc$Longitude..WGS.84.datum.)

#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum.)

Plot of spruce moss forest-fen transition area with latitude, longitude for x and y axes

#make dataframe for SMF (spruce moss forest) transitoin; with latitude, longitude for x and y axes
SMF.site.trunc<-SMF.site[,c(3:5)]
plot(SMF.site.trunc$Latitude..WGS.84.datum., SMF.site.trunc$Longitude..WGS.84.datum.)

#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum.)

Plot of entire sampling grid coloured by habitat type; latitude, longitude for x and y axes

#try polygons per habitat type
hab<-site[,c(3:5), drop=FALSE]
hab.types<-hab[,1, drop=FALSE]
fen<-hab[which(hab$Habitat.Description=="Fen"),]
smf<-hab[which(hab$Habitat.Description=="Spruce moss"),]
smf.fen<-hab[which(hab$Habitat.Description=="Spruce moss - Fen"),]

tun<-hab[which(hab$Habitat.Description=="Tundra"),]
shr<-hab[which(hab$Habitat.Description=="Shrubs"),]
tun.shr<-hab[which(hab$Habitat.Description=="Tundra - Shrubs"),]

#use convex hull function to plot habitats
find_hull <- function(df) {
    res.ch <- df[chull(hab$Latitude..WGS.84.datum., hab$Longitude..WGS.84.datum), ]
    res <- bezier(res.ch)
    res <- data.frame(x=res$Latitude..WGS.84.datum.,y=res$hab$Longitude..WGS.84.datum)
    res$z <- res$z
    res
  }

hab.ggplot<-ggplot(hab, aes(hab$Latitude..WGS.84.datum., hab$Longitude..WGS.84.datum., colour = hab$Habitat.Description)) + geom_point()
print(hab.ggplot)

***

Plot of entire sampling grid with convex hulls for habitat type; latitude, longitude for x and y axes

hulls <- ddply(hab, "hab$Habitat.Description", find_hull)
ggplot(hab,  aes(hab$Latitude..WGS.84.datum., hab$Longitude..WGS.84.datum., colour = hab$Habitat.Description,fill = hab$Habitat.Description)) +
  geom_point() + geom_polygon(data = hulls,alpha = 0.4)

***

Plot map with points showing Spruce moss forest - fen transition

par(mfrow=c(1,2),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="grey99", pch=16, cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Spruce Moss Forest - Fen Transition", cex.main=1.75)
points(fen$Latitude..WGS.84.datum., fen$Longitude..WGS.84.datum., col="grey80", pch=16, cex=1.75)
points(smf$Latitude..WGS.84.datum., smf$Longitude..WGS.84.datum., col="grey55",pch=16, cex=1.75)
points(smf.fen$Latitude..WGS.84.datum., smf.fen$Longitude..WGS.84.datum., col="black", pch=16, cex=1.75)

***

Plot map with points showing tundra-shrub transition

plot(site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum., xlab="Latitude", ylab="Longitude",
col="grey99", pch=16, cex=1.75, cex.lab=1.5, cex.axis=1.25, main="Tundra - Shrubs Transition", cex.main=1.75)
points(tun$Latitude..WGS.84.datum., tun$Longitude..WGS.84.datum., col="grey80", pch=16, cex=1.75)
points(shr$Latitude..WGS.84.datum., shr$Longitude..WGS.84.datum., col="grey55",pch=16, cex=1.75)
points(tun.shr$Latitude..WGS.84.datum., tun.shr$Longitude..WGS.84.datum., col="black", pch=16, cex=1.75)

***

Plot species richness 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))

abd.sp<-cbind(sm.abd[,c(2:58, 60:75, 78:114)],lg.abd[,c(1:4)])
abd.sp.abd<-cbind(sm.abd[,c(2:58, 60:75, 78:114)],lg.abd[,c(1:4)])
abd.sp[abd.sp>0]<-1
abd.sp.rich<-rowSums(abd.sp)
cols <- myColorRamp(c("grey99", "black"), abd.sp.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", 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)
    )

***

Richness - Histogram, shapiro test, mean, sd

hist(abd.sp.rich, col="grey50")

shapiro.test(abd.sp.rich)  #richness is right skewed; there are more low values
## 
##  Shapiro-Wilk normality test
## 
## data:  abd.sp.rich
## W = 0.9455, p-value = 2.877e-06
#Mean species richness for plots in large grid
sp.rich.mean<-mean(abd.sp.rich)
sp.rich.mean
## [1] 10.28409
#SD of species richness for plots in large grid
sp.rich.sd<-sd(abd.sp.rich)  
sp.rich.sd
## [1] 4.443805

Plot average species richness at each sampling band

plot(rich.band, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Species Richness", xlab="Sampling Band",las=1, cex.axis=1.25, cex.lab=1.5,
 main="Average Species Richness at Each Sampling Band", cex.main=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)

***

Plot histogram of underlying data for species richness at each sampling band

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

shapiro.test(rich.band)   #Shapiro test states that not normal data
## 
##  Shapiro-Wilk normality test
## 
## data:  rich.band
## W = 0.876, p-value = 0.03362

Boxplot of richness per elevation band

boxplot(richness ~ band, col="grey60")

***

Reorder based on the median of each band richness level

rich.med <- sort(tapply(richness, band, median))
boxplot(richness ~ factor(band, levels=names(rich.med)),col="grey60")

***

Visualize effect sizes with plot.design

# The best way to view the effect sizes graphically is to use plot.design()
plot.design(richness ~ band, ylab = "Species Richness", cex=0.75)

***

ANOVAs for richness based on band

# Run an ANOVA using aov()
aov.rich <- aov(sqrt(richness)~band)
summary(aov.rich) #
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## band         15  35.57  2.3713   8.351 7.22e-14 ***
## Residuals   160  45.43  0.2839                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Or, run a linear model (MaxAbund as a function of Diet)
anovrich.lm <- lm(sqrt(richness)~band)
anova(anovrich.lm) 
## Analysis of Variance Table
## 
## Response: sqrt(richness)
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## band       15 35.570 2.37132  8.3514 7.22e-14 ***
## Residuals 160 45.431 0.28394                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Show diagnostic plots for Richness per band anovas

opar <- par(mfrow=c(2,2))
plot(anovrich.lm, pch=19, col="black")

par(opar)

# Test assumption of normality of residuals
shapiro.test(resid(anovrich.lm))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(anovrich.lm)
## W = 0.9911, p-value = 0.3482
# data is normal

# Test assumption of homogeneity of variance
bartlett.test(sqrt(richness), band)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  sqrt(richness) and band
## Bartlett's K-squared = 24.2141, df = 15, p-value = 0.06155
# variance is not homogenous

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

# to make list easier to read
Tukey.ordered <- TukeyHSD(aov.rich,ordered=T)
#give only those results with p<0.05
Tukey.ordered$band[which(Tukey.ordered$band[,4] < 0.05),]
##             diff         lwr      upr        p adj
## 5-745m 0.7911863 0.000145557 1.582227 4.990434e-02
## 4-745m 0.9088356 0.117794848 1.699876 9.109216e-03
## 3-745m 0.9684484 0.177407656 1.759489 3.463054e-03
## 1-745m 1.3570851 0.566044369 2.148126 1.716598e-06
## 2-745m 1.6108627 0.819821966 2.401903 4.869711e-09
## 4-775m 0.8270147 0.035973952 1.618055 3.068738e-02
## 3-775m 0.8866275 0.095586760 1.677668 1.283876e-02
## 1-775m 1.2752642 0.484223473 2.066305 9.976517e-06
## 2-775m 1.5290418 0.738001070 2.320083 3.427531e-08
## 4-675m 0.8086868 0.017646020 1.599728 3.949905e-02
## 3-675m 0.8682996 0.077258828 1.659340 1.691645e-02
## 1-675m 1.2569363 0.465895541 2.047977 1.464862e-05
## 2-675m 1.5107139 0.719673138 2.301755 5.267555e-08
## 1-615m 1.1729398 0.381899030 1.963981 8.092737e-05
## 2-615m 1.4267174 0.635676627 2.217758 3.636729e-07
## 1-800m 1.0493627 0.258321985 1.840403 8.447851e-04
## 2-800m 1.3031403 0.512099583 2.094181 5.522128e-06
## 1-645m 1.0309210 0.239880235 1.821962 1.176090e-03
## 2-645m 1.2846986 0.493657832 2.075739 8.174547e-06
## 1-710m 1.0292823 0.238241557 1.820323 1.210867e-03
## 2-710m 1.2830599 0.492019154 2.074101 8.462952e-06
## 1-600m 0.9538715 0.162830758 1.744912 4.413182e-03
## 2-600m 1.2076491 0.416608355 1.998690 4.035659e-05
## 1-7    0.9213018 0.130261033 1.712343 7.481759e-03
## 2-7    1.1750794 0.384038631 1.966120 7.756401e-05
## 1-6    0.7935398 0.002499014 1.584581 4.837969e-02
## 2-6    1.0473174 0.256276612 1.838358 8.765900e-04
## 2-8    0.8344147 0.043373907 1.625455 2.765498e-02
## 2-5    0.8196764 0.028635665 1.610717 3.398194e-02

NRI based on presence-absence values

for (i in 1:length(phy.dist)){
abd.sp.ses.mpd<-ses.mpd(abd.sp, phy.dist[[i]], null.model="richness", abundance.weighted=FALSE, runs=999, iterations=1000)
 }
#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 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)
    )

***