##%######################################################%## # # #### Demographic analysis of PNG sharks #### # # ##%######################################################%## # Mon Apr 20 09:14:59 2020 ------------------------------ # Dr Jonathan Smart # SARDI Aquatic Sciences # Adelaide, Australia # jonathan.smart@sa.gov.au library(tidyverse) library(SharkDemography) library(directlabels)# to add text to contours library(cowplot) ### Load previously created datafiles for each species ------------------------ silky_data <- readRDS("silky.RDS") OCS_data <- readRDS("OCS.RDS") AML_data <- readRDS("AML.RDS") silvertip_data <- readRDS("silvertip.RDS") limbatus_data <- readRDS("limbatus.RDS") Tiger_shark_data <- readRDS("Tiger.RDS") Bigeye_thresher_data <- readRDS("BigEye.RDS") Mako_data <- readRDS("mako.RDS") hammerhead_data <- readRDS("hammehead.RDS") Spinner_shark_data <- readRDS("Spinner.RDS") Blue_shark_data <- readRDS("Blue.RDS") Pelagic_thresher_data <- readRDS("PelagicThresher.RDS") ### Select appropriate M estimators to be used for each species Silky_M_ests <- c("Chen.Yuan", "Jensen.mat") OCS_M_ests <- c("Pet.Wro","Jensen.mat","Chen.Yuan", "Then_pauly", "Jensen", "Chen.Want") Silvertip_M_ests <- c("Pet.Wro","Jensen" "Jensen.mat","Chen.Want") limbatus_M_ests <- c("Pet.Wro","Chen.Yuan", "Then_pauly", "Jensen","Chen.Want") Tiger_M_ests <- c("Pet.Wro","Jensen.mat","Chen.Yuan", "Then_hoenig","Then_pauly", "Jensen","Chen.Want") BigEye_M_ests <- c("Pet.Wro","Jensen.mat","Chen.Yuan", "Then_pauly", "Chen.Want") Mako_M_ests <- c("Chen.Yuan","Jensen", "Jensen.mat","Then_pauly") Spinner_shark_M_ests <- c("Pet.Wro", "Jensen.mat", "Then_pauly") Pelagic_thresher_M_ests <- c("Jensen.mat","Then_pauly") ### Conduct standard analysis for each species-------------- # perform analyses silky_results <- Simulate_demography(10000, silky_data, M.estimators = Silky_M_ests) OCS_results <- Simulate_demography(10000, OCS_data, M.estimators = OCS_M_ests) AML_results <- Simulate_demography(10000, AML_data) silvertip_results <- Simulate_demography(10000, silvertip_data, M.estimators = Silvertip_M_ests) limbatus_results <- Simulate_demography(10000, limbatus_data, M.estimators = limbatus_M_ests) Tiger_shark_results <- Simulate_demography(10000, Tiger_shark_data, M.estimators = Tiger_M_ests) Bigeye_thresher_results <- Simulate_demography(10000, Bigeye_thresher_data, M.estimators = BigEye_M_ests) Mako_results <- Simulate_demography(10000, Mako_data, M.estimators = Mako_M_ests) hammerhead_shark_results <- Simulate_demography(10000, hammerhead_data) Spinner_shark_results <- Simulate_demography(10000, Spinner_shark_data, M.estimators = Spinner_shark_M_ests) Blue_shark_results <- Simulate_demography(10000, Blue_shark_data) Pelagic_thresher_shark_results <- Simulate_demography(10000, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests) # create ouput table for MS create_results_table <- function(x){ x <- as.data.frame(x$`MonteCarlo_summary`) %>% rownames_to_column("Parameter") y <- x %>% transmute(results = paste0(round(Mean,2)," (",round(`.025`,2)," - ", round(`.975`,2),")")) %>% t() %>% as.data.frame() %>% set_names(x$Parameter) %>% rownames_to_column() %>% .[-1] return(y) } Results_table <- bind_rows(create_results_table(silky_results) %>% mutate(Species = "Silky sharks"), create_results_table(OCS_results) %>% mutate(Species = "Oceanic whitetip sharks"), create_results_table(AML_results) %>% mutate(Species = "Grey reef sharks"), create_results_table(silvertip_results) %>% mutate(Species = "Silvertip sharks"), create_results_table(limbatus_results) %>% mutate(Species = "Common blacktip sharks"), create_results_table(Tiger_shark_results) %>% mutate(Species = "Tiger sharks"), create_results_table(Bigeye_thresher_results) %>% mutate(Species = "Big eye thesher sharks"), create_results_table(Mako_results) %>% mutate(Species = "Shortfin mako sharks"), create_results_table(hammerhead_shark_results) %>% mutate(Species = "Scalloped hammerhead sharks"), create_results_table(Spinner_shark_results) %>% mutate(Species = "Spinner sharks"), create_results_table(Blue_shark_results) %>% mutate(Species = "Blue sharks"), create_results_table(Pelagic_thresher_shark_results) %>% mutate(Species = "Pelagic thresher sharks")) %>% arrange(Species) %>% .[c(9,1:8)] # xlsx::write.xlsx(Results_table,file = "Data/Demography results.xlsx", row.names = FALSE) # combine simulations in a dataset and filter out anything beyound the 99th percentiles to get rid of outliers combined_sims <- bind_rows(silky_results$simulations %>% mutate(Species = "Silky sharks"), OCS_results$simulations %>% mutate(Species = "Oceanic whitetip sharks"), AML_results$simulations %>% mutate(Species = "Grey reef sharks"), silvertip_results$simulations %>% mutate(Species = "Silvertip sharks") %>% filter(Lambda < 1.3), # Silvertips sometimes get inflated results due to wide error on growth pars. remove these for clean dists limbatus_results$simulations %>% mutate(Species = "Common blacktip sharks"), Tiger_shark_results$simulations %>% mutate(Species = "Tiger sharks"), Bigeye_thresher_results$simulations %>% mutate(Species = "Big eye thesher sharks"), Mako_results$simulations %>% mutate(Species = "Shortfin mako sharks"), hammerhead_shark_results$simulations %>% mutate(Species = "Scalloped hammerhead sharks"), Spinner_shark_results$simulations %>% mutate(Species = "Spinner sharks"), Blue_shark_results$simulations %>% mutate(Species = "Blue sharks"), Pelagic_thresher_shark_results$simulations %>% mutate(Species = "Pelagic thresher sharks")) %>% select(-M.estimator) %>% gather(Parameter, Value, -Species) %>% group_by(Parameter, Species) %>% filter(quantile(Value, 0.99,na.rm = TRUE) > Value, quantile(Value, 0.01,na.rm = TRUE) < Value) # rename factor levels for nice plot labels combined_sims$Parameter <- factor(combined_sims$Parameter, levels = c("Lambda", "R0", "G", "elast.fecund", "elast.juv.survival", "elast.adult.survival", "juv.ratio", "adult.ratio"), labels = c("lambda", "R[0]", "G", "Fecundity~elasticity", "Juvenile~survival~elasticity", "Adult~survival~elasticity", "Juvenile~survival~elasticity~ratio", "Adult~survival~elasticity~ratio")) # Change to Latin names combined_sims$Species <-factor(combined_sims$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) # main plot p1 <- ggplot(combined_sims %>% filter(Parameter %in% c("lambda", "R[0]", "G")) %>% filter(Parameter != "R[0]" | !Species %in% c("Prionace glauca", "Sphyrna lewini")), aes(x = Species, y = Value, fill = Species)) + geom_hline(data = data.frame(Parameter = "lambda", Species = unique(combined_sims$Species), Value = 1), aes(fill = NULL, yintercept = Value), col = "red", size = 1, linetype = "dashed")+ geom_violin(scale = "width", width = 0.7) + stat_summary(fun.y = mean, col = "magenta", geom="point")+ facet_wrap(~Parameter, ncol = 1,strip.position = "left", scales = "free_y", labeller = label_parsed) + scale_fill_viridis_d()+ expand_limits(y = 0.9)+ labs(y = NULL, x = NULL)+ guides(fill = FALSE)+ theme_bw()+ theme(axis.text.x = element_text(angle = 60, vjust = 1, hjust=1, size = 10,face = "italic"), strip.background = element_rect(fill = "white", colour = NA), strip.placement = "outside") # Blue sharks need their own axis for R0 which will be added over the top of the main plot p2 <- ggplot(combined_sims %>% filter(Parameter %in% c("lambda", "R[0]", "G")) %>% filter(Parameter == "R[0]", Species == "Prionace glauca"), aes(x = Species, y = Value, fill = Species)) + geom_violin(scale = "width", width = 0.7) + stat_summary(fun.y = mean, col = "magenta", geom="point")+ scale_fill_viridis_d(begin = 11/12)+ labs(y = NULL, x = NULL)+ scale_y_continuous(limits = c(0,300), expand = c(0,0))+ guides(fill = FALSE)+ annotate(geom = "line", x = -Inf, y = c(-Inf, Inf))+ theme(axis.text.x = element_blank(), strip.background = element_rect(fill = "white", colour = NA), strip.placement = "outside", panel.background = element_blank(), axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.text = element_text(size = 7)) # Same applies for scalloped hammerheads p3 <- ggplot(combined_sims %>% filter(Parameter %in% c("lambda", "R[0]", "G")) %>% filter(Parameter == "R[0]", Species == "Sphyrna lewini"), aes(x = Species, y = Value, fill = Species)) + geom_violin(scale = "width", width = 0.7) + stat_summary(fun.y = mean, col = "magenta", geom="point")+ scale_fill_viridis_d(begin = 12/12)+ labs(y = NULL, x = NULL)+ scale_y_continuous(limits = c(0,80), expand = c(0,0))+ guides(fill = FALSE)+ annotate(geom = "line", x = -Inf, y = c(-Inf, Inf))+ theme(axis.text.x = element_blank(), strip.background = element_rect(fill = "white", colour = NA), strip.placement = "outside", panel.background = element_blank(), axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.text = element_text(size = 7)) # These need to be altered as plots are produced so that blues and hammers align up # with their axes p4 <- ggdraw() + draw_plot(p1, 0, 0, 1, 1)+ draw_plot(p2, width = 0.12, height = 0.25, x = .795, y = .464) + draw_plot(p3, width = 0.12, height = 0.25, x = .87, y = .464) p4 # export ggsave("standard.demographic.analysis.plot.png",p4, width = 7, height = 8.5, units = "in") # mortality plot for appendix mort_names <- c( 'Chen.Yuan' = "Chen and Yuan (2006)", 'Jensen.mat' = "Jensen (1996) incorporating age-at-maturity", 'Pet.Wro' = "Petersen and Wroblewski (1984)", 'Then_hoenig' = "Hoenig (1983) - Updated by Then et al. (2015)" ) silky_mort_dists<- Simulate_demography(10000, silky_data) # # Code for checking mort dists of lambda ggplot(OCS_results$simulations, aes(Lambda)) + geom_histogram() + geom_vline(xintercept = 1, col = "red", linetype = "dashed", size = 1)+ labs(x = expression(lambda), y = "Count")+ facet_wrap(~M.estimator, nrow = 4)+ theme_bw() #### Determine and plot eigenvectors ----------------------------------------------------- # perform analyses silky_eigenvectors <- Estimate_eigenvectors(1000, silky_data, M.estimators = Silky_M_ests) OCS_eigenvectors <- Estimate_eigenvectors(1000, OCS_data, M.estimators = OCS_M_ests) AML_eigenvectors <- Estimate_eigenvectors(1000, AML_data) silvertip_eigenvectors <- Estimate_eigenvectors(1000, silvertip_data, M.estimators = Silvertip_M_ests) limbatus_eigenvectors <- Estimate_eigenvectors(1000, limbatus_data, M.estimators = limbatus_M_ests) Tiger_shark_eigenvectors <- Estimate_eigenvectors(1000, Tiger_shark_data, M.estimators = Tiger_M_ests) Bigeye_thresher_eigenvectors <- Estimate_eigenvectors(1000, Bigeye_thresher_data, M.estimators = BigEye_M_ests) Mako_eigenvectors <- Estimate_eigenvectors(1000, Mako_data, M.estimators = Mako_M_ests) hammerhead_shark_eigenvectors <- Estimate_eigenvectors(1000, hammerhead_data) Spinner_shark_eigenvectors <- Estimate_eigenvectors(1000, Spinner_shark_data, M.estimators = Spinner_shark_M_ests) Blue_shark_eigenvectors <- Estimate_eigenvectors(1000, Blue_shark_data) Pelagic_thresher_shark_eigenvectors <- Estimate_eigenvectors(1000, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests) combined_eigenvectors <- bind_rows(silky_eigenvectors %>% mutate(Species = "Silky sharks"), OCS_eigenvectors %>% mutate(Species = "Oceanic whitetip sharks"), AML_eigenvectors %>% mutate(Species = "Grey reef sharks"), silvertip_eigenvectors %>% mutate(Species = "Silvertip sharks"), limbatus_eigenvectors %>% mutate(Species = "Common blacktip sharks"), Tiger_shark_eigenvectors %>% mutate(Species = "Tiger sharks"), Bigeye_thresher_eigenvectors %>% mutate(Species = "Big eye thesher sharks"), Mako_eigenvectors%>% mutate(Species = "Shortfin mako sharks"), hammerhead_shark_eigenvectors %>% mutate(Species = "Scalloped hammerhead sharks"), Spinner_shark_eigenvectors %>% mutate(Species = "Spinner sharks"), Blue_shark_eigenvectors %>% mutate(Species = "Blue sharks"), Pelagic_thresher_shark_eigenvectors %>% mutate(Species = "Pelagic thresher sharks")) # Change to Latin names combined_eigenvectors$Species <-factor(combined_eigenvectors$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) Stable_age_plot <- ggplot(combined_eigenvectors %>% filter(Type == "Stable Age Distribution"), aes(x = Age, y = mean, ymin = low, ymax = upp, fill = Species)) + geom_ribbon(alpha = .8)+ geom_line(size = 1, col = "black")+ scale_fill_viridis_d(guide = FALSE) + labs(y = "Proportion of population", x = "Age (years)")+ facet_wrap(~Species, ncol = 3, scales = "free")+ theme_bw() + theme(strip.text = element_text(face = "italic")) ggsave("Stable age distribution.png",Stable_age_plot, width = 7, height = 7, units = "in") Repro_val_plot <- ggplot(combined_eigenvectors %>% filter(Type != "Stable Age Distribution"), aes(x = Age, y = mean, ymin = low, ymax = upp, fill = Species)) + geom_ribbon(alpha = .8)+ geom_line(size = 1, col = "black")+ scale_fill_viridis_d(guide = FALSE) + labs(y = "Reproductive Value", x = "Age (years)")+ scale_y_continuous(breaks = scales::pretty_breaks())+ facet_wrap(~Species, ncol = 3, scales = "free")+ theme_bw() + theme(strip.text = element_text(face = "italic")) ggsave("Reproductive value.png",Repro_val_plot, width = 7, height = 7, units = "in") #### F critical analysis ----------------------------------------------------------------- # perform analyses silky_F_Crit <- Simulate_F_critical(1000, silky_data, M.estimators = Silky_M_ests, max.F = 0.6) OCS_F_Crit <- Simulate_F_critical(1000, OCS_data, M.estimators = OCS_M_ests, max.F = 0.6) AML_F_Crit <- Simulate_F_critical(1000, AML_data, max.F = 0.6) silvertip_F_Crit <- Simulate_F_critical(1000, silvertip_data, M.estimators = Silvertip_M_ests, max.F = 0.6) limbatus_F_Crit <- Simulate_F_critical(1000, limbatus_data, M.estimators = limbatus_M_ests, max.F = 0.6) Tiger_shark_F_Crit <- Simulate_F_critical(1000, Tiger_shark_data, M.estimators = Tiger_M_ests, max.F = 0.6) Bigeye_thresher_F_Crit <- Simulate_F_critical(1000, Bigeye_thresher_data, M.estimators = BigEye_M_ests, max.F = 0.6) Mako_F_Crit <- Simulate_F_critical(1000, Mako_data, M.estimators = Mako_M_ests, max.F = 0.6) hammerhead_shark_F_Crit <- Simulate_F_critical(1000, hammerhead_data, max.F = 0.6) Spinner_shark_F_Crit <- Simulate_F_critical(1000, Spinner_shark_data, M.estimators = Spinner_shark_M_ests, max.F = 0.6) Blue_shark_F_Crit <- Simulate_F_critical(1000, Blue_shark_data, max.F = 0.6) Pelagic_thresher_shark_F_Crit <- Simulate_F_critical(1000, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests, max.F = 0.6) combined_F_crit_sims <- bind_rows(silky_F_Crit$Simulations %>% mutate(Species = "Silky sharks"), OCS_F_Crit$Simulations %>% mutate(Species = "Oceanic whitetip sharks"), AML_F_Crit$Simulations %>% mutate(Species = "Grey reef sharks"), silvertip_F_Crit$Simulations %>% mutate(Species = "Silvertip sharks"), limbatus_F_Crit$Simulations %>% mutate(Species = "Common blacktip sharks"), Tiger_shark_F_Crit$Simulations %>% mutate(Species = "Tiger sharks"), Bigeye_thresher_F_Crit$Simulations %>% mutate(Species = "Big eye thesher sharks"), Mako_F_Crit$Simulations %>% mutate(Species = "Shortfin mako sharks"), hammerhead_shark_F_Crit$Simulations %>% mutate(Species = "Scalloped hammerhead sharks"), Spinner_shark_F_Crit$Simulations %>% mutate(Species = "Spinner sharks"), Blue_shark_F_Crit$Simulations %>% mutate(Species = "Blue sharks"), Pelagic_thresher_shark_F_Crit$Simulations %>% mutate(Species = "Pelagic thresher sharks")) combined_F_crit_res <- combined_F_crit_sims %>% gather(Stat, Value, -F., - Species) %>% filter(Value >= 1) %>% group_by(Species, Stat) %>% filter( Value == min(Value)) %>% select(-Value) %>% spread(Stat, F.) %>% group_by(Species) %>% mutate_all(funs(if_else(!is.na(.), ., 0))) %>% as.data.frame() # saveRDS(combined_F_crit_sims, "All_F_Crit_1000_sims.RDS") # combined_F_crit_sims <- readRDS("All_F_Crit_1000_sims.RDS") combined_F_crit_res <- combined_F_crit_sims %>% gather(Stat, Value, -F., - Species) %>% filter(Value >= 1) %>% group_by(Species, Stat) %>% filter( Value == min(Value)) %>% select(-Value) %>% spread(Stat, F.) %>% group_by(Species) %>% mutate_all(funs(if_else(!is.na(.), ., 0))) %>% as.data.frame() # Change to Latin names combined_F_crit_sims$Species <-factor(combined_F_crit_sims$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) # Change to Latin names combined_F_crit_res$Species <-factor(combined_F_crit_res$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) F_crit_plot <- ggplot(combined_F_crit_sims, aes( x = F., y = AVG, ymin = lwr, ymax = upr, fill = Species))+ geom_hline(yintercept = 1, col = "red", size = 1)+ geom_rect(data = combined_F_crit_res, inherit.aes = FALSE, aes(xmin = lwr, xmax = upr, ymin = -Inf, ymax = Inf), alpha = .3, fill = "grey")+ geom_vline(data = combined_F_crit_res, aes(xintercept = AVG), size = 1, linetype = "dashed")+ geom_ribbon(alpha = .8)+ geom_line(size = 1, col = "black")+ labs(y = expression(lambda), x = expression(italic("F")))+ scale_fill_viridis_d(guide = F)+ facet_wrap(~Species, ncol = 3, scales = "free") + theme_bw() + theme(strip.text = element_text(face = "italic")) ggsave("F_crit_plot.png",F_crit_plot, width = 7, height = 7, units = "in") ### Age at last capture analysis --------------------------- #NOTE: These can take hours to run for each species when 1000 simulations for each F/AALC combo are used silky_AALC <- Simulate_AALC(1000, silky_data, M.estimators = Silky_M_ests, min.AALC = 20) OCS_AALC <- Simulate_AALC(1000, OCS_data, M.estimators = OCS_M_ests, min.AALC = 20) AML_AALC <- Simulate_AALC(1000, AML_data, min.AALC = 20) silvertip_AALC <- Simulate_AALC(1000, silvertip_data, M.estimators = Silvertip_M_ests, min.AALC = 20) limbatus_AALC <- Simulate_AALC(1000, limbatus_data, M.estimators = limbatus_M_ests, min.AALC = 20) Tiger_shark_AALC <- Simulate_AALC(1000, Tiger_shark_data, M.estimators = Tiger_M_ests, min.AALC = 20) Bigeye_thresher_AALC <- Simulate_AALC(1000, Bigeye_thresher_data, M.estimators = BigEye_M_ests, min.AALC = 20) Mako_AALC <- Simulate_AALC(1000, Mako_data, M.estimators = Mako_M_ests, min.AALC = 20) hammerhead_shark_AALC <- Simulate_AALC(1000, hammerhead_data, min.AALC = 20) Spinner_shark_AALC <- Simulate_AALC(1000, Spinner_shark_data, M.estimators = Spinner_shark_M_ests, min.AALC = 20) Blue_shark_AALC <- Simulate_AALC(1000, Blue_shark_data, min.AALC = 20) Pelagic_thresher_shark_AALC <- Simulate_AALC(1000, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests, min.AALC = 20) AALC_combined_outputs <- bind_rows(as.data.frame(silky_AALC[[1]]) %>% mutate(Species = "Silky sharks"), as.data.frame(OCS_AALC[[1]]) %>% mutate(Species = "Oceanic whitetip sharks"), as.data.frame(AML_AALC[[1]]) %>% mutate(Species = "Grey reef sharks"), as.data.frame(silvertip_AALC[[1]]) %>% mutate(Species = "Silvertip sharks"), as.data.frame(limbatus_AALC[[1]]) %>% mutate(Species = "Common blacktip sharks"), as.data.frame(Tiger_shark_AALC[[1]]) %>% mutate(Species = "Tiger sharks"), as.data.frame(Bigeye_thresher_AALC[[1]]) %>% mutate(Species = "Big eye thesher sharks"), as.data.frame(Mako_AALC[[1]]) %>% mutate(Species = "Shortfin mako sharks"), as.data.frame(hammerhead_shark_AALC[[1]]) %>% mutate(Species = "Scalloped hammerhead sharks"), as.data.frame(Spinner_shark_AALC[[1]]) %>% mutate(Species = "Spinner sharks"), as.data.frame(Blue_shark_AALC[[1]]) %>% mutate(Species = "Blue sharks"), as.data.frame(Pelagic_thresher_shark_AALC[[1]]) %>% mutate(Species = "Pelagic thresher sharks")) # saveRDS(AALC_combined_outputs, "Data/All_AALC_1000_sims_outputs.RDS") # AALC_combined_outputs <- readRDS("Data/All_AALC_1000_sims_outputs.RDS") AALC_combined_results <- bind_rows(as.data.frame(silky_AALC[[2]]) %>% mutate(Species = "Silky sharks"), as.data.frame(OCS_AALC[[2]]) %>% mutate(Species = "Oceanic whitetip sharks"), as.data.frame(AML_AALC[[2]]) %>% mutate(Species = "Grey reef sharks"), as.data.frame(silvertip_AALC[[2]]) %>% mutate(Species = "Silvertip sharks"), as.data.frame(limbatus_AALC[[2]]) %>% mutate(Species = "Common blacktip sharks"), as.data.frame(Tiger_shark_AALC[[2]]) %>% mutate(Species = "Tiger sharks"), as.data.frame(Bigeye_thresher_AALC[[2]]) %>% mutate(Species = "Big eye thesher sharks"), as.data.frame(Mako_AALC[[2]]) %>% mutate(Species = "Shortfin mako sharks"), as.data.frame(hammerhead_shark_AALC[[2]]) %>% mutate(Species = "Scalloped hammerhead sharks"), as.data.frame(Spinner_shark_AALC[[2]]) %>% mutate(Species = "Spinner sharks"), as.data.frame(Blue_shark_AALC[[2]]) %>% mutate(Species = "Blue sharks"), as.data.frame(Pelagic_thresher_shark_AALC[[2]]) %>% mutate(Species = "Pelagic thresher sharks")) # saveRDS(AALC_combined_results, "Data/All_AALC_1000_sims.RDS") # AALC_combined_results <- readRDS("Data/All_AALC_1000_sims.RDS") # Change to Latin names AALC_combined_results$Species <-factor(AALC_combined_results$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) ALL_AALC_plot <- ggplot(AALC_combined_results , aes(x = F., y = AALC, z = mean.lambda))+ facet_wrap(~Species, scales= "free", ncol = 3)+ geom_raster(aes(fill = mean.lambda), interpolate = TRUE, show.legend = T) + stat_contour(aes(col = ..level..), colour = "white",size = 1, breaks= seq(0.6, 1.2,.2))+ stat_contour(aes(col = ..level..), colour = "yellow",size = 1, breaks=1)+ scale_fill_viridis_c(option = "plasma",name = expression(lambda), breaks=seq(0.6, 1.6,.2)) + scale_x_continuous(breaks = seq(0,max(AALC_combined_results$F.),.2), expand = c(0,0))+ scale_y_continuous(breaks = seq(0,max(AALC_combined_results$AALC),2), expand = c(0,0))+ labs(y = "Age-at-last-capture (AALC)", x = expression(italic("F")))+ coord_cartesian(xlim = c(0,1), ylim = c(0,20))+ theme_bw() + theme(plot.title = element_text(hjust = 0.5), strip.text = element_text(face = "italic")) # library(directlabels)# to add text to contours ALL_AALC_plot <- direct.label(ALL_AALC_plot, list("far.from.others.borders", "calc.boxes", "enlarge.box",cex=0.8, fill = "mediumslateblue", rot = 0, "draw.rects", colour = "White")) ALL_AALC_plot ggsave("AALC_plot.png",ALL_AALC_plot, width = 7, height = 7, units = "in") ### Age at first capture analysis --------------------------- #NOTE: These can take hours to run for each species when 1000 simulations for each F/AAFC combo are used silky_AAFC <- Simulate_AAFC(1000, silky_data, M.estimators = Silky_M_ests, max.AAFC = 20) OCS_AAFC <- Simulate_AAFC(1000, OCS_data, M.estimators = OCS_M_ests, max.AAFC = 24) AML_AAFC <- Simulate_AAFC(1000, AML_data, max.AAFC = 20) silvertip_AAFC <- Simulate_AAFC(1000, silvertip_data, M.estimators = Silvertip_M_ests, max.AAFC = 20) limbatus_AAFC <- Simulate_AAFC(1000, limbatus_data, M.estimators = limbatus_M_ests, max.AAFC = 20) Tiger_shark_AAFC <- Simulate_AAFC(1000, Tiger_shark_data, M.estimators = Tiger_M_ests, max.AAFC = 20) Bigeye_thresher_AAFC <- Simulate_AAFC(1000, Bigeye_thresher_data, M.estimators = BigEye_M_ests, max.AAFC = 24) Mako_AAFC <- Simulate_AAFC(1000, Mako_data, M.estimators = Mako_M_ests, max.AAFC = 28) hammerhead_shark_AAFC <- Simulate_AAFC(1000, hammerhead_data, max.AAFC = 20) Spinner_shark_AAFC <- Simulate_AAFC(1000, Spinner_shark_data, M.estimators = Spinner_shark_M_ests, max.AAFC = 20) Blue_shark_AAFC <- Simulate_AAFC(1000, Blue_shark_data, max.AAFC = 16) Pelagic_thresher_shark_AAFC <- Simulate_AAFC(1000, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests, max.AAFC = 24) AAFC_combined_outputs <- bind_rows(as.data.frame(silky_AAFC[[1]]) %>% mutate(Species = "Silky sharks"), as.data.frame(OCS_AAFC[[1]]) %>% mutate(Species = "Oceanic whitetip sharks"), as.data.frame(AML_AAFC[[1]]) %>% mutate(Species = "Grey reef sharks"), as.data.frame(silvertip_AAFC[[1]]) %>% mutate(Species = "Silvertip sharks"), as.data.frame(limbatus_AAFC[[1]]) %>% mutate(Species = "Common blacktip sharks"), as.data.frame(Tiger_shark_AAFC[[1]]) %>% mutate(Species = "Tiger sharks"), as.data.frame(Bigeye_thresher_AAFC[[1]]) %>% mutate(Species = "Big eye thesher sharks"), as.data.frame(Mako_AAFC[[1]]) %>% mutate(Species = "Shortfin mako sharks"), as.data.frame(hammerhead_shark_AAFC[[1]]) %>% mutate(Species = "Scalloped hammerhead sharks"), as.data.frame(Spinner_shark_AAFC[[1]]) %>% mutate(Species = "Spinner sharks"), as.data.frame(Blue_shark_AAFC[[1]]) %>% mutate(Species = "Blue sharks"), as.data.frame(Pelagic_thresher_shark_AAFC[[1]]) %>% mutate(Species = "Pelagic thresher sharks")) # saveRDS(AAFC_combined_outputs, "Data/All_AAFC_1000_sims_outputs.RDS") # AAFC_combined_outputs <- readRDS("Data/All_AAFC_1000_sims_outputs.RDS") AAFC_combined_results <- bind_rows(as.data.frame(silky_AAFC[[2]]) %>% mutate(Species = "Silky sharks"), as.data.frame(OCS_AAFC[[2]]) %>% mutate(Species = "Oceanic whitetip sharks"), as.data.frame(AML_AAFC[[2]]) %>% mutate(Species = "Grey reef sharks"), as.data.frame(silvertip_AAFC[[2]]) %>% mutate(Species = "Silvertip sharks"), as.data.frame(limbatus_AAFC[[2]]) %>% mutate(Species = "Common blacktip sharks"), as.data.frame(Tiger_shark_AAFC[[2]]) %>% mutate(Species = "Tiger sharks"), as.data.frame(Bigeye_thresher_AAFC[[2]]) %>% mutate(Species = "Big eye thesher sharks"), as.data.frame(Mako_AAFC[[2]]) %>% mutate(Species = "Shortfin mako sharks"), as.data.frame(hammerhead_shark_AAFC[[2]]) %>% mutate(Species = "Scalloped hammerhead sharks"), as.data.frame(Spinner_shark_AAFC[[2]]) %>% mutate(Species = "Spinner sharks"), as.data.frame(Blue_shark_AAFC[[2]]) %>% mutate(Species = "Blue sharks"), as.data.frame(Pelagic_thresher_shark_AAFC[[2]]) %>% mutate(Species = "Pelagic thresher sharks")) # saveRDS(AAFC_combined_results, "Data/All_AAFC_1000_sim.RDS") AAFC_combined_results <- readRDS("Data/All_AAFC_1000_sim.RDS") # Change to Latin names AAFC_combined_results$Species <-factor(AAFC_combined_results$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) ALL_AAFC_plot <- ggplot(AAFC_combined_results , aes(x = F., y = AAFC, z = mean.lambda))+ facet_wrap(~Species, scales= "free", ncol = 3)+ geom_raster(aes(fill = mean.lambda), interpolate = TRUE, show.legend = T) + stat_contour(aes(col = ..level..), colour = "white",size = 1, breaks= c(1.3,seq(0.6, 1.2,.2)))+ stat_contour(aes(col = ..level..), colour = "yellow",size = 1, breaks=1)+ scale_fill_viridis_c(option = "plasma",name = expression(lambda), breaks=c(1.3,seq(0.6, 1.2,.2))) + scale_x_continuous(breaks = seq(0,max(AAFC_combined_results$F.),.2),expand = c(0,0))+ scale_y_continuous(breaks = scales::pretty_breaks(8), expand = c(0,0))+ labs(y = "Age-at-first-capture (AAFC)", x = expression(italic("F")))+ theme_bw() + theme(plot.title = element_text(hjust = 0.5), strip.text = element_text(face = "italic")) # library(directlabels)# to add text to contours ALL_AAFC_plot <- direct.label(ALL_AAFC_plot, list("far.from.others.borders", "calc.boxes", "enlarge.box",cex=0.6, fill = "mediumslateblue", rot = 0, "draw.rects", colour = "White")) ALL_AAFC_plot ggsave("AAFC_plot.png",ALL_AAFC_plot, width = 7, height = 7.3, units = "in") ### Harvest Slot Simulations ------------------------------------------------------------- HS_n <- 1000 silky_HS <- Simulate_harvest_slots(HS_n, silky_data, M.estimators = Silky_M_ests, Age.mid.point = 0:28, HS.width = 0:8,max.F = 1) # saveRDS(silky_HS, paste0("Data/Silky_HS_",HS_n,"_sims_outputs.RDS")) OCS_HS <- Simulate_harvest_slots(HS_n, OCS_data, M.estimators = OCS_M_ests, Age.mid.point = 0: 24, HS.width = 0:8,max.F = 1) # saveRDS(OCS_HS, paste0("Data/OCS_HS_",HS_n,"_sims_outputs.RDS")) AML_HS <- Simulate_harvest_slots(HS_n, AML_data, Age.mid.point = 0:20, HS.width = 0:8,max.F = 1) # saveRDS(AML_HS, paste0("Data/AML_HS_",HS_n,"_sims_outputs.RDS")) silvertip_HS <- Simulate_harvest_slots(HS_n, silvertip_data, M.estimators = Silvertip_M_ests, Age.mid.point = 0:20, HS.width = 0:8,max.F = 1) # saveRDS(silvertip_HS, paste0("Data/silvertip_HS_",HS_n,"_sims_outputs.RDS")) limbatus_HS <- Simulate_harvest_slots(HS_n, limbatus_data, M.estimators = limbatus_M_ests, Age.mid.point = 0:20, HS.width = 0:8,max.F = 1) # saveRDS(limbatus_HS, paste0("Data/limbatus_HS_",HS_n,"_sims_outputs.RDS")) Tiger_shark_HS <- Simulate_harvest_slots(HS_n, Tiger_shark_data, M.estimators = Tiger_M_ests,Age.mid.point = 0:33, HS.width = 0:8,max.F = 1) # saveRDS(Tiger_shark_HS, paste0("Data/Tiger_HS_",HS_n,"_sims_outputs.RDS")) Bigeye_thresher_HS <- Simulate_harvest_slots(HS_n, Bigeye_thresher_data, M.estimators = BigEye_M_ests, Age.mid.point = 0:24, HS.width = 0:8,max.F = 1) # saveRDS(Bigeye_thresher_HS, paste0("Data/Bigeye_HS_",HS_n,"_sims_outputs.RDS")) Mako_HS <- Simulate_harvest_slots(HS_n, Mako_data, M.estimators = Mako_M_ests, Age.mid.point = 0:28, HS.width = 0:8,max.F = 1) # saveRDS(Mako_HS, paste0("Data/mako_HS_",HS_n,"_sims_outputs.RDS")) hammerhead_shark_HS <- Simulate_harvest_slots(HS_n, hammerhead_data, Age.mid.point = 0:35, HS.width = 0:8,max.F = 1) # saveRDS(hammerhead_shark_HS, paste0("Data/hammerhead_HS_",HS_n,"_sims_outputs.RDS")) Spinner_shark_HS <- Simulate_harvest_slots(HS_n, Spinner_shark_data, M.estimators = Spinner_shark_M_ests, Age.mid.point = 0:21, HS.width = 0:8,max.F = 1) # saveRDS(Spinner_shark_HS, paste0("Data/Spinner_HS_",HS_n,"_sims_outputs.RDS")) Blue_shark_HS <- Simulate_harvest_slots(HS_n, Blue_shark_data, Age.mid.point = 0:20, HS.width = 0:8,max.F = 1) # saveRDS(Blue_shark_HS, paste0("Data/Blue_HS_",HS_n,"_sims_outputs.RDS")) Pelagic_thresher_shark_HS <- Simulate_harvest_slots(HS_n, Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests, Age.mid.point = 0:24, HS.width = 0:8,max.F = 1) # saveRDS(Pelagic_thresher_shark_HS, paste0("Data/Pelagic_thresher_HS_",HS_n,"_sims_outputs.RDS")) HS_combined_outputs <- bind_rows(silky_HS %>% mutate(Species = "Silky sharks"), OCS_HS %>% mutate(Species = "Oceanic whitetip sharks"), AML_HS %>% mutate(Species = "Grey reef sharks"), silvertip_HS%>% mutate(Species = "Silvertip sharks"), limbatus_HS %>% mutate(Species = "Common blacktip sharks"), Tiger_shark_HS %>% mutate(Species = "Tiger sharks"), Bigeye_thresher_HS %>% mutate(Species = "Big eye thesher sharks"), Mako_HS %>% mutate(Species = "Shortfin mako sharks"), hammerhead_shark_HS %>% mutate(Species = "Scalloped hammerhead sharks"), Spinner_shark_HS %>% mutate(Species = "Spinner sharks"), Blue_shark_HS %>% mutate(Species = "Blue sharks"), Pelagic_thresher_shark_HS %>% mutate(Species = "Pelagic thresher sharks")) # saveRDS(HS_combined_outputs, paste0("Data/All_HS_",HS_n,"_sims_outputs.RDS")) # HS_combined_outputs <- readRDS("Data/All_HS_1000_sims_outputs.RDS") # Change to Latin names HS_combined_outputs$Species <-factor(HS_combined_outputs$Species, levels = c( "Pelagic thresher sharks", "Big eye thesher sharks", "Silvertip sharks", "Grey reef sharks", "Spinner sharks", "Silky sharks", "Common blacktip sharks", "Oceanic whitetip sharks", "Tiger sharks", "Shortfin mako sharks", "Blue sharks", "Scalloped hammerhead sharks" ), labels = c( "Alopias pelagicus", "Alopias superciliosus", "Carcharhinus albimarginatus", "Carcharhinus amblyrhynchos", "Carcharhinus brevipinna", "Carcharhinus falciformis", "Carcharhinus limbatus", "Carcharhinus longimanus", "Galeocerdo cuvier", "Isurus oxyrinchus", "Prionace glauca", "Sphyrna lewini" )) HS_plot <- ggplot(HS_combined_outputs, aes(x = MinAge, y = MaxAge, z = F.))+ facet_wrap(~Species, scales= "free", ncol = 3)+ geom_raster(aes(fill = F.), interpolate = T, show.legend = T) + stat_contour(aes(col = ..level..), colour = "white",size = 1,breaks=c(.1,.3,.5,1,1.2)) + scale_fill_viridis_c(name = expression(italic("F")),breaks=c(.1,.3,.5,1,1.2)) + scale_x_continuous(expand = c(0,0), name = "Minimum legal age (years)", breaks = seq(0,50,2))+ scale_y_continuous(expand = c(0,0), name = "Maximum legal age (years)", breaks = seq(0,50,2))+ theme_bw() + theme(panel.grid = element_blank(), strip.text = element_text(face = "italic")) HS_plot <- direct.label(HS_plot, list("far.from.others.borders", "calc.boxes", "enlarge.box",cex=0.8, fill = "mediumslateblue", rot = 0, "draw.rects", colour = "White")) ggsave("Harvest_slot_plot.png",HS_plot, width = 8.5, height = 8.5, units = "in") ### Estimate mortality distributions ------------------------------------------------------ # create ouput table for MS create_age_invariant_mortality_table <- function(x){ x <- as.data.frame(x[[1]]) y <- x %>% transmute(Method = Method, results = paste0(round(AVG,2)," (",round(low,2)," - ", round(high,2),")")) return(y) } create_age_dependent_mortality_table <- function(x){ max <- length(x) if(max == 1 ) return("No age-dependent methods applied") y <- plyr::ldply(x) %>% mutate(Method = .id) %>% select(-.id) %>% filter(Method != "Age_invariant_mortality") %>% transmute(Method = Method,#id = 1:nrow(.), results = paste0(round(AVG,2)," (",round(low,2)," - ", round(high,2),")")) %>% group_by(Method) %>% mutate(Age = row_number()) %>% spread(Method, results) %>% mutate(Age = 0:(nrow(.)-1)) return(as.data.frame(y)) } Results_table <- bind_rows(create_age_invariant_mortality_table(Estimate_mortality_dists(data = silky_data, M.estimators = Silky_M_ests)) %>% mutate(Species = "Silky sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = OCS_data,M.estimators = OCS_M_ests)) %>% mutate(Species = "Oceanic whitetip sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = silvertip_data, M.estimators = Silvertip_M_ests)) %>% mutate(Species = "Silvertip sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = limbatus_data, M.estimators = limbatus_M_ests)) %>% mutate(Species = "Common blacktip sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Tiger_shark_data, M.estimators = Tiger_M_ests)) %>% mutate(Species = "Tiger sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Bigeye_thresher_data, M.estimators = BigEye_M_ests)) %>% mutate(Species = "Big eye thesher sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Mako_data, M.estimators = Mako_M_ests)) %>% mutate(Species = "Shortfin mako sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = hammerhead_data)) %>% mutate(Species = "Scalloped hammerhead sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Spinner_shark_data, M.estimators = Spinner_shark_M_ests)) %>% mutate(Species = "Spinner sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Blue_shark_data)) %>% mutate(Species = "Blue sharks"), create_age_invariant_mortality_table(Estimate_mortality_dists(data = Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests)) %>% mutate(Species = "Pelagic thresher sharks")) %>% spread(Method,results) silky_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = silky_data, M.estimators = Silky_M_ests)) OCS_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = OCS_data, M.estimators = OCS_M_ests)) silvertip_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = silvertip_data, M.estimators = Silvertip_M_ests)) limbatus_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = limbatus_data, M.estimators = limbatus_M_ests)) Tiger_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Tiger_shark_data, M.estimators = Tiger_M_ests)) Bigeye_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Bigeye_thresher_data, M.estimators = BigEye_M_ests)) Mako_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Mako_data, M.estimators = Mako_M_ests)) hammer_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = hammerhead_data)) spinner_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Spinner_shark_data, M.estimators = Spinner_shark_M_ests)) blue_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Blue_shark_data)) pelagic_mort <- create_age_dependent_mortality_table(Estimate_mortality_dists(data = Pelagic_thresher_data, M.estimators = Pelagic_thresher_M_ests)) xlsx::write.xlsx(file = "Mortality results.xlsx", Results_table,sheetName = "Age_invariant", row.names = FALSE, append = FALSE) xlsx::write.xlsx(file = "Mortality results.xlsx", silky_mort,sheetName = "Silky", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", OCS_mort,sheetName = "OCS", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", silvertip_mort,sheetName = "silvertip", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", limbatus_mort,sheetName = "limbatus", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", Tiger_mort,sheetName = "Tiger", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", Bigeye_mort,sheetName = "Bigeye", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", Mako_mort,sheetName = "Mako", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", hammer_mort,sheetName = "hammerhead", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", spinner_mort,sheetName = "Spinner", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", blue_mort,sheetName = "Blue", row.names = FALSE, append = TRUE) xlsx::write.xlsx(file = "Mortality results.xlsx", pelagic_mort,sheetName = "Pelagic", row.names = FALSE, append = TRUE) ### Export life history values ----------------------------------------------------------- create_life_history_output_table <- function(x){ # growth parameters section growth_pars <- round(as.data.frame(x[["growth"]][2]),2) colnames(growth_pars) <- gsub("pars.","",colnames(growth_pars)) growth.se <- as.data.frame(x[["growth"]][3]) growth_pars[1,]<- paste0(growth_pars[1,]," (", growth.se[1,],")") if(all(!is.na(x$`growth`$corr.matrix))){ growth_pars[1,4] <- paste( paste(as.character(x$`growth`$corr.matrix)[1:3], collapse = " "), paste(as.character(x$`growth`$corr.matrix)[4:6], collapse = " "), paste(as.character(x$`growth`$corr.matrix)[7:9] , collapse = " "), sep = "\n") colnames(growth_pars)[4] <- "Growth parameter correlation matrix" } colnames(growth_pars)[2] <- case_when( x[["growth"]][1] == "logistic" ~ "glog", x[["growth"]][1] == "Von Bertalanffy" ~ "k", x[["growth"]][1] == "Gompertz" ~ "ggom", TRUE ~ "k") # maturity if(x$maturity$`model.type` == "logistic - int/slope"){ mat_pars <- round(as.data.frame(x[["maturity"]][2]),2) colnames(mat_pars) <- gsub("pars.","",colnames(mat_pars)) mat.se <- round(as.data.frame(x[["maturity"]][3]),2) mat_pars[1,]<- paste0(mat_pars[1,]," (", mat.se[1,],")") if(all(!is.na(x$`maturity`$corr.matrix))){ mat_pars[1,3] <- paste( paste(as.character(x$`maturity`$corr.matrix)[1:2], collapse = " "), paste(as.character(x$`maturity`$corr.matrix)[3:4], collapse = " "), sep = "\n") colnames(mat_pars)[3] <- "Maturity parameter correlation matrix" } } else if(x$maturity$`model.type` == "normal"){ mat_pars <- as.data.frame(paste0(round(x$maturity$pars$`mean`,2)," (", round(x$maturity$pars$se,2),")")) colnames(mat_pars) <- "Age-at-maturity" } else { mat_pars <- as.data.frame(paste0(x$maturity$pars$`min`," - ", x$maturity$pars$max)) colnames(mat_pars) <- "Age-at-maturity" } #repro litter.size <- as.data.frame(paste0(x[["litter.size"]][1] ," (",x[["litter.size"]][2],")")) colnames(litter.size) <- "Litter size" Other_pars <- data.frame(Gestation = round(x[["gest.period"]],2), repro = round(x[["repro.cycle"]],2), max.age = x$max.age$`min`) colnames(Other_pars) <- c("Gestation period (years)", "Length of reproductive cycle (years)", "Maximum age (years)") results <- cbind(growth_pars, mat_pars, litter.size, Other_pars) return(results) } life_history_values <- bind_rows(create_life_history_output_table(silky_data) %>% mutate(Species = "Silky sharks"), create_life_history_output_table(OCS_data) %>% mutate(Species = "Oceanic whitetip sharks"), create_life_history_output_table(AML_data) %>% mutate(Species = "Grey reef sharks"), create_life_history_output_table(silvertip_data) %>% mutate(Species = "Silvertip sharks"), create_life_history_output_table(limbatus_data)%>% mutate(Species = "Common blacktip sharks"), create_life_history_output_table(Tiger_shark_data)%>% mutate(Species = "Tiger sharks"), create_life_history_output_table(Bigeye_thresher_data)%>% mutate(Species = "Big eye thesher sharks"), create_life_history_output_table(Mako_data) %>% mutate(Species = "Shortfin mako sharks"), create_life_history_output_table(hammerhead_data) %>% mutate(Species = "Scalloped hammerhead sharks"), create_life_history_output_table(Spinner_shark_data)%>% mutate(Species = "Spinner sharks"), create_life_history_output_table(Blue_shark_data) %>% mutate(Species = "Blue sharks"), create_life_history_output_table(Pelagic_thresher_data)%>% mutate(Species = "Pelagic thresher sharks") ) %>% arrange(Species) %>% select(Species, Linf, k, ggom, glog, L0, t0, `Growth parameter correlation matrix`, intercept, slope, `Age-at-maturity`, `Maturity parameter correlation matrix`, `Maximum age (years)`,`Litter size`, `Gestation period (years)`, `Length of reproductive cycle (years)` ) xlsx::write.xlsx(file = "life_history_values.xlsx", life_history_values, row.names = FALSE)