### WEST SISTER ISLAND NEST COUNTS #################################################
# Joshua Booker
# 21 Aug 2018, updated 2/3/2026 by Jessica Duez
# This analysis is adapted from Shieldcastle's SAS code, received 6/6/2016 
###################################################################################

#STEP 1: Create nest count raw data spreadsheet (csv) with only plot, year, and species for desired year
#STEP 2: run the code below (highlight, CTRL+Enter)
#STEP 3: update report from last year (WSIsummary2025.doc)
#STEP 4: Send population estimates to OHIO DNR
#STEP 5: upload data and report to ServCat (https://iris.fws.gov/APPS/ServCat/

#-------importing and inspecting data-----------------
library(readxl)
library(tidyr)
library(dplyr)

WSIdata <- read.csv("S:/biology/Birds/west sister island/nest counts/WSI_NestCount2025.csv")

###tidy and analyze the data########################

#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  
         
###Be sure to rename WSIdata_summary2025 with proper year before running the following code, as it will overwrite the existing file!###

write.csv(WSIdata_summary, "S:/Biology/Birds/west sister island/nest counts/WSIsummarys/WSIdata_summary2025.csv")