---
title: "Bycatch_sims_ALE_BBH_CJFAS"
author: "Kerry"
date: "10/27/2022"
output: html_document
---

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

#load libraries
```{r}
library(tidyverse)
library(stringr)
library(rubias)
library(dplyr)
```

**ALEWIFE SIMULATIONS**
#import genotypes for genetic baselines
```{r}
baseline_ALE <- read.csv("../Inputfiles/Alewife_baseline_final.csv", stringsAsFactors = F)
#add own path to inputs
```

#FILTERING DATA
```{r}
baseline_rem_ALE <- baseline_ALE %>%
  group_by(repunit, collection) %>%
  filter(n() > 20) %>%
  ungroup()

#View(baseline_rem_ALE)
```

#simulations rubias
```{r}
#do simultations to assess alewife baseline

arl_ALE <- assess_reference_mc(baseline_rem_ALE, 
                            gen_start_col = 5, 
                            reps = 100, 
                            mixsize = 100)

arl_ALE

rep_ale <- arl_ALE %>%
  group_by(iter, repunit) %>%
  summarize(post_mean = sum(post_mean),
            mle = sum(mle),
            true_rep_prop = sum(omega),
            n = sum(n)) %>%
  group_by(iter) %>%
  mutate(samp_frac = n / sum(n))
```

```{r}
new_labels_ALE <- c(
                    `1GLS` = "GLS",
                    `1NSC` = "NSC",
                    `2BOF` = "BOF",
                    `3NNE` = "NNE",
                    `4MB` = "MB",
                    `5NUN` = "NUN",
                    `6BIS` = "BIS",
                    `7LIS` = "LIS",
                    `8MAT` = "MAT",
                    `9ALB` = "ALB"
                    )

```

```{r}
bc_sim_ref_ALE_rem <- rep_ale %>%   filter(repunit != "1BBH") %>% mutate(repunit = sub("[1-9]_[1-9]","",repunit))

g <- ggplot(bc_sim_ref_ALE_rem, aes(x = true_rep_prop, y = post_mean, fill = repunit)) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_fill_brewer(direction = -1, palette = "Spectral") +
  geom_point(shape = 21, size = 1, stroke = 0.2) +
  xlab("True proportion") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  ylab("Estimated proportion") +
  facet_wrap(repunit~., ncol = 5, nrow = 2, labeller = as_labeller(new_labels_ALE)) +
  geom_abline(intercept = 0, slope = 1) +
  xlim(0, 0.6)+
  ylim(0, 0.6)
g

ggsave("ALE_subrep_sims_finalplot.pdf", width = 10, height = 6)
```

**BLUEBACK HERRING SIMULATIONS**
#import genotypes for genetic baselines
```{r}
baseline_BBH <- read.csv("../Inputfiles/Blueback_baseline_final.csv", stringsAsFactors = F)
#add own path to inputs
```

#FILTERING DATA
```{r}
baseline_rem_BBH <- baseline_BBH %>%
  group_by(repunit, collection) %>%
  filter(n() > 20) %>%
  ungroup()

#View(baseline_rem_BBH)
```

#simulations rubias
```{r}
#do simultations to assess alewife baseline

arl_BBH <- assess_reference_mc(baseline_rem_BBH, 
                            gen_start_col = 5, 
                            reps = 100, 
                            mixsize = 100)

arl_BBH

rep_bbh <- arl_BBH %>%
  group_by(iter, repunit) %>%
  summarize(post_mean = sum(post_mean),
            mle = sum(mle),
            true_rep_prop = sum(omega),
            n = sum(n)) %>%
  group_by(iter) %>%
  mutate(samp_frac = n / sum(n))
```

```{r}
new_labels_BBH <- c(
                    `2_1CAN` = "CAN",
                    `2_2NNE` = "NNE",
                    `3_1MNE` = "MNE",
                    `4_1SNE` = "SNE",
                    `5_1LIS` = "LIS",
                    `5_2MAT` = "MAT",
                    `5_3ALBE` = "ALB",
                    `6_1CF` = "CF",
                    `6_2SAT` = "SAT",
                    `6_3STR` = "STR"
                    )
```

```{r}
bc_sim_ref_BBH_rem <- rep_bbh %>%   filter(repunit != "1_1ALE")

g <- ggplot(bc_sim_ref_BBH_rem, aes(x = true_rep_prop, y = post_mean, fill = repunit)) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  scale_fill_brewer(direction = -1, palette = "Spectral") +
  geom_point(shape = 21, size = 1, stroke = 0.2) +
  xlab("True proportion") +
  theme(legend.title = element_blank()) +
  theme(legend.position = "none") +
  ylab("Estimated proportion") +
  facet_wrap(repunit~., ncol = 5, nrow = 2, labeller = as_labeller(new_labels_BBH)) +
  geom_abline(intercept = 0, slope = 1) +
  xlim(0, 0.6)+
  ylim(0, 0.6)
g

ggsave("BBH_subrep_sims_finalplot.pdf", width = 10, height = 6)
```

**ALEWIFE Reporting groups**
#import ALE bycatch genotypes
```{r}
bycatch_all_subrep_ALE <- read.csv("../Inputfiles/Alewife_bycatch_final.csv", stringsAsFactors = F)
```

#filt groups less than 20
```{r}
bycatch_subrep_all_rem_ALE <- bycatch_all_subrep_ALE %>%
  group_by(repunit, collection) %>%
  filter(n() > 20) %>%
  ungroup()
```

#ML estimates
```{r}
bycatch_all_subrepunit10_mix_est_ml_ALE <- infer_mixture(reference = baseline_rem_ALE, mixture = bycatch_subrep_all_rem_ALE, gen_start_col = 5)
#bycatch_all_subrepunit10_mix_est_ml_ALE
```

#PB corrected estimates
```{r}
bycatch_all_subrepunit10_mix_est_pb_ALE <- infer_mixture(reference = baseline_rem_ALE, mixture = bycatch_subrep_all_rem_ALE, gen_start_col = 5, method = "PB")
#bycatch_all_subrepunit10_mix_est_pb_ALE
```

```{r}
rep_mix_all_10est_sub_pb_ALE <- bycatch_all_subrepunit10_mix_est_pb_ALE$mixing_proportions %>%
  group_by(mixture_collection, repunit) %>%
  summarise(repprop = sum(pi)) %>%
  filter(repunit != "1BBH")

# adding mixing proportions over collections in the repunit
PB_all_ALE <- left_join(rep_mix_all_10est_sub_pb_ALE, bycatch_all_subrepunit10_mix_est_pb_ALE$bootstrapped_proportions)
```

#Plotting PB corrected mixing proportions estimates
```{r}
ALE_plot <- ggplot(PB_all_ALE, aes(x = mixture_collection, y = repprop, fill = repunit))
ALE_plot + geom_col(aes(fill = repunit), colour="black", stat="identity") +
    scale_fill_brewer(palette = "Spectral", direction = -1, labels = c("GLS", "NSC", "BOF", "NNE", "MB", "NUN", "BIS", "LIS", "MAT", "ALB")) +
  scale_x_discrete(breaks=c("615_2015_Feb_Mar","539_611_2015_Dec_Jan","539_611_2014_Feb_Mar", "539_611_2014_Dec_Jan", "539_611_2013_Feb_Mar", "539_611_2013_Dec_Jan", "539_611_2012_Feb_Mar", "539_611_2012_Dec_Jan","537_613_2015_Feb_Mar", "537_613_2013_Feb_Mar", "537_613_2013_Dec_Jan", "521_2015_Feb_Mar", "521_2015_Dec_Jan" ),
                   labels=c("NJLI_HW2_15","SNE_HW1_15", "SNE_HW2_14","SNE_HW1_14", "SNE_HW2_13","SNE_HW1_13", "SNE_HW2_12","SNE_HW1_12", "LIBI_HW2_15", "LIBI_HW2_13", "LIBI_HW1_13","CC_HW2_15", "CC_HW1_15")) +
    theme_bw() +
    scale_y_continuous(breaks=c(0, 0.5, 1)) +
  theme(text = element_text(size=15)) +
  xlab("") +
  ylab("") +
  coord_flip()
  ggsave("ALE_HWS_newlines.pdf")
```

**Blueback herring Reporting groups**
#import ALE bycatch genotypes
```{r}
bycatch_all_subrep_BBH <- read.csv("../Inputfiles/Blueback_bycatch_final.csv", stringsAsFactors = FALSE)
```

#Remove mixture pops with less than 10 individual
```{r}
bycatch_subrep_all_rem_BBH <- bycatch_all_subrep_BBH %>%
  group_by(repunit, collection) %>%
  filter(n() > 10) %>%
  ungroup()
```

#ML_all BBH
```{r}
bycatch_all_subrepunit10_mix_est_ml_BBH <- infer_mixture(reference = baseline_rem_BBH, mixture = bycatch_subrep_all_rem_BBH, gen_start_col = 5)
#bycatch_all_subrepunit10_mix_est_ml_BBH
```

#PB_all BBH
```{r}
bycatch_all_subrepunit10_mix_est_pb_BBH <- infer_mixture(reference = baseline_rem_BBH, mixture = bycatch_subrep_all_rem_BBH, gen_start_col = 5, method = "PB")
#bycatch_all_subrepunit10_mix_est_pb_BBH
```

# for mixing proportions of corrected estimates PB
# filter out BBH from analyses
```{r}
rep_mix_all_10est_sub_pb_BBH <- bycatch_all_subrepunit10_mix_est_pb_BBH$mixing_proportions %>%
  group_by(mixture_collection, repunit) %>%
  summarise(repprop = sum(pi)) %>%
  filter(repunit != "1_1ALE")

# adding mixing proportions over collections in the repunit
PB_all_BBH <- left_join(rep_mix_all_10est_sub_pb_BBH, bycatch_all_subrepunit10_mix_est_pb_BBH$bootstrapped_proportions)
```

#plotting BBH
```{r}
PB_BBH_repmix <- PB_all_BBH %>%
  filter(mixture_collection != "539_611_2014_Feb_Mar" & mixture_collection != "539_611_2014_Dec_Jan" & mixture_collection != "539_611_2013_Feb_Mar" ) 

p <- ggplot(PB_BBH_repmix, aes(x = mixture_collection, y = repprop, fill = repunit)) 
p + geom_col(aes(fill = repunit)) +
    scale_fill_brewer(palette = "Spectral",  direction = -1, labels = c("CAN", "NNE", "MNE", "SNE", "LIS", "MAT", "ALB", "CF", "SAT", "STR")) +
    geom_col(aes(fill = repunit), colour="black", stat="identity") +
   scale_x_discrete(breaks=c("615_2015_Feb_Mar","539_611_2015_Dec_Jan", "539_611_2013_Dec_Jan", "539_611_2012_Feb_Mar", "539_611_2012_Dec_Jan","537_613_2015_Feb_Mar", "537_613_2012_Dec_Jan", "521_2015_Dec_Jan", "521_2014_May", "521_2013_Dec_Jan", "521_2012_Feb_Mar"), labels=c("NJLI_HW2_15","SNE_HW1_15", "SNE_HW1_13", "SNE_HW2_12","SNE_HW1_12", "LIBI_HW2_15", "LIBI_HW1_12","CC_HW2_15", "CC_SP1_14", "CC_HW1_13", "CC_HW2_12")) +
    theme_bw() +
  scale_y_continuous(breaks=c(0, 0.5, 1)) +
  theme(text = element_text(size=15)) +
  #theme(legend.title = element_blank(), legend.position = "none") +
  xlab("") +
  ylab("") +
  coord_flip()
  ggsave("BBH_HWs_NEW_lines.pdf")
```
