# load packages
library(pscl) # for McFadden R2
library(ggplot2) # for plotting effects
library(PerformanceAnalytics) # for correlations
library(car) # for Anova and VIF
library(multcomp) # for posthoc pairwise comparisons
library(effects) # for plotting effects from models
library(irr) # for calculating ICC
library(epiDisplay) # for odds ratios
library(ggeffects) # for marginal effects from models
library(ggpubr) # for composing figures

# import and view data; species presence denoted by 0=absent 1=present
# 88 arrays (initially 90 but 2 were excluded from analysis due to no detections)
data<-read.csv("Ethier_Wilson_data.csv")
head(data)

# Variable definitions
	# Array - a unique identifier given to each array deployment
	# columns 2-33 indicate presence (1) or absence (0) of each species at the array
	# MaxDBH - max diameter at breast height (cm) derived from 5 standard habitat sampling points
	# meanDensity - # trees/m^2 in 2-m radius, averaged among 5 standard habitat sampling points
	# meanCC - % canopy cover, averaged among 5 standard habitat sampling points
	# TotalVegRich - total number of species of trees, shrubs, and ericaceous plants encounterd
	# HabitatClass - arrays were categorized as spruce/fir (SPF), birch/poplar (BIR), or tamarack (LAR)
	# meanTemp - temperature (degrees Celsius), measured with Kestrel at 20-minute intervals and then averaged across 24 hours
	# JulianDate - day of the year
	# Year - 2016 or 2017
	# Rain - rainfall amount (mm) recorded at Goose Bay
	# Rain_y_n - whether any precipitation fell during the 24-hour recording (0=no; 1=yes)
	# DateCode - date when array was deployed (yyyy-mm-dd)
	# Wind_Kestrel - windspeed (km/h) measured with Kestrel near ground level inside array (or a nearby array) at 20-minute intervals, then averaged across 24 hours
	# Wind_EC - windspeed (km/h) measured by Environment Canada in open field at Goose Bay each hour then averaged over 24 hours

# make Year and HabitatClass factor variables
data$Year <- as.factor(data$Year)
data$HabitatClass <- as.factor(data$HabitatClass)

# calculate descriptives for Kestrel and ECCC wind speed (Methods: microphone arrays: paragraph 4)
length(data$Rain[data$Rain_y_n=="Y"]) # N=23 arrays with precipitation
summary(data$Rain[data$Rain_y_n=="Y"]) # amount of precipitation among those 23 arrays
sd(data$Rain[data$Rain_y_n=="Y"]) # SD precipitation among those 23 arrays
summary(data$Wind_Kestrel) # windspeed measured at ground level by Kestrel
sd(data$Wind_Kestrel) # sd windspeed measured at ground level by Kestrel
summary(data$Wind_EC) # windspeed measured in open at Goose Bay
sd(data$Wind_EC) # SD windspeed measured by ECCC in open at Goose Bay
cor.test(data$Wind_Kestrel, data$Wind_EC) # correlation between Kestrel and ECCC windspeed
plot(data$Wind_Kestrel, data$Wind_EC) # visualize weak positive correlation

# calculate bird species richness for each site
data$BirdRich<-rowSums(data[,2:33]!=0)

# calculate proportion of the 88 sites where each species was present; results paragraph 1 & figure 2
x <- colSums(data[,2:33]!=0)/nrow(data)
(x <- data.frame(proportion=sort(x, decreasing=T)))

# get count, mean, sd, min, max for each (a)biotic variable among the 88 sites; results paragraph 2
nrow(data) # number of arrays in analysis
sapply(data[,c("MaxDBH", "meanDensity", "meanCC", "TotalVegRich", "meanTemp", "JulianDate", "Wind_Kestrel", "Wind_EC", "BirdRich", "Rain")], mean, na.rm=T)
sapply(data[,c("MaxDBH", "meanDensity", "meanCC", "TotalVegRich", "meanTemp", "JulianDate", "Wind_Kestrel", "Wind_EC", "BirdRich", "Rain")], sd, na.rm=T)
sapply(data[,c("MaxDBH", "meanDensity", "meanCC", "TotalVegRich", "meanTemp", "JulianDate", "Wind_Kestrel", "Wind_EC", "BirdRich", "Rain")], min, na.rm=T)
sapply(data[,c("MaxDBH", "meanDensity", "meanCC", "TotalVegRich", "meanTemp", "JulianDate", "Wind_Kestrel", "Wind_EC", "BirdRich", "Rain")], max, na.rm=T)


#######################################################################################################################
# Figure 2
#######################################################################################################################
# adjust common names for use in figure 2
x$Species <- factor(c("Ruby-crowned Kinglet", "Dark-eyed Junco", "Yellow-rumped Warbler", "Pine Siskin", "Swainson's Thrush", "American Robin", "White-throated Sparrow", "Hermit Thrush", "Boreal Chickadee", "Fox Sparrow", "Orange-crowned Warbler", "Black-throated Green Warbler", "Tennessee Warbler", "Cape May Warbler", "Lincoln's Sparrow", "Pine Grosbeak", "Canada Jay", "Yellow-bellied Flycatcher", "Red-breasted Nuthatch", "Northern Waterthrush", "Alder Flycatcher", "Magnolia Warbler", "Brown Creeper", "Common Raven", "Bohemian Waxwing", "Canada Goose", "Common Loon", "American Crow", "American Redstart", "Blackpoll Warbler", "Common Nighthawk", "Winter Wren"))

# create figure
(fig2 <- ggplot(x, aes(x=reorder(Species, -proportion), y=proportion)) +
	geom_bar(stat = "identity", color="black", fill="white", width=0.4) +
	scale_y_continuous(expand = c(0, 0), limits = c(0, 0.85)) +
	theme(axis.title=element_text(size=10), axis.text=element_text(size=8, colour="black"), legend.position = "none", axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), axis.text.x=element_text(angle=90, vjust=0.5, hjust=1), axis.ticks.x=element_blank()) +
	labs(x="Species", y="Proportion of locations with detections"))

ggsave("Ethier & Wilson Figure 2.pdf", device="pdf", width=6, height=5)


#######################################################################################################################
##Species Richness Analysis##
#######################################################################################################################
# run linear model and assess terms using ANOVA and Summary functions; Results: Species richness: paragraph 1; Table 2
BirdRichglm<-lm(BirdRich~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+Wind_EC+JulianDate+HabitatClass, data = data)
car::Anova(BirdRichglm, type=3) # ANOVA table
summary(BirdRichglm) # coefficients table, R2, and overall model stats
plot(allEffects(BirdRichglm)) # plot all effects

# validate model using diagnostic plots and by calculating correlations and VIFs among predictors; looks good
chart.Correlation(as.matrix(data[, c("MaxDBH", "meanDensity", "meanCC", "TotalVegRich", "meanTemp", "Wind_EC", "JulianDate")]))
car::vif(BirdRichglm) # Methods: Statistical analysis: Species richness: paragraph 1
plot(data$meanTemp, fitted(BirdRichglm))
plot(data$Wind_EC, fitted(BirdRichglm))
plot(BirdRichglm)
hist(resid(BirdRichglm))


#######################################################################################################################
## Calculate ICC values for avian richness in 2016 and 2017 to determine the consistency in richness between years##
## Results: Species richness: paragraph 2
#######################################################################################################################
data.rep <- read.csv("Ethier_Wilson_repdata.csv") # read repeated measures data
data.rep <- data.rep[order(data.rep$Pair),] # order data.rep by 'Pair'
head(data.rep)

# Variable definitions
	# Site - a unique identifier given to each array deployment
	# Pair - a unique identifier given to each array location (there were two arrays deployed per location for 20 locations)
	# Year - 2016 or 2017
	# Julian - day of the year
	# columns 5-26 indicate presence (1) or absence (0) of each species at the array

# calculate species richness at each of the 20 sites in each year; species presence indicated as absent=0 present=1
data.rep$BirdRich <- rowSums(data.rep[,5:26]!=0)
data2016 <- subset(data.rep, Year==2016) # create 2016 subset
data2017 <- subset(data.rep, Year==2017) # create 2017 subset

# get descriptives for species richness per year
mean(data2016$BirdRich)
sd(data2016$BirdRich)
min(data2016$BirdRich)
max(data2016$BirdRich)
mean(data2017$BirdRich)
sd(data2017$BirdRich)
min(data2017$BirdRich)
max(data2017$BirdRich)

# Calculate ICC values for avian richness in 2016 and 2017 to determine the consistency in richness between years
IRRdata <- as.matrix(cbind(data2016$BirdRich, data2017$BirdRich))
icc(IRRdata, model = "twoway", type = "consistency") # using 2-way ANOVA, which assumes both subjects and years are
# chosen from larger pools; use 'consistency' if only interested in correlations, 'agreement' if also interested in
# differences among means among years; 'consistency' and 'two-way' are most appropriate here, though 'agreement'
# yields similar results.


#######################################################################################################################
# Figure 3
#######################################################################################################################
# Figure 3a - species richness vs temperature vs date
(fig3a <- ggplot()
	+ geom_point(data=data, aes(x=JulianDate, y=meanTemp, size=BirdRich, fill=Year, colour=Year), shape=21, alpha=0.4, show.legend=T)
	+ guides(fill="none")
	+ guides(colour="none")
#	+ scale_fill_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_fill_manual(values=c("2016"="black", "2017"="white")) # black and white
#	+ scale_colour_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_colour_manual(values=c("2016"="black", "2017"="black")) # black and white
	+ theme(axis.title=element_text(size=11), axis.text=element_text(size=9, colour="black"), axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank(), legend.position=c(0.65,0.08), legend.direction="horizontal", legend.spacing.x=unit(0, 'cm'), legend.title=element_text(colour="black", size=9), legend.text=element_text(colour="black", size=9), legend.key=element_rect(fill="white"))
	+ labs(x="Day of year", y=expression(Mean~daily~temperature~"("*degree*C*")"), size="# species"))

# Figure 3b - species richness vs mean daily temperature
fig3b.EMM <- ggeffect(BirdRichglm, terms = "meanTemp [all]", ci.lvl = 0.95)
(fig3b <- ggplot(fig3b.EMM, aes(x=x, y=predicted))
	+ geom_line()
	+ geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill="gray", alpha=0.6)
	+ geom_point(data=data, aes(x=meanTemp, y=BirdRich, fill=Year, colour=Year), shape=21, size=2.5, alpha=0.4)
#	+ scale_fill_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_fill_manual(values=c("2016"="black", "2017"="white")) # black and white
#	+ scale_colour_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_colour_manual(values=c("2016"="black", "2017"="black")) # black and white
	+ theme(axis.title=element_text(size=11), axis.text=element_text(size=9, colour="black"), legend.position = "none", axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
	+ labs(x=expression(Mean~daily~temperature~"("*degree*C*")"), y="Number of avian species"))

# Figure 3c - species richness vs day of year
fig3c.EMM <- ggeffect(BirdRichglm, terms = "JulianDate [all]", ci.lvl = 0.95)
(fig3c <- ggplot(fig3c.EMM, aes(x=x, y=predicted))
	+ geom_line()
	+ geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill="gray", alpha=0.6)
	+ geom_point(data=data, aes(x=JulianDate, y=BirdRich, fill=Year, colour=Year), shape=21, size=2.5, alpha=0.4)
#	+ scale_fill_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_fill_manual(values=c("2016"="black", "2017"="white")) # black and white
#	+ scale_colour_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_colour_manual(values=c("2016"="black", "2017"="black")) # black and white
	+ theme(axis.title=element_text(size=11), axis.text=element_text(size=9, colour="black"), legend.position = "none", axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
	+ labs(x="Day of year", y="Number of avian species"))

# Figure 3d - species richness vs mean daily wind speed at ECCC station
fig3d.EMM <- ggeffect(BirdRichglm, terms = "Wind_EC [all]", ci.lvl = 0.95)
(fig3d <- ggplot(fig3d.EMM, aes(x=x, y=predicted))
	+ geom_line()
	+ geom_ribbon(aes(ymin = conf.low, ymax = conf.high), fill="gray", alpha=0.6)
	+ geom_point(data=data, aes(x=Wind_EC, y=BirdRich, fill=Year, colour=Year), shape=21, size=2.5, alpha=0.4)
#	+ scale_fill_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_fill_manual(values=c("2016"="black", "2017"="white")) # black and white
#	+ scale_colour_manual(values=c("2016"="purple", "2017"="green")) # colour
	+ scale_colour_manual(values=c("2016"="black", "2017"="black")) # black and white
	+ theme(axis.title=element_text(size=11), axis.text=element_text(size=9, colour="black"), legend.position = "none", axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank())
	+ labs(x="Mean daily wind speed (km/h)", y="Number of avian species"))

# arrange the 3 figures in a multipanel plot with 1 rows and 3 columns
ggarrange(fig3a, fig3b, fig3c, fig3d, nrow=2, ncol=2, align="hv", labels="auto")
# ggsave("Ethier & Wilson Figure 3.pdf", device="pdf", width=7.08, height=7.08) # use for colour version
ggsave("Ethier & Wilson Figure 3 bw.pdf", device="pdf", width=7.08, height=7.08) # use for black-and-white version


#######################################################################################################################
## Species-specific detection models##
## Results: Species-specific detection: paragraph 1; Table 3; Table S2
#######################################################################################################################

# get numbers of arrays per habitat type per species for Table S2
table(data$HabitatClass) # across species
table(data$HabitatClass[data$Yellow.bellied.Flycatcher==1])
table(data$HabitatClass[data$Whiskey.Jack==1])
table(data$HabitatClass[data$Boreal.Chickadee==1])
table(data$HabitatClass[data$Ruby.crowned.Kinglet==1])
table(data$HabitatClass[data$Swainson.s.Thrush==1])
table(data$HabitatClass[data$Hermit.Thrush==1])
table(data$HabitatClass[data$American.Robin==1])
table(data$HabitatClass[data$Pine.Grosbeak==1])
table(data$HabitatClass[data$Pine.Siskin==1])
table(data$HabitatClass[data$Fox.Sparrow==1])
table(data$HabitatClass[data$Dark.eyed.Junco==1])
table(data$HabitatClass[data$Lincoln.s.Sparrow==1])
table(data$HabitatClass[data$White.throated.Sparrow==1])
table(data$HabitatClass[data$Tennessee.Warbler==1])
table(data$HabitatClass[data$Orange.crowned.Warbler==1])
table(data$HabitatClass[data$Cape.May.Warbler==1])
table(data$HabitatClass[data$Yellow.rumped.Warbler==1])
table(data$HabitatClass[data$Black.throated.Green.Warbler==1])

# limit analysis to species that were present in at least 10% of sites (9 or more of the 88 sites)
# this excludes 14 of 32 species, including Alder Flycatcher, American Crow, American Redstart, Blackpoll Warbler, Bohemian Waxwing, Brown Creeper, Canada Goose, Common Loon, Common Nighthawk, Common Raven, Magnolia Warbler, Northern Waterthrush, Red-breasted Nuthatch, and Winter Wren
# this leaves 18 of 32 species remaining
colSums(data[,2:33]!=0) # show number of trials in which each species was detected

## the following analyses follow a consistent pattern of running a GLM, getting a summary, determining pseudo-R2 and testing model fit##
# for models with significant predictors, calculate descriptives (odds ratios and their 95% CIs) for the predictor

##Yellow-bellied Flycatcher##
YBFLglm1<-glm(Yellow.bellied.Flycatcher~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(YBFLglm1)
summary(multcomp::glht(YBFLglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(YBFLglm1, crude=F, simplified=T)
pR2(YBFLglm1)
plot(allEffects(YBFLglm1)) # plot all effects

##Whiskey Jack/Grey Jay##
WHJAglm1<-glm(Whiskey.Jack~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(WHJAglm1)
summary(multcomp::glht(WHJAglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(WHJAglm1, crude=F, simplified=T)
pR2(WHJAglm1)
plot(allEffects(WHJAglm1)) # plot all effects

##Boreal Chickadee##
BOCHglm1<-glm(Boreal.Chickadee~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(BOCHglm1)
summary(multcomp::glht(BOCHglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(BOCHglm1, crude=F, simplified=T)
pR2(BOCHglm1)
plot(allEffects(BOCHglm1)) # plot all effects

##Ruby-crowned Kinglet##
RCKIglm1<-glm(Ruby.crowned.Kinglet~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(RCKIglm1)
summary(multcomp::glht(RCKIglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(RCKIglm1, crude=F, simplified=T)
pR2(RCKIglm1)
plot(allEffects(RCKIglm1)) # plot all effects

##Swainson's Thrush##
SWTHglm1<-glm(Swainson.s.Thrush~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(SWTHglm1)
summary(multcomp::glht(SWTHglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(SWTHglm1, crude=F, simplified=T)
pR2(SWTHglm1)
plot(allEffects(SWTHglm1)) # plot all effects

##Hermit Thrush##
HETHglm1<-glm(Hermit.Thrush~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(HETHglm1)
summary(multcomp::glht(HETHglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(HETHglm1, crude=F, simplified=T)
pR2(HETHglm1)
plot(allEffects(HETHglm1)) # plot all effects

##American Robin##
AMROglm1<-glm(American.Robin~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(AMROglm1)
summary(multcomp::glht(AMROglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(AMROglm1, crude=F, simplified=T)
pR2(AMROglm1)
table(data$American.Robin, data$HabitatClass) # distribution among forest types
plot(allEffects(AMROglm1)) # plot all effects

##Pine Grosbeak##
PIGRglm1<-glm(Pine.Grosbeak~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(PIGRglm1)
summary(multcomp::glht(PIGRglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(PIGRglm1, crude=F, simplified=T)
pR2(PIGRglm1)
plot(allEffects(PIGRglm1)) # plot all effects

##Pine Siskin##
PISIglm1<-glm(Pine.Siskin~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(PISIglm1)
summary(multcomp::glht(PISIglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(PISIglm1, crude=F, simplified=T)
pR2(PISIglm1)
plot(allEffects(PISIglm1)) # plot all effects

##Fox Sparrow##
FOSPglm1<-glm(Fox.Sparrow~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(FOSPglm1)
summary(multcomp::glht(FOSPglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(FOSPglm1, crude=F, simplified=T)
pR2(FOSPglm1)
plot(allEffects(FOSPglm1)) # plot all effects

##Dark-eyed Junco##
DEJUglm1<-glm(Dark.eyed.Junco~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(DEJUglm1)
summary(multcomp::glht(DEJUglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(DEJUglm1, crude=F, simplified=T)
pR2(DEJUglm1)
plot(allEffects(DEJUglm1)) # plot all effects

##Lincoln's Sparrow##
LISPglm1<-glm(Lincoln.s.Sparrow~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(LISPglm1)
summary(multcomp::glht(LISPglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(LISPglm1, crude=F, simplified=T)
pR2(LISPglm1)
plot(allEffects(LISPglm1)) # plot all effects

##White-throated Sparrow##
WTSPglm1<-glm(White.throated.Sparrow~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(WTSPglm1)
summary(multcomp::glht(WTSPglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(WTSPglm1, crude=F, simplified=T)
pR2(WTSPglm1)
plot(allEffects(WTSPglm1)) # plot all effects

##Tennessee Warbler##
TEWAglm1<-glm(Tennessee.Warbler~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(TEWAglm1)
summary(multcomp::glht(TEWAglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(TEWAglm1, crude=F, simplified=T)
pR2(TEWAglm1)
plot(allEffects(TEWAglm1)) # plot all effects

##Orange-crowned Warbler##
OCWAglm1<-glm(Orange.crowned.Warbler~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(OCWAglm1)
summary(multcomp::glht(OCWAglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(OCWAglm1, crude=F, simplified=T)
pR2(OCWAglm1)
plot(allEffects(OCWAglm1)) # plot all effects

##Cape May Warbler##
CMWAglm1<-glm(Cape.May.Warbler~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(CMWAglm1)
summary(multcomp::glht(CMWAglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(CMWAglm1, crude=F, simplified=T)
pR2(CMWAglm1)
plot(allEffects(CMWAglm1)) # plot all effects

##Yellow-rumped Warbler##
YRWAglm1<-glm(Yellow.rumped.Warbler~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(YRWAglm1)
summary(multcomp::glht(YRWAglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(YRWAglm1, crude=F, simplified=T)
pR2(YRWAglm1)
plot(allEffects(YRWAglm1)) # plot all effects

##Black-throated Green Warbler##
BTGWglm1<-glm(Black.throated.Green.Warbler~MaxDBH+meanDensity+meanCC+TotalVegRich+meanTemp+JulianDate+Wind_EC+HabitatClass, data=data,family=binomial(link='logit'))
summary(BTGWglm1)
summary(multcomp::glht(BTGWglm1, linfct = mcp(HabitatClass = "Tukey"))) # pairwise comparisons for forest type
logistic.display(BTGWglm1, crude=F, simplified=T)
pR2(BTGWglm1)
table(data$Black.throated.Green.Warbler, data$HabitatClass) # show presence by forest type
plot(allEffects(BTGWglm1)) # plot all effects


#######################################################################################################################
## Spatial consistency in species composition##
## Results: Spatial and temporal consistency in avian community composition: Paragraph1
#######################################################################################################################

# create duplicate of data, but excluding species present in less than 10% of sites; also exclude BirdRich and re-calculate it based on remaining species
# this leaves 18 of 32 species remaining
data2 <- subset(data, select=-c(Alder.Flycatcher, American.Crow, American.Redstart, Blackpoll.Warbler, Bohemian.Waxwing, Brown.Creeper, Canada.Goose, Common.Loon, Common.Nighthawk, Common.Raven, Magnolia.Warbler, Northern.Waterthrush, Red.breasted.Nuthatch, Winter.Wren, BirdRich))
data2$BirdRich2 <-rowSums(data2[,2:19]!=0)
head(data2) # view data

data2 <- data2[order(data2$Array),] # order data2 by array number
data.birch <- data2[data2$HabitatClass == "BIR", ] # create subset of data2 for birch sites
data.tamarack <- data2[data2$HabitatClass == "LAR", ] # create subset of data2 for tamarack sites
data.spruce <- data2[data2$HabitatClass == "SPF", ] # create subset of data2 for spruce sites

# create an empty matrix for each pair of habitat varieties (e.g., bvt = birch versus tamarack),
# including control in which sites of a given habitat are compared to each other
diversity.matrix.bvt <- matrix(nrow = nrow(data.birch), ncol = nrow(data.tamarack))
diversity.matrix.svt <- matrix(nrow = nrow(data.spruce), ncol = nrow(data.tamarack))
diversity.matrix.bvs <- matrix(nrow = nrow(data.birch), ncol = nrow(data.spruce))
diversity.matrix.bvb <- matrix(nrow = nrow(data.birch), ncol = nrow(data.birch))
diversity.matrix.tvt <- matrix(nrow = nrow(data.tamarack), ncol = nrow(data.tamarack))
diversity.matrix.svs <- matrix(nrow = nrow(data.spruce), ncol = nrow(data.spruce))

# populate all cells in bvt matrix
for (i in 1:nrow(data.birch)) {
	for (j in 1:nrow(data.tamarack)) {
		
		birch.site.vector <- as.vector(t(data.birch[i, 2:19]))
		tamarack.site.vector <- as.vector(t(data.tamarack[j, 2:19]))
		site.diffs <- birch.site.vector + tamarack.site.vector
		
		a <- sum(site.diffs == 2)
		bc <- sum(site.diffs == 1)
		
		diversity.matrix.bvt[i,j] <- 2*a / (2*a + bc)
				
	}
}

# populate all cells in svt matrix
for (i in 1:nrow(data.spruce)) {
	for (j in 1:nrow(data.tamarack)) {
		
		spruce.site.vector <- as.vector(t(data.spruce[i, 2:19]))
		tamarack.site.vector <- as.vector(t(data.tamarack[j, 2:19]))
		site.diffs <- spruce.site.vector + tamarack.site.vector
		
		a <- sum(site.diffs == 2)
		bc <- sum(site.diffs == 1)
		
		diversity.matrix.svt[i,j] <- 2*a / (2*a + bc)
				
	}
}

# populate all cells in bvs matrix
for (i in 1:nrow(data.birch)) {
	for (j in 1:nrow(data.spruce)) {
		
		birch.site.vector <- as.vector(t(data.birch[i, 2:19]))
		spruce.site.vector <- as.vector(t(data.spruce[j, 2:19]))
		site.diffs <- birch.site.vector + spruce.site.vector
		
		a <- sum(site.diffs == 2)
		bc <- sum(site.diffs == 1)
		
		diversity.matrix.bvs[i,j] <- 2*a / (2*a + bc)
				
	}
}

# populate all cells in bvb matrix
for (i in 1:nrow(data.birch)) {
  for (j in 1:nrow(data.birch)) {
    
    birch.site.vector1 <- as.vector(t(data.birch[i, 2:19]))
    birch.site.vector2 <- as.vector(t(data.birch[j, 2:19]))
    site.diffs <- birch.site.vector1 +  birch.site.vector2
    
    a <- sum(site.diffs == 2)
    bc <- sum(site.diffs == 1)
    
    diversity.matrix.bvb[i,j] <- 2*a / (2*a + bc)
    
  }
}

# populate all cells in tvt matrix
for (i in 1:nrow(data.tamarack)) {
  for (j in 1:nrow(data.tamarack)) {
    
    tamarack.site.vector1 <- as.vector(t(data.tamarack[i, 2:23]))
    tamarack.site.vector2 <- as.vector(t(data.tamarack[j, 2:23]))
    site.diffs <- tamarack.site.vector1 +  tamarack.site.vector2
    
    a <- sum(site.diffs == 2)
    bc <- sum(site.diffs == 1)
    
    diversity.matrix.tvt[i,j] <- 2*a / (2*a + bc)
    
  }
}

# populate all cells in svs matrix
for (i in 1:nrow(data.spruce)) {
  for (j in 1:nrow(data.spruce)) {
    
    spruce.site.vector1 <- as.vector(t(data.spruce[i, 2:23]))
    spruce.site.vector2 <- as.vector(t(data.spruce[j, 2:23]))
    site.diffs <- spruce.site.vector1 +  spruce.site.vector2
    
    a <- sum(site.diffs == 2)
    bc <- sum(site.diffs == 1)
    
    diversity.matrix.svs[i,j] <- 2*a / (2*a + bc)
    
  }
}

# remove diagonal containing 1s from the 3 control matrices
diag(diversity.matrix.bvb) <- NA
diag(diversity.matrix.tvt) <- NA
diag(diversity.matrix.svs) <- NA

# calculate mean of all cells in a given matrix (Cs values)
mean(diversity.matrix.bvt)
mean(diversity.matrix.svt)
mean(diversity.matrix.bvs)
mean(diversity.matrix.bvb, na.rm=T) # need to ignore diagonal containing NAs
mean(diversity.matrix.tvt, na.rm=T) # need to ignore diagonal containing NAs
mean(diversity.matrix.svs, na.rm=T) # need to ignore diagonal containing NAs

# calculate sd of all cells in a given matrix
sd(diversity.matrix.bvt)
sd(diversity.matrix.svt)
sd(diversity.matrix.bvs)
sd(diversity.matrix.bvb, na.rm=T) # need to ignore diagonal containing NAs
sd(diversity.matrix.tvt, na.rm=T) # need to ignore diagonal containing NAs
sd(diversity.matrix.svs, na.rm=T) # need to ignore diagonal containing NAs


#######################################################################################################################
##Temporal consistency in species composition##
## Results: Spatial and temporal consistency in avian community composition: Paragraph1
#######################################################################################################################

# adjust data2016 and data2017 by excluding species that were present in less than 10% of the original 88 sites (i.e., < 9)
# this leaves 18 of 32 species remaining; also exclude BirdRich column
data2016_10 <- subset(data2016, select=-c(Alder.Flycatcher, Magnolia.Warbler, Northern.Waterthrush, Red.breasted.Nuthatch, BirdRich))
data2017_10 <- subset(data2017, select=-c(Alder.Flycatcher, Magnolia.Warbler, Northern.Waterthrush, Red.breasted.Nuthatch, BirdRich))

# create a dataframe to store the 20 Ct values for the 20 sites sampled in 2016 and 2017
diversity.matrix.16v17 <- cbind(data2016_10[,1:4], data2017_10[,1:4])
diversity.matrix.16v17$Ct <- NA
colnames(diversity.matrix.16v17)<-c("Site.2016", "Pair.2016", "Year.2016", "Julian.2016", "Site.2017", "Pair.2017", "Year.2017", "Julian.2017", "Ct")
diversity.matrix.16v17$Julian.mean <-rowMeans(diversity.matrix.16v17[c("Julian.2016", "Julian.2017")])

# populate all cells in 16v17 matrix
for (i in 1:nrow(data2016_10)) {
		
	vector2016 <- as.vector(t(data2016_10[i, 5:22]))
	vector2017 <- as.vector(t(data2017_10[i, 5:22]))
	year.diffs <- vector2016 + vector2017
		
	a <- sum(year.diffs == 2)
	bc <- sum(year.diffs == 1)
		
	diversity.matrix.16v17[i,"Ct"] <- 2*a / (2*a + bc)

}

# calculate mean and sd of all cells in diversity.matrix.16v17 (Cs values)
mean(diversity.matrix.16v17$Ct)
sd(diversity.matrix.16v17$Ct)
diversity.matrix.16v17 # show all values
range(diversity.matrix.16v17$Ct) # get range of values

# test for relationship between temporal turnover and when in the season the location was sampled
plot(Ct~Julian.mean, data=diversity.matrix.16v17)
summary(lm(Ct~Julian.mean, data=diversity.matrix.16v17))