# en caso no se tengan las dependencias
# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
# install.packages("caTools")
# library("caTools")
# # install.packages("r4ss")
library(r4ss)
library(ss3diags)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags) # diagnosticos del modelo
library(ss3sim) # evaluación de sesgo
library(tidyverse)
library(kableExtra)
library(egg)
library(reshape2)
library(ggridges)Evaluación de stock Almeja (Venus antiqua; Ameghinomya antiqua)
Análisis de sensibilidad y desempeño del los modelos como insumo para Programa de Seguimiento de Pesquerías bajo Planes de Manejo. Año 2023
1. ANTECEDENTES
1.1. Descripción y objetivo del documento
Este documento contiene el flujo de análisis y modelación de los distintos escenarios de Venus antiqua almeja para los stocl de Ancud y la zona sur de la Regiñon de Los Lagos como parte de la asesoría técnica que lleva a cabo el IFOP, en el marco del programa de seguimiento de pesquerías bajo planes de Manejo.
1.2. Modelo de evaluación
El modelo de dinámica poblacional de la almeja, corresponde a un enfoque de evaluación del tipo estadístico con estructura de edad, donde la dinámica progresa avanzando en el tiempo t, y las capturas son causantes de la mortalidad por pesca F, la mortalidad natural es constante M = 0,28. La relación entre la población y las capturas responde a la base de la ecuación de Baranov, y se consideran para el modelo y estimaciones el rango de edad entre 2 a 12+ (años). Sin embargo, las estimaciones del modelo tienen su origen en la edad cero sobre la base de una condición inicial estado estable. La dinámica esta modelada por un reclutamiento tipo Beverton y Holt.
A continuación se presentan los escenarios que mas influencia reconocida tienen sobre la dinámica poblacional de almeja, entre ellos el nivel de escarpamiento de la relación stock-recluta (stepness) y la motralidad natural M. Para ello se configuran 7 escenarios de distintas combinaciones y evaluar el impacto en las distintas variables poblacionales como ; biomasas, F y reclutamiento.
En ambos procesos de evaluación (unidades poblacionales de Ancud y X Sur) se consideraron los siguientes escenarios:
| Escenario | Descripción |
|---|---|
| s1 | Modelo base (stepness = 0.8; M= 0.28) |
| s2 | Stepness = 0.8, M = 0.2 |
| s3 | Stepness = 0.7 |
| s4 | Stepness = 0.9 |
| s5 | M = 0.3 |
| s6 | Stepness = 0.7 y M = 0.3 |
| s7 | Stepness = 0.7 y M = 0.4 |
1.3. Plataforma de modelación
Los modelos implementados fueron configurados utilizando Stock Synthesis (SS3 de aqui en mas)(https://vlab.noaa.gov/web/stock-synthesis), que es un modelo de evaluación de stock edad y talla estrucuturado, en la clase de modelo denominado “Modelo de análisis integrado”. SS tiene un sub-modelo poblacional de stock que simula crecimiento, madurez, fecundidad, reclutamiento, movimiento, y procesos de mortalidad, y sub-modelos de observation y valores esperados para diferentes tipos de datos. El modelo es codificado en C++ con parámetros de estimación activados por diferenciación automática (ADMB) (Methot & Wetzel, 2013). El análisis de resultados y salidas emplea herramientas de R e interfase gráfica de la librería r4ss (https://github.com/r4ss/r4ss) (Taylor, 2019).
Se realiza una modelación con la plataforma SS3 (V.3.30.19) y sus outputs leídos con la librería “r4ss” (Taylor, 2019)
2. Lectura y Run
model ## 2.1. Librerías
Instalación de las librerías necesarias
2.1. Run models
Los componentes de verosimilitud, además de los análisis de residuales permiten identificar entre los bloques de modelos cuales de las configuraciones presenta un desempeño adecuado en términos estadísticos de ajuste a la información.
Estos modelos, son los seleccionados para presentar en sus principales salidas para fines informativos de indicadores, puntos biológicos de referencia y estimaciones poblacionales.
Identifico las carpetas;
# Defino directorios Para Ancud
dirmain <- here("SA_2023", "Ancud")
dir1<-here("SA_2023", "Ancud", "s1") # Modelo todas las flotas
dir2<-here("SA_2023", "Ancud", "s2")
dir3<-here("SA_2023", "Ancud", "s3")
dir4<-here("SA_2023", "Ancud", "s4")
dir5<-here("SA_2023", "Ancud", "s5")
dir6<-here("SA_2023", "Ancud", "s6")
dir7<-here("SA_2023", "Ancud", "s7")
# Defino Directorios para X Sur
dirmains <- here("SA_2023", "XSur")
dir1s<-here("SA_2023", "XSur", "s1")
dir2s<-here("SA_2023", "XSur", "s2")
dir3s<-here("SA_2023", "XSur", "s3")
dir4s<-here("SA_2023", "XSur", "s4")
dir5s<-here("SA_2023", "XSur", "s5")
dir6s<-here("SA_2023", "XSur", "s6")
dir7s<-here("SA_2023", "XSur", "s7")Ahora corro todos los modelos juntos de Ancud y X Sur
dirancud <- c(
file.path(dir1),
file.path(dir2),
file.path(dir3),
file.path(dir4),
file.path(dir5),
file.path(dir6),
file.path(dir7)
)
# Corro todos los modelos de Ancud
r4ss::run_SS_models(dirvec = dirancud, model = './ss_osx',
skipfinished = FALSE)
# Considerar que al correr el modelo dentro de cada `chunk` se cambia el directorio.
# Para X Sur
dirxsur <- c(
file.path(dir1s),
file.path(dir2s),
file.path(dir3s),
file.path(dir4s),
file.path(dir5s),
file.path(dir6s),
file.path(dir7s)
)
# Corro todos los modelos de XSur
r4ss::run_SS_models(dirvec = dirxsur, model = './ss_osx',
skipfinished = FALSE)2.3. Leyendo salidas
Leo las salidas de todos los escenarios (7 por stock; Ancud y X Sur) de los modelos con la función SS_output.
base.model1 <- SS_output(dir=dir1,covar=F,forecast=F)
base.model2 <- SS_output(dir=dir2,covar=F,forecast=F)
base.model3 <- SS_output(dir=dir3,covar=F,forecast=F)
base.model4 <- SS_output(dir=dir4,covar=F,forecast=F)
base.model5 <- SS_output(dir=dir5,covar=F,forecast=F)
base.model6 <- SS_output(dir=dir6,covar=F,forecast=F)
base.model7 <- SS_output(dir=dir7,covar=F,forecast=F)
#Leo modelos de XSur
base.model1s <- SS_output(dir=dir1s,covar=F,forecast=F)
base.model2s <- SS_output(dir=dir2s,covar=F,forecast=F)
base.model3s <- SS_output(dir=dir3s,covar=F,forecast=F)
base.model4s <- SS_output(dir=dir4s,covar=F,forecast=F)
base.model5s <- SS_output(dir=dir5s,covar=F,forecast=F)
base.model6s <- SS_output(dir=dir6s,covar=F,forecast=F)
base.model7s <- SS_output(dir=dir7s,covar=F,forecast=F)2.4. Datos Disponibles
Los datos y tablas de parametrós son obtenidas del modelo base de cada stock, es decir s1.
2.4.1. Ancud
Data disponible para este escenario. se identifica como flota buzodado que es la principal fuente de información de la flota. El rendimiento es medido como;
CPUE = Captura /(HoraBuceo * Numero de Buzos)
donde las capturas son en kilogramos.
SSplotData(base.model1, subplot = 2,
fleetnames = "Buzo",
fleetcol = 4)
Leo los datos de los templates de SS3 de Ancud modelo base s1
Tallas composiciones Ancud
Debo extrear el dato del .daty trasnformarlo a un objeto geom
# Utilizo solo el modelo base "s1"
leancud <- as.data.frame(dat$lencomp[,c(1,7:44)])
# Aplicar la función melt()
df_long <- melt(leancud, id.vars = "Yr",
variable.name = "Variable",
value.name = "Valor")
df_long1 <- gsub("l", "", df_long$Variable)
df_long$Variable <- df_long1
# Desagrupo
long2 <- df_long %>%
drop_na(Variable) %>%
type.convert(as.is = TRUE) %>%
uncount(Valor)Un grafico simple
nbco <- ggplot(long2, aes(x=Variable, y = as.factor(Yr)))+
geom_density_ridges(stat = "binline", bins = 39,
scale = 3,
draw_baseline = FALSE,
alpha=.5,
fill=9)+
geom_vline(xintercept = 55, color = "red")+
scale_y_discrete(breaks = seq(from = 1996, to = 2022, by = 2))+
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none")+
xlab("Longitud (mm.)")+
ylab("")
nbco
CPUE
cpueancud <- as.data.frame(dat$CPUE[,c(1,4)])cpuea <- ggplot(cpueancud,
aes(x = year, y = obs)) +
geom_point(shape = 16, size = 2, alpha = 0.6) +
stat_smooth(method = "loess")+
theme(panel.background = element_blank(),
panel.border = element_rect(fill = NA, size = 1),
strip.background = element_rect(fill = "white",
color = "white", size = 1)) +
theme(text = element_text(size=13)) +
theme_bw()+
xlab("") +
ylab("CPUE (Captura/N Buzos*HOras de Buceo)")+
ylim(10,100)Capturas Ancud
CapAncud <- dat$catch
kbl(CapAncud, booktabs = T,format = "html",
caption = "Estadísticas de Pesca almeja Ancud") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) %>%
landscape()| year | seas | fleet | catch | catch_se |
|---|---|---|---|---|
| -999 | 3 | 1 | 0.00 | 0.05 |
| 1965 | 3 | 1 | 105.00 | 0.05 |
| 1966 | 3 | 1 | 229.00 | 0.05 |
| 1967 | 3 | 1 | 600.00 | 0.05 |
| 1968 | 3 | 1 | 255.00 | 0.05 |
| 1969 | 3 | 1 | 2383.00 | 0.05 |
| 1970 | 3 | 1 | 792.00 | 0.05 |
| 1971 | 3 | 1 | 824.00 | 0.05 |
| 1972 | 3 | 1 | 2686.00 | 0.05 |
| 1973 | 3 | 1 | 855.00 | 0.05 |
| 1974 | 3 | 1 | 67.00 | 0.05 |
| 1975 | 3 | 1 | 300.00 | 0.05 |
| 1976 | 3 | 1 | 403.00 | 0.05 |
| 1977 | 3 | 1 | 2590.00 | 0.05 |
| 1978 | 3 | 1 | 8105.00 | 0.05 |
| 1979 | 3 | 1 | 11591.00 | 0.05 |
| 1980 | 3 | 1 | 10572.00 | 0.05 |
| 1981 | 3 | 1 | 10023.00 | 0.05 |
| 1982 | 3 | 1 | 9661.00 | 0.05 |
| 1983 | 3 | 1 | 10950.00 | 0.05 |
| 1984 | 3 | 1 | 14137.00 | 0.05 |
| 1985 | 3 | 1 | 17751.00 | 0.05 |
| 1986 | 3 | 1 | 20191.00 | 0.05 |
| 1987 | 3 | 1 | 17980.00 | 0.05 |
| 1988 | 3 | 1 | 18644.00 | 0.05 |
| 1989 | 3 | 1 | 16806.00 | 0.05 |
| 1990 | 3 | 1 | 9797.00 | 0.05 |
| 1991 | 3 | 1 | 11041.00 | 0.05 |
| 1992 | 3 | 1 | 7390.00 | 0.05 |
| 1993 | 3 | 1 | 3429.00 | 0.05 |
| 1994 | 3 | 1 | 1474.00 | 0.05 |
| 1995 | 3 | 1 | 2019.00 | 0.05 |
| 1996 | 3 | 1 | 4359.00 | 0.05 |
| 1997 | 3 | 1 | 1975.76 | 0.05 |
| 1998 | 3 | 1 | 7604.84 | 0.05 |
| 1999 | 3 | 1 | 1621.77 | 0.05 |
| 2000 | 3 | 1 | 1159.00 | 0.05 |
| 2001 | 3 | 1 | 5971.00 | 0.05 |
| 2002 | 3 | 1 | 2054.00 | 0.05 |
| 2003 | 3 | 1 | 2043.00 | 0.05 |
| 2004 | 3 | 1 | 7936.00 | 0.05 |
| 2005 | 3 | 1 | 980.00 | 0.05 |
| 2006 | 3 | 1 | 4466.00 | 0.05 |
| 2007 | 3 | 1 | 2160.00 | 0.05 |
| 2008 | 3 | 1 | 788.00 | 0.05 |
| 2009 | 3 | 1 | 993.00 | 0.05 |
| 2010 | 3 | 1 | 1152.00 | 0.05 |
| 2011 | 3 | 1 | 981.00 | 0.05 |
| 2012 | 3 | 1 | 3022.00 | 0.05 |
| 2013 | 3 | 1 | 114.00 | 0.05 |
| 2014 | 3 | 1 | 208.00 | 0.05 |
| 2015 | 3 | 1 | 657.00 | 0.05 |
| 2016 | 3 | 1 | 855.00 | 0.05 |
| 2017 | 3 | 1 | 870.00 | 0.05 |
| 2018 | 3 | 1 | 1496.00 | 0.05 |
| 2019 | 3 | 1 | 1496.00 | 0.05 |
| 2020 | 3 | 1 | 1352.00 | 0.05 |
| 2021 | 3 | 1 | 2566.00 | 0.05 |
| 2022 | 3 | 1 | 2426.00 | 0.05 |
2.4.2. X Sur
SSplotData(base.model1s, subplot = 2,
fleetnames = "Buzo",
fleetcol = 3)
Leo los datos de los templates de SS3 de X Sur modelo base s1
Tallas composiciones ZSur
Debo extrear el dato del .daty trasnformarlo a un objeto geom
# Utilizo solo el modelo base "s1"
lexsur <- as.data.frame(dats$lencomp[,c(1,7:44)])
# Aplicar la función melt()
df_longxsur <- melt(lexsur, id.vars = "Yr",
variable.name = "Variable",
value.name = "Valor")
df_longxsur1 <- gsub("l", "", df_longxsur$Variable)
df_longxsur$Variable <- df_longxsur1
# Desagrupo
longsur2 <- df_longxsur %>%
drop_na(Variable) %>%
type.convert(as.is = TRUE) %>%
uncount(Valor)Un grafico simple
nbsur <- ggplot(longsur2, aes(x=Variable, y = as.factor(Yr)))+
geom_density_ridges(stat = "binline", bins = 38,
scale = 3,
draw_baseline = FALSE,
alpha=.5,
fill="blue")+
geom_vline(xintercept = 55, color = "red")+
scale_y_discrete(breaks = seq(from = 1996, to = 2022, by = 2))+
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none", legend.title = element_blank())+
xlab("Longitud (mm.)")+
ylab("")
nbsur
CPUE
cpuexsur <- as.data.frame(dats$CPUE[,c(1,4)])cpuexs <- ggplot(cpuexsur,
aes(x = year, y = obs)) +
geom_point(shape = 16, size = 2, alpha = 0.6) +
stat_smooth(method = "loess",
color="red")+
theme(panel.background = element_blank(),
panel.border = element_rect(fill = NA, size = 1),
strip.background = element_rect(fill = "white",
color = "white", size = 1)) +
theme(text = element_text(size=13)) +
theme_bw()+
xlab("") +
ylab("CPUE (Captura/N Buzos*HOras de Buceo)")+
ylim(10,100)Ambas zonas
#Plot as a grid
ggarrange(cpuea, cpuexs, nrow = 1)
Capturas X Sur
CapXSur <- dats$catch
kbl(CapXSur, booktabs = T,format = "html",
caption = "Estadísticas de Pesca almeja X Sur") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) %>%
landscape()| year | seas | fleet | catch | catch_se |
|---|---|---|---|---|
| -999 | 3 | 1 | 0 | 0.05 |
| 1967 | 3 | 1 | 59 | 0.05 |
| 1968 | 3 | 1 | 256 | 0.05 |
| 1969 | 3 | 1 | 343 | 0.05 |
| 1970 | 3 | 1 | 194 | 0.05 |
| 1971 | 3 | 1 | 267 | 0.05 |
| 1972 | 3 | 1 | 530 | 0.05 |
| 1973 | 3 | 1 | 302 | 0.05 |
| 1974 | 3 | 1 | 341 | 0.05 |
| 1975 | 3 | 1 | 176 | 0.05 |
| 1976 | 3 | 1 | 303 | 0.05 |
| 1977 | 3 | 1 | 413 | 0.05 |
| 1978 | 3 | 1 | 3000 | 0.05 |
| 1979 | 3 | 1 | 7617 | 0.05 |
| 1980 | 3 | 1 | 1491 | 0.05 |
| 1981 | 3 | 1 | 1762 | 0.05 |
| 1982 | 3 | 1 | 2766 | 0.05 |
| 1983 | 3 | 1 | 4069 | 0.05 |
| 1984 | 3 | 1 | 4698 | 0.05 |
| 1985 | 3 | 1 | 3295 | 0.05 |
| 1986 | 3 | 1 | 7094 | 0.05 |
| 1987 | 3 | 1 | 4939 | 0.05 |
| 1988 | 3 | 1 | 11826 | 0.05 |
| 1989 | 3 | 1 | 9980 | 0.05 |
| 1990 | 3 | 1 | 8762 | 0.05 |
| 1991 | 3 | 1 | 19350 | 0.05 |
| 1992 | 3 | 1 | 18502 | 0.05 |
| 1993 | 3 | 1 | 11294 | 0.05 |
| 1994 | 3 | 1 | 10805 | 0.05 |
| 1995 | 3 | 1 | 8440 | 0.05 |
| 1996 | 3 | 1 | 10771 | 0.05 |
| 1997 | 3 | 1 | 5655 | 0.05 |
| 1998 | 3 | 1 | 14040 | 0.05 |
| 1999 | 3 | 1 | 8095 | 0.05 |
| 2000 | 3 | 1 | 10652 | 0.05 |
| 2001 | 3 | 1 | 18054 | 0.05 |
| 2002 | 3 | 1 | 2802 | 0.05 |
| 2003 | 3 | 1 | 5725 | 0.05 |
| 2004 | 3 | 1 | 13937 | 0.05 |
| 2005 | 3 | 1 | 11116 | 0.05 |
| 2006 | 3 | 1 | 11322 | 0.05 |
| 2007 | 3 | 1 | 7856 | 0.05 |
| 2008 | 3 | 1 | 13706 | 0.05 |
| 2009 | 3 | 1 | 12758 | 0.05 |
| 2010 | 3 | 1 | 18325 | 0.05 |
| 2011 | 3 | 1 | 17375 | 0.05 |
| 2012 | 3 | 1 | 9202 | 0.05 |
| 2013 | 3 | 1 | 4883 | 0.05 |
| 2014 | 3 | 1 | 7282 | 0.05 |
| 2015 | 3 | 1 | 8515 | 0.05 |
| 2016 | 3 | 1 | 8776 | 0.05 |
| 2017 | 3 | 1 | 7525 | 0.05 |
| 2018 | 3 | 1 | 5111 | 0.05 |
| 2019 | 3 | 1 | 5111 | 0.05 |
| 2020 | 3 | 1 | 4578 | 0.05 |
| 2021 | 3 | 1 | 5278 | 0.05 |
| 2022 | 3 | 1 | 4400 | 0.05 |
2.4.3. Parámetros de historia de vida
Vizualizo los parámetros biológicos utilizados en el modelo base. Cabe señalar que se utilizan los mismos parámetros biológicos para ambos stock.
parbio <-round(ctl$MG_parms,2)
kbl(parbio, booktabs = T,format = "html",
caption = "Parámetros biológicos") %>%
kable_styling(latex_options = c("striped",
"condensed","scale_down"),
full_width = FALSE) %>%
landscape()| LO | HI | INIT | PRIOR | PR_SD | PR_type | PHASE | env_var&link | dev_link | dev_minyr | dev_maxyr | dev_PH | Block | Block_Fxn | PType | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NatM_p_1_Fem_GP_1 | 0.10 | 0.80 | 0.28 | 0.30 | 99.0 | 3 | -2 | 0 | 0 | 0 | 0 | 0.5 | 0 | 0 | 1 |
| L_at_Amin_Fem_GP_1 | 20.00 | 60.00 | 25.00 | 25.00 | 99.0 | 0 | -2 | 0 | 0 | 0 | 0 | 0.5 | 0 | 0 | 2 |
| L_at_Amax_Fem_GP_1 | 20.00 | 150.00 | 96.45 | 98.00 | 110.0 | 0 | -2 | 0 | 0 | 0 | 0 | 0.5 | 0 | 0 | 2 |
| VonBert_K_Fem_GP_1 | 0.01 | 1.00 | 0.13 | 0.15 | 99.0 | 0 | -3 | 0 | 0 | 0 | 0 | 0.5 | 0 | 0 | 2 |
| CV_young_Fem_GP_1 | 0.03 | 0.15 | 0.07 | 0.15 | 99.0 | 0 | 5 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 2 |
| CV_old_Fem_GP_1 | 0.03 | 0.15 | 0.07 | 0.15 | 99.0 | 0 | 5 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 2 |
| Wtlen_1_Fem_GP_1 | -3.00 | 3.00 | 0.01 | 0.01 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 3 |
| Wtlen_2_Fem_GP_1 | -3.00 | 3.00 | 3.06 | 3.06 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 3 |
| Mat50%_Fem_GP_1 | -3.00 | 70.00 | 43.50 | 43.50 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 4 |
| Mat_slope_Fem_GP_1 | -3.00 | 3.00 | -0.15 | -0.15 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 4 |
| Eggs/kg_inter_Fem_GP_1 | -3.00 | 3.00 | 1.00 | 1.00 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 5 |
| Eggs/kg_slope_wt_Fem_GP_1 | -3.00 | 3.00 | 0.00 | 0.00 | 99.0 | 0 | -50 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 5 |
| CohortGrowDev | 0.10 | 10.00 | 1.00 | 1.00 | 1.0 | 0 | -1 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 11 |
| FracFemale_GP_1 | 0.00 | 1.00 | 0.50 | 0.50 | 0.5 | 0 | -99 | 0 | 0 | 0 | 0 | 0.0 | 0 | 0 | 14 |
3. Resultados
3.1. Ancud
SSplotCohortCatch(base.model1, subplots = 1)
AJuste de tallas por flota
SSplotComps(base.model1, subplots = 1,
maxrows = 6,
maxcols = 5,
linescol = 4)
Otros plots
SSplotDynamicB0(base.model1, uncertainty = F)
#SSplotSPR(base.model3)SSplotPars(base.model1)
SSplotTimeseries(base.model1, subplot = 7)
SSplotSPR(base.model1, subplots = 4)
SSplotSummaryF(base.model1)
SSplotTimeseries(base.model1, subplot = 9)
3.2. x Sur
SSplotCohortCatch(base.model1s, subplots = 1)
AJuste de tallas por flota
SSplotComps(base.model1s, subplots = 1,
maxrows = 6,
maxcols = 5,
linescol = 3)
Otros plots
SSplotDynamicB0(base.model1s, uncertainty = F)
#SSplotSPR(base.model3)
SSplotPars(base.model1s)
SSplotTimeseries(base.model1s, subplot = 7)
SSplotSPR(base.model1s, subplots = 4)
SSplotSummaryF(base.model1s)
SSplotTimeseries(base.model1s, subplot = 9)
4. Diagnóstico
4.1. Ancud
4.1.1. Comparación entre escenarios
#PLOT labels, name of each model run
legend.labels <- c('s1', 's2' , 's3','s4', 's5', 's6', 's7')
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglist <- SSgetoutput(keyvec = NULL,
dirvec = c(dir1, dir2, dir3, dir4, dir5, dir6, dir7),
getcovar = F)#create summary of model runs from list above
summaryoutput <- SSsummarize(biglist)col=c('#b2182b','#ef8a62','#fddbc7','#ffffff','#e0e0e0','#999999','#4d4d4d')
SSplotComparisons(summaryoutput,
subplots= c(3,7),
legendlabels = legend.labels,
plot = TRUE,
col = col)Tablas comparativas para Ancud
tablecom <- SStableComparisons(summaryoutput,
likenames = c("TOTAL", "Survey",
"Length_comp",
"Age_comp", "priors",
"Size_at_age"),
names = c("Recr_Virgin", "R0",
"steep", "NatM",
"L_at_Amax", "VonBert_K",
"SSB_Virg", "Bratio_2022",
"SPRratio_2022"),
digits = NULL,
modelnames = c('s1', 's2' , 's3','s4',
's5', 's6', 's7'),
csv = FALSE,
csvdir = dirmain,
csvfile = "parameter_comparison_table.csv",
verbose = TRUE,
mcmc = FALSE)kbl(tablecom,
longtable = F,
booktabs = T,
caption = "Outputs por escenario Ancud") %>%
kable_styling(latex_options = c("scale_down",
"condensed",
"hold_position"))| Label | s1 | s2 | s3 | s4 | s5 | s6 | s7 |
|---|---|---|---|---|---|---|---|
| TOTAL_like | 117.662000 | 224.566000 | 105.9690000 | 256.6120000 | 256.6120000 | 142.3520000 | 142.2920000 |
| Survey_like | -2.031080 | 20.598500 | -2.5733800 | 29.7023000 | 29.7023000 | 15.7690000 | 12.0602000 |
| Length_comp_like | 109.034000 | 128.479000 | 107.1580000 | 155.3840000 | 155.3840000 | 141.7280000 | 135.3750000 |
| Parm_priors_like | 0.000000 | 0.000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 | 0.0000000 |
| Recr_Virgin_thousands | 20.085500 | 20.085500 | 20.1059000 | 20.0855000 | 20.0855000 | 23.2868000 | 36.1705000 |
| SR_LN(R0) | 3.000000 | 3.000000 | 3.0010100 | 3.0000000 | 3.0000000 | 3.1478900 | 3.5882400 |
| SR_BH_steep | 0.800000 | 0.800000 | 0.7000000 | 0.9000000 | 0.9000000 | 0.7000000 | 0.7000000 |
| NatM_p_1_Fem_GP_1 | 0.280000 | 0.200000 | 0.2800000 | 0.2800000 | 0.2800000 | 0.3000000 | 0.4000000 |
| L_at_Amax_Fem_GP_1 | 96.450000 | 96.450000 | 96.4500000 | 96.4500000 | 96.4500000 | 96.4500000 | 96.4500000 |
| VonBert_K_Fem_GP_1 | 0.130000 | 0.130000 | 0.1300000 | 0.1300000 | 0.1300000 | 0.1300000 | 0.1300000 |
| SSB_Virgin_thousand_mt | 20.493000 | 58.002000 | 20.5960000 | 43.8000000 | 43.8000000 | 40.6120000 | 22.0430000 |
| Bratio_2022 | 0.116453 | 0.100168 | 0.0920532 | 0.0630878 | 0.0630878 | 0.0190572 | 0.0428671 |
| SPRratio_2022 | 0.834316 | 0.855079 | 0.8494910 | 0.9158110 | 0.9158110 | 0.9451810 | 0.9059800 |
4.1.2. Retrospectivo
#do retrospective model runs
SS_doRetro(dir1, '', newsubdir = "retrospectives",
subdirstart = "retro",
years = 0:-5,
overwrite = TRUE, exefile = "ss",
extras = "-nox -nohess",
intern = FALSE, CallType = "system",
RemoveBlocks = FALSE)retroModels <- SSgetoutput(dirvec=file.path(dir1,
"retrospectives",
paste("retro",0:-5,sep="")))
retroSummary <- SSsummarize(retroModels, verbose = F)
endyrvec <- retroSummary$endyrs + 0:-5SSplotComparisons(retroSummary,
endyrvec=endyrvec,
legendlabels=paste("Data",0:-5,"years"),
plotdir=dirmain,
plot=TRUE,print=T)
TableCompare <- SStableComparisons(retroSummary,
names=names,
modelnames=c('B','-1','-2','-3','-4','-5'),
csv=TRUE,
csvfile="RetroRunsAncud.csv",verbose=FALSE)
writing table to:
/Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/RetroRunsAncud.csvkbl(TableCompare,
longtable = F,
booktabs = T,
caption = "Comparación de parámetros en retrospectivo Ancud") %>%
kable_styling(latex_options = c("scale_down",
"condensed",
"hold_position"))| Label | B | -1 | -2 | -3 | -4 | -5 |
|---|---|---|---|---|---|---|
| TOTAL_like | 124.288000 | 111.23600 | 96.92950 | 88.80960 | 78.55280 | 67.10840 |
| Survey_like | 0.685869 | -2.20614 | -5.93946 | -7.59118 | -6.69037 | -6.37544 |
| Length_comp_like | 112.033000 | 101.57000 | 89.40540 | 83.22400 | 75.38980 | 68.50600 |
| Parm_priors_like | 0.000000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 |
Otro tipo de diagnostico con ss3diag
SSplotRunstest(base.model1,
add = T,
legendcex = 0.8,
subplot = "cpue",
verbose = F)
Index runs.p test sigma3.lo sigma3.hi type
1 Fishery 0.003 Failed -0.478438 0.478438 cpueParámetro de Rho
sspar(mfrow = c(1, 2), plot.cex = 0.8)
rb = SSplotRetro(retroSummary,
add = T, forecast = F,
legend = F, verbose = F)
rf = SSplotRetro(retroSummary, add = T, subplots = "F",
ylim = c(0, 1), forecast = F,
legendloc = "topleft",
legendcex = 0.8, verbose = F)
4.2. X Sur
4.2.1. Comparación entre escenarios
#PLOT labels, name of each model run
legend.labels <- c('s1', 's2' , 's3','s4', 's5', 's6', 's7')
#read in all model runs
#note if cover=T you need a hessian; if covar=F you do not need a hessian
biglistx <- SSgetoutput(keyvec = NULL,
dirvec = c(dir1s, dir2s, dir3s, dir4s,
dir5s, dir6s, dir7s),
getcovar = F)#create summary of model runs from list above
summaryoutputx <- SSsummarize(biglistx)SSplotComparisons(summaryoutputx,
subplots= c(3,7),
legendlabels = legend.labels,
col=col)Tablas comparativas para X Sur
tablecomx <- SStableComparisons(summaryoutputx,
likenames = c("TOTAL", "Survey",
"Length_comp",
"Age_comp", "priors",
"Size_at_age"),
names = c("Recr_Virgin", "R0",
"steep", "NatM",
"L_at_Amax", "VonBert_K",
"SSB_Virg", "Bratio_2022",
"SPRratio_2022"),
digits = NULL,
modelnames = c('s1', 's2' , 's3','s4',
's5', 's6', 's7'),
csv = FALSE,
csvdir = dirmains,
csvfile = "parameter_comparison_table.csv",
verbose = TRUE,
mcmc = FALSE)kbl(tablecomx,
longtable = F,
booktabs = T,
caption = "Outputs por escenario X Sur") %>%
kable_styling(latex_options = c("scale_down",
"condensed",
"hold_position"))| Label | s1 | s2 | s3 | s4 | s5 | s6 | s7 |
|---|---|---|---|---|---|---|---|
| TOTAL_like | 17.347400 | 18.990100 | 17.734000 | 17.16510 | 17.216200 | 17.491400 | 16.987100 |
| Survey_like | -19.922900 | -18.328900 | -20.091500 | -19.80100 | -20.333700 | -20.486200 | -22.457100 |
| Length_comp_like | 61.835800 | 60.239400 | 61.949300 | 61.76850 | 62.256200 | 62.354000 | 64.432700 |
| Parm_priors_like | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 |
| Recr_Virgin_thousands | 44.634300 | 21.643700 | 46.466400 | 43.27820 | 53.337300 | 55.244800 | 132.201000 |
| SR_LN(R0) | 3.798500 | 3.074720 | 3.838730 | 3.76765 | 3.976640 | 4.011770 | 4.884330 |
| SR_BH_steep | 0.800000 | 0.800000 | 0.700000 | 0.90000 | 0.800000 | 0.700000 | 0.700000 |
| NatM_p_1_Fem_GP_1 | 0.280000 | 0.200000 | 0.280000 | 0.28000 | 0.300000 | 0.300000 | 0.400000 |
| L_at_Amax_Fem_GP_1 | 96.450000 | 96.450000 | 96.450000 | 96.45000 | 96.450000 | 96.450000 | 96.450000 |
| VonBert_K_Fem_GP_1 | 0.130000 | 0.130000 | 0.130000 | 0.13000 | 0.130000 | 0.130000 | 0.130000 |
| SSB_Virgin_thousand_mt | 45.116000 | 61.502000 | 46.959000 | 43.74800 | 42.499000 | 44.014000 | 35.023000 |
| Bratio_2022 | 0.396936 | 0.246470 | 0.374195 | 0.41595 | 0.438471 | 0.416022 | 0.642663 |
| SPRratio_2022 | 0.585145 | 0.717544 | 0.587563 | 0.58288 | 0.552999 | 0.555374 | 0.405268 |
4.2.2. Retrospectivo
#do retrospective model runs
SS_doRetro(dir1s, '', newsubdir = "retrospectives",
subdirstart = "retro",
years = 0:-5,
overwrite = TRUE, exefile = "ss",
extras = "-nox -nohess", intern = FALSE, CallType = "system",
RemoveBlocks = FALSE)retroModelsx <- SSgetoutput(dirvec=file.path(dir1s,
"retrospectives",
paste("retro",0:-5,sep="")))
length(dirvec) as input to SSgetoutput: 6
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro0/Report.sso
added element 'replist1' to list
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro-1/Report.sso
added element 'replist2' to list
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro-2/Report.sso
added element 'replist3' to list
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro-3/Report.sso
added element 'replist4' to list
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro-4/Report.sso
added element 'replist5' to list
reading output from /Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/SA_2023/XSur/s1/retrospectives/retro-5/Report.sso
added element 'replist6' to list
retroSummaryx <- SSsummarize(retroModelsx)
endyrvecx <- retroSummaryx$endyrs + 0:-5SSplotComparisons(retroSummaryx,
endyrvec=endyrvecx,
legendlabels=paste("Data",0:-5,"years"),
plotdir=dirmains,
plot=TRUE,print=T)
Parámetro de Rho
sspar(mfrow = c(1, 2), plot.cex = 0.8)
rb = SSplotRetro(retroSummaryx,
add = T, forecast = F,
legend = F, verbose = F)
rf = SSplotRetro(retroSummaryx, add = T, subplots = "F",
ylim = c(0, 1), forecast = F,
legendloc = "topleft",
legendcex = 0.8, verbose = F)
TableComparex <- SStableComparisons(retroSummaryx,
names=names,
modelnames=c('B','-1','-2','-3','-4','-5'),
csv=TRUE,
csvfile="RetroRunsXsur.csv",verbose=FALSE)
writing table to:
/Users/mauriciomardones/IFOP/Prog_Seg_PM_Bentónicos_2023/Almeja_X_XI/RetroRunsXsur.csvkbl(TableComparex,
longtable = F,
booktabs = T,
caption = "Comparación de parámetros en retrospectivo Ancud") %>%
kable_styling(latex_options = c("scale_down",
"condensed",
"hold_position"))| Label | B | -1 | -2 | -3 | -4 | -5 |
|---|---|---|---|---|---|---|
| TOTAL_like | 17.3474 | 13.7264 | 8.06341 | 8.29356 | 8.30107 | 8.40243 |
| Survey_like | -19.9229 | -19.7624 | -18.90180 | -17.56710 | -16.43520 | -15.25970 |
| Length_comp_like | 61.8358 | 57.9637 | 52.60980 | 51.67460 | 50.73320 | 49.84500 |
| Parm_priors_like | 0.0000 | 0.0000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 |
5. REFERENCIAS
Methot, R. D., & Wetzel, C. R. (2013). Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research, 142, 86–99. https://doi.org/10.1016/j.fishres.2012.10.012
Taylor, I. (2019). Using R for Stock Synthesis Installing R and getting R4SS. Fisheries Science, November.