The following shows the locations of the 176 plots 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.
site.trunc<-site[,c(3:5)]
plot(site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum., lwd=2, xlab="Latitude", ylab="Longitude",col="grey40", pch=16, cex=1.5, cex.lab=1.5, cex.axis=1.25, main="Location of Grid Plots", cex.main=1.5)
#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum)
These two plots show the locations of the plots in the two intensively sampled transition areas. Each intensively sampled transition area has 11 transects with 10 plots per transect. The plot layout within each of the two transects consisted of five pairs of plots within a length of approximately 100 metres.
par(mfrow=c(1,2))
#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., lwd=2, xlab="Latitude", ylab="Longitude",col="grey40", pch=16, cex=0.5, cex.lab=1.5, cex.axis=1.25, main="Location of Tundra-Shrub Transition Plots", cex.main=0.7)
#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum)
#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., lwd=2, xlab="Latitude", ylab="Longitude",col="grey40", pch=16, cex=0.5, cex.lab=1.5, cex.axis=1.25, main="Location of Spruce Moss Forest - Fen Transition Plots", cex.main=0.7)
#text( site.trunc$Latitude..WGS.84.datum., site.trunc$Longitude..WGS.84.datum)
These are my subjective habitat classifications for the 176 plots in the large sampling grid.
#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(size=4) +
theme_bw() + theme(panel.border = element_rect(colour="black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ xlab("Latitude") + ylab("Longitude") +
ggtitle("Habitat Types for Plots in Large Sampling Grid")
hab.ggplot+theme(axis.title.y = element_text(size=16, vjust=1),axis.title.x = element_text(size=16), axis.text.x = element_text(size=12), axis.text.y= element_text(size=12,angle=90,hjust=0.5))+
theme(plot.title = element_text(vjust = 2, size=15, face="bold")) + theme(legend.title = element_text(size=12))+
scale_linetype_discrete("Habitat Type") +
scale_shape_discrete("Habitat Type") +
scale_colour_discrete("Habitat Type")
The following was an unsuccessful attempt to create convex hulls for each habitat type.
hulls <- ddply(hab, "hab$Habitat.Description", find_hull)
hab.ggplot<-ggplot(hab, aes(hab$Latitude..WGS.84.datum., hab$Longitude..WGS.84.datum.,colour = hab$Habitat.Description)) +
geom_point(size=3) + geom_polygon(data = hulls,alpha = 0.4, colour= hab$Habitat.Description) +
theme_bw() + theme(panel.border = element_rect(colour="black"), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))+ xlab("Latitude") + ylab("Longitude") +
ggtitle("Convex Hulls of Habitat Types \nfor Plots in Large Sampling Grid")
hab.ggplot+theme(axis.title.y = element_text(size=16, vjust=1),axis.title.x = element_text(size=16), axis.text.x = element_text(size=12), axis.text.y= element_text(size=12,angle=90,hjust=0.5))+
theme(plot.title = element_text(vjust = 2, size=16, face="bold")) + theme(legend.title = element_text(size=12)) +
scale_linetype_discrete("Habitat Type") +
scale_shape_discrete("Habitat Type") +
scale_colour_discrete("Habitat Type")
The points in these two plots are coloured based on species richness.
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(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)
The following uses a colour ramp to show species richness across the large sampling grid. Plots with the highest species richness are indicated by black circles, whereas those with the lowest species richness are indicated by light gray.
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)
)
The following shows the mean, standard deviation, frequency distribution, and model diagnostics for species richness across the entire sampling grid. Species richness is not normally distributed as indicated by the histogram and Shapiro Wilks test.
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
The following plot shows the average species richness across each of the 16 sampling bands. Species richness is highest near the lake, which is nearest sampling band #1.
plot(rich.band, type="l",axes=FALSE, col="black", lwd=2, ylab="Average Species Richness", xlab="Sampling Band",las=1, cex.axis=1, 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)
box(bty="l", lwd=4)
The data underlying the species richness at each sampling band are not normally distributed.
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
Species richness generally decreases as you go up Mont Irony.
boxplot(richness ~ band, col="grey60")
This is the same boxplot, but reordered so that the bands with highest average richness values are on the right-hand side of the plot. Species richness is lowest towards the top of Mont Irony.
rich.med <- sort(tapply(richness, band, median))
boxplot(richness ~ factor(band, levels=names(rich.med)),col="grey60")
Plot showing effect sizes of species richness per band, ordered so that the highest species richness values are on the top.
# The best way to view the effect sizes graphically is to use plot.design()
plot.design(richness ~ band, ylab = "Species Richness", cex=0.75)
The following ANOVAs show that sampling bands differ significantly in average species richness.
# 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
Based on the following results, species richness per sampling band meets the assumptions of normality and homogeneous variance.
The Tukey test results show that species richness is significantly different between several of the sampling bands.
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 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