Index_Recruit.Rmd

---
title: "Index recruit Antarctic krill"
subtitle: "Suplementary information HeatWaves Antarctic Paper"
author:
  - Mardones Mauricio^[Instituo de Fomento Pesquero, mardones.mauricio@gmail.com]
  - Lucas Kruger^[Instituto Antártico Chileno]
  - Lorena Rebolledo^[Instituto Antártico Chileno]
date:  "`r format(Sys.time(), '%d %B, %Y')`"
bibliography: seaice.bib
csl: apa.csl
link-citations: yes
linkcolor: blue
always_allow_html: yes
output:
  pdf_document
---

```{=tex}
\fontsize{12}{16}
\selectfont{}
```
\newpage

```{=latex}
\setcounter{tocdepth}{3}
\tableofcontents
```
\newpage

```{r setup1}
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE,
                      message = FALSE,
                      warning = FALSE,
                      fig.align = 'center',
                      dev = 'jpeg',
                      dpi = 300,
                      out.width='120%')
#XQuartz is a mess, put this in your onload to default to cairo instead
options(bitmapType = "cairo") 
# (https://github.com/tidyverse/ggplot2/issues/2655)
# Lo mapas se hacen mas rapido
```


```{r lib, message=F, echo= TRUE}
library(here)
#statistics
library(ggsignif)
library(lubridate)
library(easystats) # multiples unciones analiticas
library(readxl)
# vizualizacion
library(ggrepel)
library(ggpubr)
library(ggridges)
library(sf)
library(GGally)
library(tidyverse, quietly = TRUE)
library(knitr, quietly = TRUE)
library(kableExtra)
library(raster)
library(egg)
library(car) #Variance inflation Factor
library(ggthemes)
library(sjPlot)
library(CCAMLRGIS)
```

# Contexto

Este análisis tiene como objetivo calcular un índice de reclutamiento (IR) del krill *Euphausia superba* a través de los datos monitoreados por la flota pesquera entre los años 1980 y 2020

# Metodología

Los pasos para calcular el indice son los siguentes:

- Calcular el cuantíl 90% de la REF como proporción de individuos que ingresan a la pesquería. 

- Con este valor del 90%, se calculan los individuos ajo esta talla provenientes de la pesquería. Luego se calcula la cantidad de inviduos bajo esa talla agrupados por distintas covariables.

- Se calcula la proporción respecto al total. Se identifica la distribución de los datos, luego se normalizan para que los datos tengan una distribución normal, y posteriormente,  se estandarizan para llevarlos a un indice e tre -1 y 1, lo cual indica reclutamientos negativos y positivos respectivamente.

- Despliegue en diferentes plot y agrupaciones temporales y espaciales.

## Load data

## Data exploratory analysis

The object `ohbio2` come from data exploration analysis in data request
CCAMLR data. This objetc have bio information from krill.

```{r  echo=FALSE}
#datos entregados por secretaria de CCMLAR
metadata <- load("~/DOCAS/Data/565_C1_KRI_2021-10-01/DATA_PRODUCT_565.RData")
# Data procesada por MMardones
#ohbio <- load("DataLengthKrill.RData")
#ohbio
#metadata
#son lo mismo
```

```{r}
#cargo objeto
meta <- get("METADATA")
c1 <- get("C1")
ohbio <- get("OBS_HAUL_BIOLOGY")
names(ohbio)
dim(c1)
dim(ohbio)
```


Join data set with master as `c1` set. This join is trought
`obs_haul_id` variable to get geoposition variables

```{r warning=FALSE}
ohbio2 <- left_join(c1, ohbio, by="obs_haul_id")
dim(ohbio2)
```

Firsts glance. Test how many register have by year. In this case,
`length_total_cm` by season ccamlr. Same exercise in date period
`date_catchperiod_start` to separate dates.

```{r}
ohbio3 <- ohbio2 %>%
  mutate(Year = year(date_catchperiod_start),
         Month = month(date_catchperiod_start),
         Day = day(date_catchperiod_start))
```

filter necesary data to further analysis
```{r}
length481 <-ohbio3 %>% 
  dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47, 43) %>% 
  filter(asd_code=="481")
#save(length481, file = "length481.RData")
```

```{r}
ohbio4 <- ohbio3 %>% 
dplyr::select(7, 9, 11, 12, 14, 24, 25, 29, 42, 44, 46, 47)
names(ohbio4)

```

## Maps works

First thing is get different rater layer to join krill data length
according different porpoises.

```{r raster}

# Cargo linea de costa
coast <- load_Coastline()
coast1<- st_as_sf(coast) 
coast2 = st_transform(coast1, "+proj=latlong +ellps=WGS84")

# con Statistical Areas con foco en 48.1
suba <- load_ASDs()
suba1 <- subset(suba[(3),])
suba1a<- st_as_sf(suba1) 
suba1aa = st_transform(suba1a, "+proj=latlong +ellps=WGS84")

# Uso las agrupaciones de Strata
strata <- st_read("~/DOCAS/Mapas/Antarctic_SHPfiles/Strata.shp",
                quiet=T)
strata=st_transform(strata, "+proj=latlong +ellps=WGS84")
strata <- strata %>% 
  dplyr::filter(ID != "Outer")

```

## Strata maps

Show strata agregation to join length data (Figure\@ref(fig:maptest).

```{r maptest, fig.cap="Strata Maps in 48.1"}
# y testeo el mapa
ssmap <- ggplot(strata)+
  geom_sf(data = strata, aes(fill=strata$ID, 
                           alpha=0.5))+
  geom_sf(data = suba1aa,alpha=0.3,
          colour="red")+
  geom_sf(data = coast2, colour="black", fill=NA)+
  scale_fill_viridis_d(option = "F",
                       name="Stratum")+
  ylim(230000, 2220000)+
  xlim(-3095349 , -1858911)+
  coord_sf(crs = 6932)+
  scale_alpha(guide="none")+
  theme_bw()
ssmap
```

## Grouping Length data into Strata

```{r ssmu1}
names(ohbio4)
ohbio5 <- ohbio4 %>% 
  drop_na() %>% 
  filter(asd_code==481) %>% 
  dplyr::select(6, 7, 8, 9, 10, 11, 12) 
ohbio6 <- st_as_sf(ohbio5 %>% 
                     drop_na(latitude_set_end), 
                   coords = c("longitude_set_end", 
                              "latitude_set_end"),  
                  crs = "+proj=latlong +ellps=WGS84")
```


Comprobar si tengo datos duplicados

```{r}
# comoprobar si tengo datos duplicados
strata2 <- st_make_valid(strata)
ohbio7 <- st_make_valid(ohbio6)
krill.strata <- st_join(strata2, ohbio7)
#saveRDS(krill.strata, "KrillData.Rdata")
colSums(is.na(krill.strata))
dim(krill.strata)

```

## Exploración primaria

```{r message=FALSE, warning=FALSE}
jz3 <- ggplot(krill.strata %>% 
                drop_na(), 
              aes(x=length_total_cm, y = as.factor(Month), 
                      fill = factor(stat(quantile))))+
  stat_density_ridges(
    geom = "density_ridges_gradient",
    calc_ecdf = TRUE,
     quantiles = c(0.10, 0.90)) +
  scale_fill_manual(
    name = "Probability", 
    values = c("#de2d26", "#fee0d2", "#de2d26"),
    labels = c("[0 - 0.10]", 
               "[0.10 - 0.90]", 
               "[0.90 - 1]"))+
  facet_grid(ID~Year) +   
  geom_vline(xintercept = 3.6, color = "red")+
  scale_x_continuous(breaks = seq(from = 3, 
                                  to = 12, 
                                  by = 2))+
  scale_y_discrete(breaks = seq(from = 1, 
                                to = 12, 
                                by = 3))+
  #scale_fill_viridis_d(name="SubArea")+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        legend.position = "bottom")+
  #xlim(10,120)+
  xlab("Total Length (cm.)")+
  ylab("")
jz3
```
Ahora identifico los distintos cuantiles de los datos de pesquería y estaciones


```{r warning=FALSE}
# Calcular cuantiles por grupo
cuantiles_por_strata <- tapply(krill.strata$length_total_cm, 
                              krill.strata$ID,
                              function(x) quantile(x, 
                                                   c(0.10,
                                                     0.5, 
                                                     0.90)))

```

Calculo el índice del reclutamiento de [@Maschette2020]

```{r warning=FALSE}
indice_reclutamiento <- krill.strata %>%
  filter(length_total_cm< 3.6 ) %>% 
  group_by(Year, Month, ID) %>%
  summarize(PROP = n() / nrow(krill.strata)) %>% 
  mutate(PROPLOG =log(PROP))
# Crear gráficos en facet_wrap de barras para representar el índice de reclutamiento
```

ahora estandarizo los datos entre -1 y 1.

```{r warning=FALSE}
a <- -1  # Límite inferior del rango objetivo
b <- 1   # Límite superior del rango objetivo

# Calcular el valor mínimo y máximo de tus datos
min_x <- min(indice_reclutamiento$PROPLOG)
max_x <- max(indice_reclutamiento$PROPLOG)

# Aplicar la fórmula de normalización
indice_reclutamiento$PROPLOG2 <- ((indice_reclutamiento$PROPLOG- min_x) / (max_x - min_x)) * (b - a) + a
```


veo la distribucion de las variables

```{r warning=FALSE, message=FALSE}
nor <- ggplot(indice_reclutamiento, aes(PROP)) + 
  geom_histogram(fill="grey")+
  theme_few()
log <- ggplot(indice_reclutamiento, aes(PROPLOG)) + 
  geom_histogram(fill="grey")+
  theme_few()

ggarrange(nor, log, ncol = 2)
```


```{r warning=FALSE}
indrec1 <- ggplot(indice_reclutamiento , 
       aes(x = factor(Year), 
           y = PROP,
           fill=ID)) +
  geom_boxplot() +
  facet_wrap(ID~., ncol=5) +
  scale_fill_viridis_d(option = "F",
                       name="Stratum")+
  scale_x_discrete(breaks = seq(from = 1996, 
                                to = 2022, 
                                by = 4))+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  labs(x = "ANO", 
        y = "Índice de Reclutamiento")+
  ylim(0, 0.001)
indrec1
```
ahora como columnas por meses pero cambio los nombres

```{r}
# Crear una nueva columna con los nombres de los meses
indice_reclutamiento2 <- indice_reclutamiento %>%
  mutate(Month = case_when(
    Month == 1 ~ "January",
    Month == 2 ~ "February",
    Month == 3 ~ "March",
    Month == 4 ~ "April",
    Month == 5 ~ "May",
    Month == 6 ~ "June",
    Month == 7 ~ "July",
    Month == 8 ~ "August",
    Month == 9 ~ "September",
    Month == 10 ~ "October",
    Month == 11 ~ "November",
    Month == 12 ~ "December")) %>% 
  mutate(Month = factor(Month, levels = c("January", 
                                          "February", 
                                          "March",
                                          "April",
                                          "May",
                                          "June",
                                          "July",
                                          "August",
                                          "September",
                                          "October",
                                          "November", 
                                          "December")))


#ahora trimestres

# Definir una función para asignar trimestres
# Definir reglas para asignar trimestres
indice_reclutamiento2<- indice_reclutamiento2 %>%
  mutate(quarter = case_when(
    Month %in% c("January", "February", "March") ~ "Q1",
    Month %in% c("April", "May", "June") ~ "Q2",
    Month %in% c("July", "August", "September") ~ "Q3",
    Month %in% c("October", "November", "December") ~ "Q4"))

```



```{r warning=FALSE, message=FALSE}
indseg3 <- ggplot(indice_reclutamiento2  %>% 
  group_by(Year,ID) %>%
  summarise(PROPLOG3=mean(PROPLOG2)), 
       aes(x = factor(Year), 
           y = PROPLOG3,
           fill=PROPLOG3 > 0)) +
  geom_bar(stat = "identity")  +
  scale_fill_manual(values = c("black", "red"),
                    labels = c("Negativo", "Positivo"),
                    name="Index Recruit Krill") +
  facet_wrap(.~ID, ncol = 5) +
  geom_hline(yintercept = 0, color = "red")+
  scale_x_discrete(breaks = seq(from = 1996, to = 2022, by = 1))+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text = element_text(size=6),
        legend.position = "bottom")+
  labs(x = "",
       y ="")+
  coord_flip()
indseg3
```


```{r warning=FALSE, message=FALSE}
indrec4 <- ggplot(indice_reclutamiento2  %>% 
  group_by(Year,ID,Month) %>%
  summarise(PROPLOG3=mean(PROPLOG2)), 
       aes(x = Month, 
           y = PROPLOG3,
           fill=PROPLOG3 > 0)) +
  geom_bar(stat = "identity")  +
  scale_fill_manual(values = c("black", "red"),
                    labels = c("Negative", "Positive"),
                    name="Index Recruit Krill") +
  facet_grid(ID~Year) +
  geom_hline(yintercept = 0, color = "red")+
  #scale_x_discrete(breaks = seq(from = 1, to = 12, by = 2))+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text = element_text(size=6),
        legend.position = "bottom")+
  labs(x = "", 
        y = "")+
  coord_flip()
indrec4
```

Por trimestres

```{r warning=FALSE, message=FALSE}
indrec5 <- ggplot(indice_reclutamiento2  %>% 
  group_by(Year,ID,quarter) %>%
  summarise(PROPLOG3=mean(PROPLOG2)), 
       aes(x = quarter, 
           y = PROPLOG3,
           fill=PROPLOG3 > 0)) +
  geom_bar(stat = "identity")  +
  scale_fill_manual(values = c("black", "red"),
                    labels = c("Negative", "Positive"),
                    name="Index Recruit Krill") +
  facet_grid(ID~Year) +
  geom_hline(yintercept = 0, color = "red")+
  #scale_x_discrete(breaks = seq(from = 1, to = 12, by = 2))+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text = element_text(size=6),
        legend.position = "bottom")+
  labs(x = "", 
        y = "")+
  coord_flip()
indrec5
```

Grafico como Oscilación por Strata

```{r warning=FALSE, message=FALSE}
recosc <- ggplot(indice_reclutamiento2  %>% 
  group_by(Year,ID) %>%
  summarise(PROPLOG3=mean(PROPLOG2)), 
       aes(x = Year, y = PROPLOG3)) +
  geom_ribbon(aes(ymin = pmax(PROPLOG3, 0), 
                  ymax = 0), 
              fill = "#de2d26", 
              alpha = 0.8) +  # Área por encima de la línea
  geom_ribbon(aes(ymin = pmin(PROPLOG3, 0), 
                  ymax = 0), 
              fill = "black", 
              alpha = 0.8) +  # Área por debajo de la línea
  geom_hline(yintercept = 0, color = "red")+
  geom_line(color = "black") +  # Línea de anomalías
  geom_point( alpha=0.2,
              size= 0.9)+
  labs(x = "AÑO", y = "Index Krill recruit") +
   facet_wrap(.~ID, ncol = 5)+
  theme_few()+
  theme(axis.text.x = element_text(angle = 90, hjust = 2))
recosc

```
# References