###DATA MANIPULATION CODE###
#Code for reproducing figures and tables from de-identified angler survey data
#Drum Repository:
#Paper: 
#Date last updated:8/20/2021

#Note: Please see data manipulation.R for non-reproducible data cleaning and analysis process. Below is the code necessary for replicating figures from the provided de-identified datasets.

## Reproducible code for figures using deidentified data:

#Load packages
require(data.table)
require(tidyverse)

#Load data
#Note: change file path to appropriate path on your computer, or set the working directory to the folder in which the files are saved
survey.responses<-read.csv("Mailed_survey_data_de-identified.csv")
survey.responses.baitusers<-read.csv("Mailed_survey_data_de-identified_baitusers.csv")
##Data preparation

#For bait users, categorize the method of disposal
survey.responses.baitusers$disposal_method<- ifelse(survey.responses.baitusers$Q16_1_release==TRUE,"Release in water",
ifelse(survey.responses.baitusers$Q16_2_dispose==TRUE,"Properly disposed",
ifelse(survey.responses.baitusers$Q16_3_giveaway==TRUE, "Gave away",
ifelse(survey.responses.baitusers$Q16_4_other==TRUE,"Other",
ifelse(survey.responses.baitusers$Q16_5_noneleft==TRUE,"None leftover",                                                                               
       
      
  #code for harvest disposal methods:
ifelse(survey.responses.baitusers$Q9=="Released into different water from which it was caught","Release in water",
ifelse(survey.responses.baitusers$Q9=="Released into the same water from which it was caught", "Back to source",
ifelse(survey.responses.baitusers$Q9=="Disposed on land/in garbage","Properly disposed",
ifelse(survey.responses.baitusers$Q9=="Gave to a friend", "Gave away",
ifelse(survey.responses.baitusers$Q9=="Other (please describe)","Other",
ifelse(survey.responses.baitusers$Q9=="I did not have any leftover baitfish","None leftover",

  #code for no response
ifelse(survey.responses.baitusers$Q9==""&survey.responses.baitusers$Q16=="","Did not respond", ""))))))))))))

## Code for reproducing Figure 2:
Figure2<- survey.responses.baitusers %>%
  mutate(admin = fct_relevel(admin, "Central", "South","Northeast","Northwest", "out of state")) %>%
  mutate(disposal_method = factor(disposal_method, levels=c("Release in water", "Properly disposed","Gave away","Back to source", "None leftover", "Other", "Did not respond"))) %>%
  ggplot( aes(x=admin, fill=disposal_method)) +
  geom_bar(position = "fill")  + xlab("Region" ) + ylab("proportion") +
  theme_bw()+
  theme(axis.text.x = element_text(angle=0,vjust = .5)) +
  theme(axis.text = element_text(size=15)) +
  theme(axis.title = element_text(size=20))+
  theme(legend.text = element_text(size = 15)) +
  theme(legend.title = element_text(size = 20)) +
  scale_fill_grey()+
  labs(fill="Disposal method")
Figure2

##Code for reproducing Figure 3:
median(survey.responses.baitusers$trips.per.year)
Figure3<- survey.responses.baitusers %>% 
     mutate(Q16_1_release = replace(Q16_1_release, Q16_1_release == FALSE, "no")) %>%
     mutate(Q16_1_release = replace (Q16_1_release, Q16_1_release == TRUE, "yes")) %>%
     ggplot(aes(x=trips.per.year,fill=Q16_1_release)) + 
     geom_histogram(bins = 100) + stat_bin(bins = 50) + 
     ggtitle("Number of trips per angler per year") +
     xlab("Number of trips")+
     theme_bw()+
     geom_vline(xintercept =  20) +
     theme(text = element_text(size=24), axis.text.x = element_text(size=24), 
           axis.text.y = element_text(size=24),
           legend.title = element_text(size = 24))+
     labs(fill="Release")+
     #scale_fill_discrete(name = "Release", labels = c("no","yes"))+
     scale_fill_grey()

Figure3

##Code for reproducing Table 1:
table(survey.responses.baitusers%>%group_by(live.bait, Q16) %>% summarize(n=n()))

##Code for reproducing Table 2: 

#Calculate weights for each age 
#weight=1/(sampleprop/populationprop)
weight1839<-1/(0.45/0.66)
weight4049<-1/(0.094/0.086)
weight5059<-1/(0.14/0.09)
weight6069<-1/(0.19/0.10)
weight70plus<-1/(0.11/0.05)
weightsarray<-c(weight1839,weight4049,weight5059,weight6069,weight70plus)

#Weighted release rate
sbdt <- data.table(survey.responses.baitusers) 

sbdt[, .N, by = .(age_cat)] 
sbdt[ , .N, by = .(age_cat, Q16_1_release)] 
sbdt2 <- merge(sbdt[ , .N, by = .(age_cat, Q16_1_release)],sbdt[, .N, by = .(age_cat)], by = "age_cat") #Creates one dataframe with number of respondents by disposal method per location as well as total number of respondents per location
sbdt2[,prop := N.x/N.y , ]

release.by.age<-sbdt2[Q16_1_release==TRUE,prop] #array of proportions of anglers who released per age category

#weighted average release rate: 
sum(release.by.age*weightsarray)/sum(weightsarray)

#weighted use of live bait
survey.responses %>% group_by(age_cat,live.bait) %>% summarize (n=n())
srdt<-data.table(survey.responses)
srdt2<- merge(srdt[ , .N, by = .(age_cat, live.bait)],srdt[, .N, by = .(age_cat)], by = "age_cat") #Creates one dataframe with number of respondents by disposal method per location as well as total number of respondents per location
srdt2[,prop := N.x/N.y , ]
srdt2

use.by.age<-srdt2[live.bait=="Did not use live bait",prop] #array of proportions of anglers who did not use live bait by age
sum(use.by.age*weightsarray)/sum(weightsarray)