---
title: "Drought"
author: "Emma Smith"
date: "`r Sys.Date()`"
output:
  html_document: default
  word_document: default
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r custom theme}
#create custom theme that will be used in all figures
theme_custom <- function(base_size = 12, base_family = "")
{
  theme_bw(base_size = base_size, base_family = base_family) %+replace% 
    theme(line = element_line(colour = "black", size = 0.5, 
                              linetype = 1, lineend = "butt"), 
          text = element_text(family = base_family, face = "plain",
                              colour = "black", size = base_size,
                              lineheight = 0.9,  hjust = 0.5,
                              vjust = 0.5, angle = 0, 
                              margin = margin(), debug = FALSE), 
          axis.text = element_text(size = rel(0.8),color = "black"),
          axis.ticks = element_line(colour = "black"), 
          axis.line = element_line(color = "black"),
          legend.background = element_rect(fill = NA, color = NA),
          legend.key = element_rect(color = NA), 
          plot.background = element_rect(fill = "white"),
          plot.title = element_text(face = "bold"),
          panel.background = element_rect(fill = "white", colour = "white"), 
          panel.border = element_rect(fill = NA, colour = "white"), 
          panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(), 
          strip.background = element_rect(fill = NA, colour = NA, size = 0.2),
          strip.text = element_text(color = "black"))
}

#color scale
colorvars   <- c(2014,2015,2016,2017,2018,2019)
colorvarscols <- c("gray21","gray21","gray21","gray21","red","gray21")
colors <- data.frame(colorvars,colorvarscols)
colors <- as.vector(colors$colorvarscols)
colors
```

```{r packages}
#load required packages
library(ggplot2) #used to create all figures
library(readxl) #to import data from Excel files
library(reshape2) #used to reformat data, useful for some analyses
library(tidyr) #used to reformat data, useful for some analyses
library(tidyverse) #large group of packages, useful for many things
library(gridExtra) #used to combine separate figures into a single image
library(circular) #used to run circular statistics
library(lubridate) #used to reformat dates
library(RMark) #used to run nest survival models
library(effectsize) #used to calculate effect sizes
```

```{r climate}
#upload data directly from the LTER data portal
climate <- read.csv("http://lter.konza.ksu.edu/sites/default/files/data/AWE012.csv")

#subset data to only include 2014-2019, and only May-July
climate <- climate[(climate$RECYEAR == "2014" | climate$RECYEAR == "2015" | climate$RECYEAR == "2016" | climate$RECYEAR == "2017" | climate$RECYEAR == "2018" | climate$RECYEAR == "2019"),]
climate <- climate[(climate$RECMONTH == "5" | climate$RECMONTH == "6" | climate$RECMONTH == "7"),]

#remove unneeded columns to reduce clutter
climate <- climate[, !(colnames(climate) %in% c("SMAX","SMIN","S_AVE","WAVE","DATACODE","RECTYPE","WATERSHED", "RECMONTH", "RECDAY", "DAYOFYEAR","TMIN","TAVE","DHUMID","DSRAD"))]

#reclassify columns that are incorrect class
climate$DPPT <- as.numeric(as.character(climate$DPPT))
climate$RECYEAR <- as.factor(climate$RECYEAR)
climate$TMAX <- as.numeric(as.character(climate$TMAX))

#remove NA values
climate <- na.omit(climate)

#add in a group column to separate drought and non-drought years
climate$group[climate$RECYEAR != "2018"] <- "nondrought"
climate$group[climate$RECYEAR == "2018"] <- "drought"

#calculate means, standard deviations, and standard error
climatesummary <- climate %>% group_by(RECYEAR) %>%
    summarise_each(funs(mean,sd,se=sd(.)/sqrt(n())))

#add in a group column to separate drought and non-drought years
climatesummary$group[climatesummary$RECYEAR != "2018"] <- "nondrought"
climatesummary$group[climatesummary$RECYEAR == "2018"] <- "drought"

#calculate means by group, to differentiate drought and nondrought years
climatesummary <- climatesummary %>% group_by(group) %>% summarize_each(funs(mean))

```

```{r precip}
#extract only the year and daily precip amount columns
precip <- climate %>% select(1,3)

#calculate sum total precip for each year
precipsummary <- precip %>% group_by(RECYEAR) %>% summarize(SumPrecip = sum(DPPT))

#add in a group column to separate drought and non-drought years
precipsummary$group[precipsummary$RECYEAR != "2018"] <- "nondrought"
precipsummary$group[precipsummary$RECYEAR == "2018"] <- "drought"

#calculate mean for non-drought years
precipsummary <- precipsummary %>% group_by(group) %>% summarize(precip = mean(SumPrecip), sd = sd(SumPrecip))
```

```{r temperature}
#ANOVA of mean daily maximum temperature ~ Year
anova <- aov(TMAX ~ RECYEAR, data = climate)
summary(anova)

effectsize(anova)
```

```{r nests}
#import data
nests <- read_xlsx("D://Users/inkli/Documents/Data/Drought paper/Nests.xlsx", col_names = TRUE)

#add in a group column to separate drought and non-drought years
nests$group <- " nondrought"
nests$group[nests$Year == "2018"] <- "drought"

#add a date (mm-dd) column by extracting from date 1st egg column 
nests$Date <- nests$Date1stEgg
nests$Date <- substr(nests$Date,6,10)
nests$Date <- as.Date(nests$Date,"%m-%d")

#rename species columns to spell out entire name, and add a space to sparrows and meadowlarks to force them to be in order
nests$Species[nests$Species  ==  "GRSP"] <- " Grasshopper\nSparrow"
nests$Species[nests$Species  ==  "DICK"] <- "Dickcissel"
nests$Species[nests$Species  ==  "EAME"] <- "  Eastern\nMeadowlark"

#number of young fledged per nest
numfledgelings <- as.data.frame(cbind(nests$Species, as.numeric(nests$Year), nests$NHostFledge))
names(numfledgelings) <- c("Species", "Year", "NHostFledge")
numfledgelings$NHostFledge <- as.numeric(numfledgelings$NHostFledge)
numfledgelings <- numfledgelings %>% na.omit()
numfledgelings <- numfledgelings %>% group_by(Year, Species) %>% summarize(nfledge = mean(NHostFledge), sd = sd(NHostFledge), se = (sd(NHostFledge)/sqrt(length(NHostFledge))))

#add in a group column to separate drought and non-drought years
numfledgelings$group[numfledgelings$Year != "2018"] <- "nondrought"
numfledgelings$group[numfledgelings$Year == "2018"] <- "drought"


nests$Year <- as.factor(nests$Year)
nests$JulianDate <- as.numeric(yday(nests$Date1stEgg))
EAMEdata <- as.data.frame(subset(nests, Species == "  Eastern Meadowlark"))
GRSPdata <- as.data.frame(subset(nests, Species == " Grasshopper Sparrow"))
DICKdata <- as.data.frame(subset(nests, Species == "Dickcissel"))

```

```{r nest initiation}
nestinitiation <- nests

nestinitiation$Year <- factor(nestinitiation$Year, levels = c("2014","2015","2016","2017","2019","2018"))

colorvarsini   <- c("2014","2015","2016","2017","2018","2019")
colorvarscolsini <- c("gray21","gray21","gray21","gray21","red", "gray21")
colorsini <- data.frame(colorvarsini,colorvarscolsini)
colorsini <- as.vector(colorsini$colorvarscols)
colorsini
#nest initiation smoothed curve of counts
ggplot() +
  geom_density(data = subset(nestinitiation, as.Date(Date1stEgg) >= "2014-1-1"), aes(Date, fill = Year), alpha = .5, adjust = 1/3, trim = TRUE) +
  labs(x = "Julian Date" , y = "Number of Nests\nInitiated Per Day") + 
  theme_custom() +
  scale_fill_manual(values = c(colorsini)) +
  facet_grid(rows = vars(Species)) +
  theme(legend.position = "none")
```

```{r nest success, message=FALSE}

# Nest success models, adapted from killdeer2RMark.R

# input file
nest_model <- read.csv("~/Data/Drought paper/nest_model_data.csv")
nest_model <- nest_model %>% select(1,3:9)
names(nest_model)[1] <- "NestID"

#GRSP----
nest_model_GRSP <- filter(nest_model, Species == "GRSP")
# look at top 25 lines of file to confirm input is okay
head(nest_model_GRSP, 25)

# look at first and last occasions
# needs to be 1 for dataset
(min = min(nest_model_GRSP$FirstFound))
# this gives maximum count for number of occasions
(occ = max(nest_model_GRSP$LastChecked))

# check input for errors
nest_model_GRSP$Diff = nest_model_GRSP$LastChecked - nest_model_GRSP$LastPresent
nest_model_GRSP$Error = "No"
nest_model_GRSP$Error[(nest_model_GRSP$LastPresent - nest_model_GRSP$FirstFound) < 0] <- "Yes"
nest_model_GRSP$Error[(nest_model_GRSP$LastChecked - nest_model_GRSP$LastPresent) < 0] <- "Yes"
nest_model_GRSP$Error[nest_model_GRSP$Fate == 1 & nest_model_GRSP$Diff <= 0] <- "Yes"
nest_model_GRSP$Error[nest_model_GRSP$Fate == 0 & nest_model_GRSP$Diff != 0] <- "Yes"
(Error <- subset(nest_model_GRSP, Error == "Yes"))
# drop errors if present
nest_model_GRSP <- nest_model_GRSP[which(nest_model_GRSP$Error == 'No'), ]

## Check for multicollinearity by comparing r values
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  # correlation coefficient
  r <- cor(x, y)
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste("r= ", txt, sep = "")
  text(0.5, 0.6, txt)
  
  # p-value calculation
  p <- cor.test(x, y)$p.value
  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
  txt2 <- paste("p= ", txt2, sep = "")
  if ( p < 0.01) txt2 <- paste("p= ", "<0.01", sep = "")
  text(0.5, 0.4, txt2)
}

nest_model_GRSP$Year = factor(nest_model_GRSP$Year)
nest_model_GRSP$Species = factor(nest_model_GRSP$Species)
nest_model_GRSP$FirstFound = as.numeric(nest_model_GRSP$FirstFound)

# Write a function to build the model
run.nest_model_GRSP.year = function()
{
  # Year
  Year = mark(nest_model_GRSP, nocc = occ, model = "Nest", model.parameters = list(S = list(formula = ~Year)), groups = c("Year"))
  return(collect.models())
  
}

nest_model_GRSP.results.year = run.nest_model_GRSP.year() # This runs the models above and takes a minute or 2
nest_model_GRSP.results.year = adjust.chat(1,nest_model_GRSP.results.year)

GRSP_DSR <- as.data.frame(nest_model_GRSP.results.year$Year$results$real)

GRSP_DSR$year <- c("2014", "2015", "2016", "2017", "2018", "2019")

#DICK----
nest_model_DICK <- subset(nest_model, Species == "DICK")
# look at top 25 lines of file to confirm input is okay
head(nest_model_DICK, 25)

#adjust to correct dates
nest_model_DICK$FirstFound <- nest_model_DICK$FirstFound - 14
nest_model_DICK$LastPresent <- nest_model_DICK$LastPresent - 14
nest_model_DICK$LastChecked <- nest_model_DICK$LastChecked - 14
# look at first and last occasions
# needs to be 1 for dataset
(min = min(nest_model_DICK$FirstFound))  
# this gives maximum count for number of occasions
(occ = max(nest_model_DICK$LastChecked))

# check input for errors
nest_model_DICK$Diff = nest_model_DICK$LastChecked - nest_model_DICK$LastPresent
nest_model_DICK$Error = "No"
nest_model_DICK$Error[(nest_model_DICK$LastPresent - nest_model_DICK$FirstFound) < 0] <- "Yes"
nest_model_DICK$Error[(nest_model_DICK$LastChecked - nest_model_DICK$LastPresent) < 0] <- "Yes"
nest_model_DICK$Error[nest_model_DICK$Fate == 1 & nest_model_DICK$Diff <= 0] <- "Yes"
nest_model_DICK$Error[nest_model_DICK$Fate == 0 & nest_model_DICK$Diff != 0] <- "Yes"
(Error <- subset(nest_model_DICK, Error == "Yes"))
# drop errors if present
nest_model_DICK <- nest_model_DICK[which(nest_model_DICK$Error == 'No'), ]

## Check for multicollinearity by comparing r values
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  # correlation coefficient
  r <- cor(x, y)
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste("r= ", txt, sep = "")
  text(0.5, 0.6, txt)
  
  # p-value calculation
  p <- cor.test(x, y)$p.value
  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
  txt2 <- paste("p= ", txt2, sep = "")
  if ( p < 0.01) txt2 <- paste("p= ", "<0.01", sep = "")
  text(0.5, 0.4, txt2)
}

nest_model_DICK$Year = factor(nest_model_DICK$Year)
nest_model_DICK$Species = factor(nest_model_DICK$Species)
nest_model_DICK$FirstFound = as.numeric(nest_model_DICK$FirstFound)


# Write a function for evaluating a set of competing models for the year hypothesis
run.nest_model_DICK.year = function()
{
  # Year
  Year = mark(nest_model_DICK, nocc = occ, model = "Nest", model.parameters = list(S = list(formula = ~Year)), groups = c("Year"))
  return(collect.models())
  
}

nest_model_DICK.results.year = run.nest_model_DICK.year() # This runs the models above and takes a minute or 2
nest_model_DICK.results.year = adjust.chat(1, nest_model_DICK.results.year)
DICK_DSR <- as.data.frame(nest_model_DICK.results.year$Year$results$real)
DICK_DSR$year <- c("2017", "2018", "2019")
#EAME----
nest_model_EAME <- subset(nest_model, Species == "EAME")
# look at top 25 lines of file to confirm input is okay
head(nest_model_EAME, 25)

# look at first and last occasions
# needs to be 1 for dataset
(min = min(nest_model_EAME$FirstFound))  
# this gives maximum count for number of occasions
(occ = max(nest_model_EAME$LastChecked))

# check input for errors
nest_model_EAME$Diff = nest_model_EAME$LastChecked - nest_model_EAME$LastPresent
nest_model_EAME$Error = "No"
nest_model_EAME$Error[(nest_model_EAME$LastPresent - nest_model_EAME$FirstFound) < 0 ] <- "Yes"
nest_model_EAME$Error[(nest_model_EAME$LastChecked - nest_model_EAME$LastPresent) < 0 ] <- "Yes"
nest_model_EAME$Error[nest_model_EAME$Fate == 1 & nest_model_EAME$Diff <= 0] <- "Yes"
nest_model_EAME$Error[nest_model_EAME$Fate == 0 & nest_model_EAME$Diff != 0] <- "Yes"
(Error <- subset(nest_model_EAME, Error == "Yes"))
# drop errors if present
nest_model_EAME <- nest_model_EAME[which(nest_model_EAME$Error == 'No'), ]

## Check for multicollinearity by comparing r values
panel.cor <- function(x, y, digits = 2, cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  # correlation coefficient
  r <- cor(x, y)
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste("r= ", txt, sep = "")
  text(0.5, 0.6, txt)
  
  # p-value calculation
  p <- cor.test(x, y)$p.value
  txt2 <- format(c(p, 0.123456789), digits = digits)[1]
  txt2 <- paste("p= ", txt2, sep = "")
  if ( p < 0.01) txt2 <- paste("p= ", "<0.01", sep = "")
  text(0.5, 0.4, txt2)
}

## Generate matrix of correlation coefficients to check for multicollinearity
#corr.matrix <- pairs(~ Year  +  , lower.panel = panel.cor, data=nest_model_EAME)

nest_model_EAME$Year = factor(nest_model_EAME$Year)
nest_model_EAME$Species = factor(nest_model_EAME$Species)
nest_model_EAME$FirstFound = as.numeric(nest_model_EAME$FirstFound)


# Write a function for evaluating a set of competing models for the year hypothesis
run.nest_model_EAME.year = function()
{
  # Year
  Year = mark(nest_model_EAME, nocc = occ, model = "Nest", model.parameters = list(S = list(formula = ~Year)), groups = c("Year"))
  return(collect.models())
  
}

nest_model_EAME.results.year = run.nest_model_EAME.year() # This runs the models above and takes a minute or 2
nest_model_EAME.results.year = adjust.chat(1, nest_model_EAME.results.year)
EAME_DSR <- as.data.frame(nest_model_EAME.results.year$Year$results$real)
EAME_DSR$year <- c("2017", "2018", "2019")

```

```{r DSR plot}
#combine all 3 species into a single matrix
DSR <- rbind(GRSP_DSR, DICK_DSR, EAME_DSR)

#rename columns to make them more clear
colnames(DSR) <- c("Estimate", "StandardError", "LCL", "UCL", "Species", "Note",
                   "Year")

#manually add a species column
DSR$Species <- paste(c(" Grasshopper\nSparrow", " Grasshopper\nSparrow", " Grasshopper\nSparrow", " Grasshopper\nSparrow", " Grasshopper\nSparrow", " Grasshopper\nSparrow", "Dickcissel", "Dickcissel", "Dickcissel", "  Eastern\nMeadowlark", "  Eastern\nMeadowlark", "  Eastern\nMeadowlark"))

#make sure that it's classified as a dataframe
DSR <- as.data.frame(DSR)

#add in a group column to separate drought and non-drought years
DSR$group[DSR$Year != "2018"] <- "nondrought"
DSR$group[DSR$Year == "2018"] <- "drought"

#plot all nest survival estimates against one another with upper and lower confidence limits as error bars
DSRplot <- ggplot(DSR) + geom_errorbar(aes(ymin = LCL, ymax = UCL, x = Year), width = .1) + 
  geom_point(aes(x = Year, y = Estimate, fill = group, color = group), size = 4) + 
  theme_custom() + 
  ylab(label = "Daily Survival Rate") + 
  theme(legend.position = "none") + 
  scale_fill_manual(values = c(colors)) +
  facet_grid(cols = vars(Species), scales = "free", space = "free")
DSRplot

DSRsummary <- DSR %>% group_by(group, Species) %>% summarize(DSR = mean(Estimate), se = mean(StandardError))
#extrapolate daily survival rate over entire nesting period

#length of the full nesting period for each species
DSRsummary$length <- c(23,18,18,23,18,18)

#calculate estimated probability of surviving over the entire nesting period
DSRsummary$EstSuccess <- DSRsummary$DSR^DSRsummary$length
```

```{r nest productivity}
#create figure of mean number of fledgelings with 95% CI
nfledgeplot <- ggplot(data = numfledgelings,aes(x = Year,y = nfledge,fill = group, color = group)) + 
  geom_errorbar(aes(ymin = nfledge - 2*se,ymax = nfledge + 2*se,x = Year),color = "black",width = .1,position = position_dodge(width = .9)) + 
  geom_point(size = 4) + 
  theme_custom() + 
  ylab("Mean Number of Young\nFledged per Nest") + 
  theme(legend.position = "none") + 
  scale_fill_manual(values = c(colors)) +
  facet_grid(cols = vars(Species), space = "free", scales = "free")

#combine nest survival and number of young fledged plots
grid.arrange(DSRplot, nfledgeplot, nrow = 2)
```

```{r orientation}
#add column with cardinal directions
#first make all of them "N" since it's the most complicated
nests$Direction <- "N"
#then rename the other rows based on the values in orientation column (45 degrees for each direction)
nests$Direction[nests$NestOrientation > 23 & nests$NestOrientation < 69] <- "NE"
nests$Direction[nests$NestOrientation > 68 & nests$NestOrientation < 114] <- "E"
nests$Direction[nests$NestOrientation > 113 & nests$NestOrientation < 158] <- "SE"
nests$Direction[nests$NestOrientation > 157 & nests$NestOrientation < 203] <- "S"
nests$Direction[nests$NestOrientation > 202 & nests$NestOrientation < 247] <- "SW"
nests$Direction[nests$NestOrientation > 246 & nests$NestOrientation < 292] <- "W"
nests$Direction[nests$NestOrientation > 291 & nests$NestOrientation < 336] <- "NW"

nests$Month <- month(nests$Date1stEgg)
nests <- subset(nests, Species != "Dickcissel" & Month >= 5 & Month <= 7)

#read in wind data and subset to only years of interest
wind <- read.csv("http://lter.konza.ksu.edu/sites/default/files/data/AWE011.csv")
wind <- subset(wind, RECYEAR >= 2014 & RECYEAR <= 2019)

#add group column
wind$group <- "nondrought"
wind$group[wind$RECYEAR == 2018] <- "drought"

#first make all of them "N" since it's the most complicated
wind$Direction <- "N"
#then rename the other rows based on the values in orientation column (45 degrees for each direction)
wind$Direction[wind$WDIR > 23 & wind$WDIR < 69] <- "NE"
wind$Direction[wind$WDIR > 68 & wind$WDIR < 114] <- "E"
wind$Direction[wind$WDIR > 113 & wind$WDIR < 158] <- "SE"
wind$Direction[wind$WDIR > 157 & wind$WDIR < 203] <- "S"
wind$Direction[wind$WDIR > 202 & wind$WDIR < 247] <- "SW"
wind$Direction[wind$WDIR > 246 & wind$WDIR < 292] <- "W"
wind$Direction[wind$WDIR > 291 & wind$WDIR < 336] <- "NW"

#remove unnecessary columns
wind <- wind[c(3,15:16)]

#add in a single dummy value for North in the drought year so that plot works properly
wind_dummy <- as.data.frame(c("2018"))
wind_dummy$group <- "drought"
wind_dummy$Direction <- "N"
names(wind_dummy) <- c("RECYEAR","group","Direction")
wind <- rbind(wind, wind_dummy)

#create individual graphs for each species in drought and nondrought, as well as for ambient wind in drought and nondrought
EAMEnondroughtWind <-  ggplot() + geom_bar(data = subset(nests, Species == "  Eastern\nMeadowlark" & Year != "2018"), aes(x = Direction, y = (..count..)/sum(..count..)), fill = "gray") + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) + 
  ggtitle("EAME") +
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))

EAMEdroughtWind <-  ggplot() + geom_bar(data = subset(nests, Species == "  Eastern\nMeadowlark" & Year == "2018"), aes(x = Direction, y = (..count..)/sum(..count..)), fill = "#FF9999") + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) + 
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))

GRSPnondroughtWind <- ggplot() + geom_bar(data = subset(nests, Species == " Grasshopper\nSparrow" & Year != "2018"), aes(x = Direction, y = (..count..)/sum(..count..)), fill = "gray") + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) + 
  ggtitle("GRSP") +
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))

GRSPdroughtWind <- ggplot() + geom_bar(data = subset(nests, Species == " Grasshopper\nSparrow" & Year == "2018"), aes(x = Direction, y = (..count..)/sum(..count..)), fill = "#FF9999") + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) +
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))

nondroughtWind <- ggplot(data = subset(wind, group == "nondrought")) + geom_bar(aes(x = Direction, y = (..count..)/sum(..count..), fill = "gray")) + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) + 
  ggtitle("Wind") +
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))

droughtWind <- ggplot(data = subset(wind, group == "drought")) + geom_bar(aes(x = Direction, y = (..count..)/sum(..count..), fill = "#FF9999")) + 
  scale_x_discrete(drop = FALSE, labels = c("N", "NE", "E", 
                                       "SE", "S", "SW", "W", "NW"))  + 
  coord_polar(theta = "x", start = -(22.5 * (22/7)/180)) + 
  theme_classic() +
  theme(axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.line = element_blank(),
        text = element_text(size = 20))
  
#combine all individual graphs into a single figure
grid.arrange(EAMEnondroughtWind, GRSPnondroughtWind, nondroughtWind, EAMEdroughtWind, GRSPdroughtWind, droughtWind, ncol = 3)


#split data by species to compare between groups
orientation_GRSP <- subset(nests, Species == " Grasshopper\nSparrow")
orientation_EAME <- subset(nests, Species == "  Eastern\nMeadowlark")

#Watson-Wheeler test of homogeneity of means
#GRSP
orientation_GRSP_drought <- subset(orientation_GRSP, group == "drought" )
orientation_GRSP_nondrought <- subset(orientation_GRSP, group != "drought")

orientation_GRSP_drought_circular <- as.circular(orientation_GRSP_drought$NestOrientation, units = "degrees", rotation = "clock")
orientation_GRSP_nondrought_circular <- as.circular(orientation_GRSP_nondrought$NestOrientation, units = "degrees", rotation = "clock")

#data does not fit the assumptions for watson-williams test, so use a non-parametric watson-wheeler test
watson.wheeler.test(list(orientation_GRSP_drought_circular, orientation_GRSP_nondrought_circular))

#EAME
orientation_EAME_drought <- subset(orientation_EAME, group == "drought")
orientation_EAME_nondrought <- subset(orientation_EAME, group != "drought")

orientation_EAME_drought_circular <- as.circular(orientation_EAME_drought$NestOrientation, units = "degrees", rotation = "clock")
orientation_EAME_nondrought_circular <- as.circular(orientation_EAME_nondrought$NestOrientation, units = "degrees", rotation = "clock")

#data does not fit the assumptions for watson-williams test, so use a non-parametric watson-wheeler test
watson.wheeler.test(list(orientation_EAME_drought_circular, orientation_EAME_nondrought_circular))

```

```{r small mammals}
#import small mammal data
smammals <- read.csv("~/Data/Drought paper/small_mammal_data.csv")

#subset small mammal data to only include the first capture
smammals <- subset(smammals, Recap == "N")
#select only the 4 columns of interest: year, date, grid number, and ear tag number, and rename them to make them more understandable
smammals <- select(smammals, c(1:3,5))
names(smammals) <- c("Year", "Date", "Grid", "ID")
#group by year and grid, then summarize by counting how many individuals were caught in each grid in each year
smammals <- smammals %>% group_by(Year, Grid)
smammalssummary <- summarize(smammals, count = length(ID))

#group summaries by year only and calculate mean, standard deviation, and total number of observations
smammalssummary <- smammalssummary %>% group_by(Year) %>% summarize(mean = mean(count), sd = sd(count), n = sum(count))
#calculate standard error from sd and total n
smammalssummary$se <- smammalssummary$sd/sqrt(smammalssummary$n)

#add group column to distinguish between drought and nondrought years
smammalssummary$group[smammalssummary2$Year != "2018"] <- "nondrought"
smammalssummary$group[smammalssummary2$Year == "2018"] <- "drought"

#plat figure, with errorbars
ggplot() + 
  geom_errorbar(data = smammalssummary, aes(ymin = mean-se, ymax = mean+se, x = Year), width = .1) + 
  geom_point(data = smammalssummary, aes(x = Year, y = mean), size = 5) + 
  theme_custom() +  
  scale_fill_manual(values = c(colors)) +
  ylab(label = "Mean n small mammals")
  
#ANOVA for LTER data
anova <- aov(count ~ Year, data = smammalssummary)
summary(anova)

#calculate effect size from the ANOVA
effectsize(anova)
```

```{r arthropods}
#import data, and create a year column from the DateSampled column
Arthropod_Data <- read.csv("http://lter.konza.ksu.edu/sites/default/files/data/CBS041.csv")
Arthropod_Data$Year <- as.numeric(substr(Arthropod_Data$DateSampled,6,10))

#make vector of unused columns to drop from dataframe
Arthropod_Data <- Arthropod_Data %>% select(c(13, 19))

#add in a group column to separate drought and non-drought years
Arthropod_Data$group[Arthropod_Data$Year != "2018"] <- "nondrought"
Arthropod_Data$group[Arthropod_Data$Year == "2018"] <- "drought"

#separate the drought and nondrought rows then run a t-test to determine if they are different
arthropod_data_nondrought <- subset(Arthropod_Data, group == "nondrought")
arthropod_data_drought <- subset(Arthropod_Data, group == "drought")
t.test(arthropod_data_nondrought$EdibleMass, arthropod_data_drought$EdibleMass)

#calculate means, standard deviation, and standard error
arthropodsummary <- Arthropod_Data %>% group_by(group) %>%
    summarise_each(funs(mean,sd,se=sd(.)/sqrt(n())))
```

```{r arthropod plots}
#boxplot of edible mass for drought vs nondrought
ggplot() +
  geom_boxplot(data = Arthropod_Data, aes(x = group, y = EdibleMass, fill = group), width = .25) +
  theme_custom() +
  labs(y = "Mean Edible Mass")
```

```{r nest summary table}
#calculations for Table 1 in the manuscript

#read in fresh nest data, without the processing used earlier
nests <- read_xlsx("D://Users/inkli/Documents/Data/Drought paper/Nests.xlsx", col_names = TRUE)

#import egg measurement data and remove any false zeroes
eggs <- read_xlsx("D://Users/inkli/Documents/EggMeasurements.xlsx", col_names = TRUE)
eggs <- subset(eggs, EggLength > 0 & EggWidth > 0 & EggMass > 0)

#calculate mean egg length, width, and mass per nest
eggs <- eggs %>% group_by(NestID) %>% dplyr::summarize(length = mean(EggLength), width = mean(EggWidth), mass = mean(EggMass))

#fliter nest data to only those nests where the eggs were measured
nesteggs <- subset(nests, nests[17] == TRUE)

#merge egg measurement data with nest data
eggsmeasured <- merge(nests, eggs, by = "NestID")

#summarize egg measurement data by calculating mean length, width, and mass for each species in each year
eggssummary <- eggsmeasured %>% group_by(Year, Species) %>% dplyr::summarize(length = mean(length), width = mean(width), mass = mean(mass))


#build nest summary table----
#add a year column to the nest data
nests$year <- paste(substr(nests$DateFound, 1,4))

#add in a group column to separate drought and non-drought years
nests$group[nests$year != "2018"] <- "nondrought"
nests$group[nests$year == "2018"] <- "drought"

#calculate max number of cowbird eggs in each nest by subtracting the number of host eggs from the total clutch size in each nest
nests$NBHCOEgg <- nests$TotalClutchSize-nests$MaxNHostEggs

#subset nest data to only include nests where inferred fate is known
nests <- subset(nests, InferredFateField != "no")

#create a new dataframe that only includes the number of hosts fledged, the year, and the species within each nest
NHostFledge <- nests %>% filter(c("NHostFledge", "year", "Species"))

#any nests with an NA for NHostFledge are assumed to be zero
NHostFledge[is.na(NHostFledge)] <- 0

#calculate total Nfledged per species per year
NHostFledge <- NHostFledge %>% group_by(year, species)
NHostFledgeSummary <- NHostFledge  %>% dplyr::summarize(NFledged = sum(NFledged))

#add a column with the total njmber of nests per species per year
NHostFledgeSummary$TotalN <- c(117,107,51,65,51,45,74,67,110,35,52,56)

#calculate average number of hosts fledged per nest by dividing the total number of fledgelings by the total number of nests
NHostFledgeSummary$PerNest <- NHostFledgeSummary$NFledged/NHostFledgeSummary$TotalN

#calculate the standard error for the number of hosts fledged
NHostFledgeSE <- nests %>% group_by(year, Species) %>% dplyr::summarize(se = sd(NHostFledge)/nrow(NHostFledge))

#calculate the number of hosts fledged per nest for each species, but only for nests that were successful
NHostFledgeSU <- subset(nests, InferredFateField == "SU")
NHostFledgeSU <- NHostFledgeSU %>% group_by(year, Species)
#separate out each species
NHostFledgeSUGRSP <- subset(NHostFledgeSU, Species == "GRSP")
NHostFledgeSUDICK <- subset(NHostFledgeSU, Species == "DICK")
NHostFledgeSUEAME <- subset(NHostFledgeSU, Species == "EAME")
#calculate total number of successful nests for each species
tGRSP <- as.data.frame(table(NHostFledgeSUGRSP$year))
tDICK <- as.data.frame(table(NHostFledgeSUDICK$year))
tEAME <- as.data.frame(table(NHostFledgeSUEAME$year))
#merge all 3 species back inot a single dataframe
t <- rbind(tGRSP, tDICK, tEAME)

#add the total number of successful nests to the summary table created above
NHostFledgeSummary$NSuccess <- t$Freq

#add a column with the number of young fledged per nest
NHostFledgeSummary$NFledged <- c(120,85,38,39,46,54,35,36,82,36,50,64)
NHostFledgeSummary$PerSUNest <- NHostFledgeSummary$NFledged/NHostFledgeSummary$NSuccess

#get nest initiation data----
#obtain Julian dates that the first egg was laid per nest
nests$JulianDate <- as.numeric(yday(nests$Date1stEgg))
#separate by species to obtain the total range of dates for each species
EAMEdata <- as.data.frame(subset(nests, Species == "EAME")) #range = 82, 11.7 weeks
GRSPdata <- as.data.frame(subset(nests, Species == "GRSP")) #range = 95, 13.5 weeks
DICKdata <- as.data.frame(subset(nests, Species == "DICK")) #range = 67, 9.5 weeks

#Count of partial clutch loss and partial brood loss per total N nests
#filter to only the nests that experienced partial clutch loss
clutch <- nests %>% filter(`PartialEggLoss?` == TRUE) %>% select(c(4,14,40))

#calculate total number of nests per year that experienced partial clutch loss
clutch <- clutch %>% group_by(Species, year) %>% dplyr::summarize(clutchloss = sum(clutchloss))

#do the same calculations as above but for partial brood loss
brood <- nests %>% filter(`PartialBroodLoss?` == TRUE) %>% select(c(4,14,43))
brood <- brood %>% group_by(Species, year) %>% dplyr::summarize(broodloss = sum(broodloss))

#obtain the species, year, and total number of nests from the summary table above
partialloss <- select(NHostFledgeSummary,c(2,1,4))

#merge the species, year, and total number of nests columns with the clutch and brood loss summaries (must do this in two separate steps as merging three dataframes at once doesn't work)
partialloss <- merge(partialloss, clutch)
partialloss <- merge(partialloss, brood)

#calculate clutch and brood loss per nest by dividing each by the total number of nests
partialloss$clutchlosspernest <- partialloss$clutchloss/partialloss$TotalN
partialloss$broodlosspernest <- partialloss$broodloss/partialloss$TotalN

#add a group column to distinguish between drought and nondrought years
partialloss$group[partialloss$year != "2018"] <- "nondrought"
partialloss$group[partialloss$year == "2018"] <- "drought"

#summarize partial clutch and brood loss by species and by drought/nondrought, calculating mean and standard deviation of each
partialloss <- partialloss %>% group_by(Species, group) %>% summarize(clutchlosssd = sd(clutchlosspernest), broodlosssd = sd(broodlosspernest), clutchlosspernest = mean(clutchlosspernest), broodlosspernest = mean(broodlosspernest))
```
