### WEST SISTER ISLAND NEST COUNTS ################################################# # Joshua Booker # 21 Aug 2018, updated 10/28/2021 # This analysis is adapted from Shieldcastle's SAS code, received 6/6/2016 ################################################################################### #-------importing and inspecting data----------------- library(readxl) WSIdata <- read.csv("S:/biology/Birds/west sister island/nest counts/WSIdata_allyrs_corrected.csv") #In 2019, data entry errors were founded and corrected, based on original paper data sheets #import new year's data (replace year as needed) nestcounts2024 <- read_excel("S:/biology/Birds/west sister island/nest counts/WSINWR 2024 nest counts.xlsx") #import original legacy estimates from H:\OhioNWRS\My Stations\Ottawa\West Sister Island nest counts\west sister breeding pairs 2016.pdf. This data is uncorrected and doesn't match legacy extrapolation below legacy_join <- read.csv("S:/biology/Birds/west sister island/nest counts/legacy_WSI_summary1991-2016.csv") #-----------------tidy the data---------------------------------- library(tidyr) library(dplyr) #renaming columns (replace year as needed) colnames(nestcounts2024)[1] <- "YEAR" colnames(nestcounts2024)[2] <- "PLOT" colnames(nestcounts2024)[13] <- "Canopy" colnames(nestcounts2024)[17] <- "Survey.notes" #remove plots not surveyed (*update year!*) nestcounts2024 <- nestcounts2024[!is.na(nestcounts2024$GBHE), ] write.csv(nestcounts2024, "S:/biology/Birds/west sister island/nest counts/nestcounts2024.csv") write.csv(nestcounts2024, "C:/Users/jbooker/Desktop") #add to collective dataset and export nestcounts2024$HerbCover <- as.numeric(nestcounts2024$HerbCover) nestcounts2024$ShrubCover <- as.numeric(nestcounts2024$ShrubCover) nestcounts2024$Canopy <- as.numeric(nestcounts2024$Canopy) WSIdata$X <- NULL WSIdata <- rbind(WSIdata, nestcounts2024) #re-format columns as needed to bind WSIdata <- WSIdata[-c(6675), ] #removing totals row write.csv(WSIdata, "S:/biology/Birds/west sister island/nest counts/WSIdata_allyrs_corrected.csv") #count # of plots in each year count <- count(WSIdata, YEAR) #tidy corrected counts WSIdata <- as_tibble(WSIdata) WSIdata_gathered <- WSIdata %>% gather(Species, NestCount, GBHE:DCCO, na.rm=TRUE) #combining the 8 species columns into a single column #Estimating population using Mark's extrapolation (15% of island) and corrected extrapolation (10% of island) and corrected SE (5/23/2019)-------------------------- WSIdata_gathered %>% group_by(Species, YEAR) %>% #grouping by species and year summarize(NumPtsSampled=n(), TotalCount=sum(NestCount), MeanCount=mean(NestCount, na.rm=TRUE), SD_sample=sd(NestCount, na.rm=TRUE), SE_sample=(sd(NestCount, na.rm=TRUE))/sqrt(n()), #uncorrected fpc=sqrt(((n()*10)-n())/((n()*10)-1)), #fpc=finite population correction formula. assuming N=n*10 SE_mean_corrected=SE_sample*fpc, CI95_lower_mean=MeanCount-(2.0452*SE_mean_corrected), #confidence intervals for the pop mean CI95_higher_mean=MeanCount+(2.0452*SE_mean_corrected), Pop_extrapolated=TotalCount/0.1) -> WSIdata_summary #extrapolating the Total count, assuming that 10% is counted each year no matter how many points were actually visited WSIdata_summary %>% mutate(pop_extrapolated_legacy=((MeanCount*196)*100)/15, #calculating the population estimate using the mean and 15% of island (legacy method). this controls for the differing # of points counted each year SE_pop_legacy=((SE_sample*196)*100)/10, #legacy method didn't use FPC Pop_extrapolated2=((MeanCount*196)*100)/10, #calculating the population estimate using the mean and 10% of island. this controls for the differing # of points counted each year SE_pop2=((SE_mean_corrected*196)*100)/10, CI95_lower_pop=Pop_extrapolated2-(2.0452*SE_pop2), #confidence intervals for the pop total (10% extrapolation) CI95_higher_pop=Pop_extrapolated2+(2.0452*SE_pop2)) -> WSIdata_summary write.csv(WSIdata_summary, "H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/WSIdata_summary2.csv") WSIdata_summary <- read.csv("H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/WSIdata_summary2.csv") #importing and joining GIS interp estimates##################################### #IDW estimates from current year need to be added manually to the spreadsheet below IDW_summary <- read.csv("H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/GIS/IDWsummary.csv") #Formatting to join IDW with extrap IDW_summary$Species <- as.character(IDW_summary$Species) WSIdata_simple <- WSIdata_summary[ ,-c(4:10)] #OPTIONAL: might only want Year, Species, pop estimates WSIdata_simple$X <- NULL #removing extra column WSIdata_simple$Species <- as.character(WSIdata_simple$Species) #Joining allpops <- inner_join(WSIdata_simple, IDW_summary) allpops$X <- NULL allpops$X.1 <- NULL write.csv(allpops, "H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/GIS/allpops.csv") allpops20192020 <- subset(allpops, YEAR==2019 | YEAR==2020) #PLOT TRENDS---------------------------------------------------------- #Canopy cover WSIdata_gathered %>% group_by(YEAR) %>% #grouping by year summarize(AverageCanopy=mean(Canopy, na.rm=TRUE)) -> WSIdata_Canopy plot(WSIdata_Canopy) #declining trend? #still need to break this data down into tidy format allpop_gathered <- allpops %>% gather(Method, value, c(Pop_extrapolated2, IDW)) allpop_gathered <- allpop_gathered %>% rename(Breeding_Pairs=value) #rename "value" to Breeding Pairs allpop_gathered$SE_pop2[allpop_gathered$Method=="IDW"] <- NA #deleting SE and CI for IDW estimates allpop_gathered$CI95_lower_pop[allpop_gathered$Method=="IDW"] <- NA allpop_gathered$CI95_higher_pop[allpop_gathered$Method=="IDW"] <- NA allpop_gathered$CI95_lower_pop <- ifelse(allpop_gathered$CI95_lower_pop < 0 & !is.na(allpop_gathered$CI95_lower_pop), 0, allpop_gathered$CI95_lower_pop) #need to change negative values to 0 write.csv(allpop_gathered, "H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/GIS/allpops_gathered.csv") allpop_gathered <- read.csv("H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/GIS/allpops_gathered.csv") library(ggplot2) library(readr) #extrapolated estimates trend ex_sub <- subset(allpop_gathered, Method=="Pop_extrapolated2") ggplot(ex_sub, aes(x=YEAR, y=Breeding_Pairs, col=Species))+ geom_line(aes(group=Species))+geom_point() #With DCCO management DCCOmgmt <- read_excel("H:/OhioNWRS/My Stations/Ottawa/West Sister Island nest counts/cormorant take under PRDO summary.xls", sheet = "WSI summary") #importing file DCCOmgmt<-DCCOmgmt[-c(29:32), ] #removing extra rows DCCOmgmt$DCCOpop <-NULL #removing pop numbers DCCOmgmt %>% mutate(Species="DCCO") -> DCCOmgmt ggplot(ex_sub, aes(x=YEAR, y=Breeding_Pairs, col=Species))+ #color by Species geom_line(aes(group=Species))+ #add line for each species geom_point(aes(shape=Species))+ #add diff shaped point for each species scale_shape_manual(values=c(0, 15, 1, 16, 2, 17, 4, 8))+ #setting point shapes manually scale_x_continuous(breaks=c(1991:2020)) #update year if neccessary ggplot(DCCOmgmt, aes(x=Year, y=Take, group=1))+ geom_line(size=2)+ geom_point(size=4) #plot IDW trend IDW_sub <- subset(allpop_gathered, Method=="IDW") ggplot(IDW_sub, aes(x=YEAR, y=Breeding_Pairs, col=Species))+ geom_line(aes(group=Species))+geom_point() #calculating % change in pop est between IDW and corrected extrapolation methods------------------------- IDW_sub <- as.tibble(IDW_sub) IDW_sub <- add_column(IDW_sub, change=((IDW_sub$Breeding_Pairs - ex_sub$Breeding_Pairs) / ex_sub$Breeding_Pairs)*100) simple_IDW_sub <- IDW_sub[, c(2:4, 13)] #extracting relevant columns #comparing mean change of past 2 years mean(simple_IDW_sub$change[simple_IDW_sub$YEAR==2020], na.rm=TRUE) mean(simple_IDW_sub$change[simple_IDW_sub$YEAR==2019], na.rm=TRUE) #checking relationship between sample size and change plot(simple_IDW_sub$NumPtsSampled, simple_IDW_sub$change) lm(simple_IDW_sub$change~simple_IDW_sub$NumPtsSampled) #regression abline(-99.9555, 0.4247) #plotting regression line. no strong relationship #checking if any IDW estimates fall outside of extrapolation confidence interval. Nope! (8/2019) join <- inner_join(ex_sub, IDW_sub, by=c("Species", "YEAR")) check <- subset(join, join$Breeding_Pairs.y>join$CI95_higher_pop.x|join$Breeding_Pairs.y