Modelos avanzados en evaluación de recursos pesqueros: Día 1

Dr. Giancarlo M. Correa

Cousteau Consultant Group

Modelos estructurados por edad

Sistema pesquero

¿Para qué nos sirven los modelos de evaluación?

Conceptos básicos

Stockholm University Baltic Sea Centre

Conceptos básicos

Diferencias claves en enfoques para estudiar poblaciones:

Haddon (2011)

Enfoque integrado

Combina diferentes tipos de datos en un único análisis. El ajuste a los datos se realiza por medio de una rutina de minimización de alguna función objetivo (D. Fournier and Archibald 1982).

Tipos de datos más comunes:

• Capturas
• Índices de abundancia
• Composición por tallas o edades
• Marcaje y recaptura

Enfoque integrado

Haddon (2011)

Conceptos básicos

En un modelo de evaluación estructurado por edades tenemos los siguientes componentes principalmente:

• Reclutamiento
• Mortalidad (natural y por pesca)
• Selectividad
• Crecimiento somático

Plataforma de evaluación

Software o librería orientado a la implementación de modelos de evaluación de stocks (en general).

Ventajas

• No hay necesidad de codificar
• Comunicación más sencilla

Desventajas

• No pueden ser aplicables a todos los stocks
• Efecto de caja negra

Plataforma de evaluación

Ejemplos:

• Stock Synthesis (Methot and Wetzel 2013)
• MULTIFAN-CL (Daid A. Fournier, Hampton, and Sibert 1998)
• SAM (Nielsen and Berg 2014)
• WHAM (Stock and Miller 2021)

Introducción a modelos state-space

Modelos state-space

Definición

• SSM permiten diferenciar y separar variaciones en procesos de un sistema de los errores de observación.

Podemos tener dos series temporales:

1. Serie de procesos, o estados, que representan los cambios verdaderos en el sistema.
2. Serie de observaciones que realizamos del sistema.

Los procesos cambian en el tiempo estocásticamente, lo que se conoce como variaciones en procesos (process variation).

Modelos state-space

• \(z_{t}\): estado en tiempo \(t\) (e.g., abundancia de la población)
• \(y_{t}\): observaciones

Auger-Méthé et al. (2021)

Modelos state-space

Normalmente se asume un proceso Markov (de primer orden) para modelar los estados:

\[z_t = \beta z_{t-1}+\epsilon_t\]

• \(\beta\): autocorrelación de estados en el tiempo
• \(\epsilon\): variación de procesos, donde:

\[\epsilon_t \sim N(0,\sigma_p^2)\]

\(z_0\) puede ser tratado como un parámetro adicional.

Modelos state-space

Para el caso de las observaciones, se asumen que estas son independientes, y se relacionan con los estados:

\[y_t = \alpha z_{t}+\eta_t\]

• \(\alpha\): parámetro que regula el sesgo entre observación y estado
• \(\eta\): error de observación, donde:

\[\eta_t \sim N(0,\sigma_o^2)\]

Modelos state-space

Podemos caracterizar este modelo en términos de distribuciones de probabilidad:

\[g(y_t \mid z_{t}, \boldsymbol{\theta}_o)\]

\[f(z_t \mid z_{t-1}, \boldsymbol{\theta}_p)\]

para \(t=1,...,T\), donde \(f\) y \(g\) son dos funciones de densidad de probabilidad normales, \(\boldsymbol{\theta}_p\) y \(\boldsymbol{\theta}_o\) son vectores de parámetros asociados con cada ecuación:

\[\boldsymbol{\theta}_p = (\beta, \sigma_p^2), \boldsymbol{\theta}_o = (\alpha,\sigma_o^2)\]

Modelos state-space

Los efectos fijos vienen a ser:

\[\boldsymbol{\theta} = (\boldsymbol{\theta}_p, \boldsymbol{\theta}_o, z_0)\]

Y los estados \(z_t\) son tratados como variables aleatorias.

Por lo tanto, los SSM se consideran estos modelos jerárquicos:

1. Primer nivel: ecuación de observación
2. Segundo nivel: ecuación de proceso

Modelos state-space

• Aplicaciones a una gran variedad de campos de investigación.

• Enfoque utilizado para analizar series de tiempo en ecología:
  • movimiento de animales
  • biologging
  • datos de captura-recaptura
  • epidemiología
  • evaluación de stocks

SSM en evaluación de stocks

Enfoque presentado por Sullivan (1992) y Gudmundsson (1994).

Usualmente los estados \(z_t\) vienen a ser variables como (Aeberhard, Flemming, and Nielsen 2018):

• Abundancia a la edad
• Reclutamiento
• Mortalidad por pesca
• Mortalidad natural
• Selectividad

SSM en evaluación de stocks

En modelos de evaluación tradicionales se emplea el enfoque de máxima verosimilitud penalizada (PML) para modelar variaciones temporales en algún estado.

Ejemplo

Para modelar mortalidad por pesca, se separa la mortalidad por pesca anual y la selectividad a la edad:

\[F_{a,t} = f_t S_a \exp(\delta_{a,t})\]

Asumiendo que \(\delta_{a,t} \sim N(0,\sigma_\delta^2)\). El enfoque PML fija \(\sigma_\delta^2\) (penalty term) en el modelo a un valor subjetivo.

Principal ventaja de SSM: estimación de \(\sigma_\delta^2\).

SSM en evaluación de stocks

En la actualidad, TMB (Kristensen et al. 2016) es la principal plataforma para implementar SSM debido a su eficiencia para modelar efectos aleatorios.

• Desarrollado a partir de ADMB (David A. Fournier et al. 2012)
• Usa la maquinaria de C++
• Sencilla manipulación desde R

Laboratorio

Software requerido

• MS Excel
• R instalar
• RStudio instalar
• Rtools instalar

Instalación de librerías

Instalar TMB:

install.packages('TMB')

Para instalar WHAM, utilizaremos la versión que permite incluir datos de tallas:

remotes::install_github(repo = 'GiancarloMCorrea/wham', ref='growth', 
                        INSTALL_opts = c("--no-docs", "--no-multiarch", "--no-demo"))

También les sugiero explorar la versión base de WHAM.

Modelo sencillo en TMB

Objetivo: Conocer los pasos para implementar un modelo en TMB

Ejemplo obtenido de Sean Anderson (Gompertz SSM).

Modelo sencillo en TMB

1. Tener en claro el sistema a modelar:

\[u_t = a + b u_{t-1} + \epsilon_t\] \[y_t = ln(u_t)\]

\[\epsilon_t \sim N(0,\sigma_p^2)\] \[y_t \sim N(ln(u_t),\sigma_o^2)\]

Donde \(u_t\) es el estado (abundancia) en la población.

Modelo sencillo en TMB

2. Implementar tu modelo en TMB (ver Ejercicio_1.cpp)

// Header
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {

// Input data:
DATA_VECTOR(y);

// Parameters (fixed and random variables):
PARAMETER(a); 
PARAMETER(b); 
PARAMETER(log_sigma_proc); 
PARAMETER(log_sigma_obs); 
PARAMETER_VECTOR(u); 

// Transformed parameters if required
Type sigma_proc = exp(log_sigma_proc);
Type sigma_obs = exp(log_sigma_obs);

// Reports on transformed parameters:
ADREPORT(sigma_proc)
ADREPORT(sigma_obs)

int n = y.size(); // get time series length

Type nll = 0.0; // initialize negative log likelihood

// process model:
for(int i = 1; i < n; i++){
  Type m = a + b * u[i - 1]; // Gompertz
  nll -= dnorm(u[i], m, sigma_proc, true); // nll function
}

// observation model:
for(int i = 0; i < n; i++){
  nll -= dnorm(y[i], u[i], sigma_obs, true); // nll function
}

return nll;
}

Modelo sencillo en TMB

2. Implementar tu modelo en TMB (ver Ejercicio_1.cpp)

// Header
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {

// Input data:
DATA_VECTOR(y);

// Parameters (fixed and random variables):
PARAMETER(a); 
PARAMETER(b); 
PARAMETER(log_sigma_proc); 
PARAMETER(log_sigma_obs); 
PARAMETER_VECTOR(u); 

// Transformed parameters if required
Type sigma_proc = exp(log_sigma_proc);
Type sigma_obs = exp(log_sigma_obs);

// Reports on transformed parameters:
ADREPORT(sigma_proc)
ADREPORT(sigma_obs)

int n = y.size(); // get time series length

Type nll = 0.0; // initialize negative log likelihood

// process model:
for(int i = 1; i < n; i++){
  Type m = a + b * u[i - 1]; // Gompertz
  nll -= dnorm(u[i], m, sigma_proc, true); // nll function
}

// observation model:
for(int i = 0; i < n; i++){
  nll -= dnorm(y[i], u[i], sigma_obs, true); // nll function
}

return nll;
}

Modelo sencillo en TMB

2. Implementar tu modelo en TMB (ver Ejercicio_1.cpp)

// Header
#include <TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() () {

// Input data:
DATA_VECTOR(y);

// Parameters (fixed and random variables):
PARAMETER(a); 
PARAMETER(b); 
PARAMETER(log_sigma_proc); 
PARAMETER(log_sigma_obs); 
PARAMETER_VECTOR(u); 

// Transformed parameters if required
Type sigma_proc = exp(log_sigma_proc);
Type sigma_obs = exp(log_sigma_obs);

// Reports on transformed parameters:
ADREPORT(sigma_proc)
ADREPORT(sigma_obs)

int n = y.size(); // get time series length

Type nll = 0.0; // initialize negative log likelihood

// process model:
for(int i = 1; i < n; i++){
  Type m = a + b * u[i - 1]; // Gompertz
  nll -= dnorm(u[i], m, sigma_proc, true); // nll function
}

// observation model:
for(int i = 0; i < n; i++){
  nll -= dnorm(y[i], u[i], sigma_obs, true); // nll function
}

return nll;
}

Modelo sencillo en TMB

3. Compilar el modelo (en R):

require(TMB)
compile("Ejercicio_1.cpp")
dyn.load(dynlib("Ejercicio_1"))

Modelo sencillo en TMB

4. En R, leer tus datos. Imaginemos que y son tus observaciones:

print(y)

  [1] 3.857919 4.139282 4.166160 4.593082 4.487591 4.691392 4.876311 4.681899
  [9] 4.581784 4.706909 4.361975 4.900340 4.493509 4.841002 4.922498 4.722094
 [17] 5.041822 4.885910 4.346571 4.497004 4.620252 4.214646 4.341192 4.251093
 [25] 4.634612 4.130690 4.092495 4.414942 4.317848 4.315265 4.970341 4.852596
 [33] 4.682839 4.766740 4.560845 5.270630 4.776654 5.206920 5.310690 4.500237
 [41] 4.816574 4.481971 4.218075 4.016740 4.537255 4.935705 4.412348 4.750295
 [49] 4.955589 4.473484 4.852497 4.890937 4.776721 4.486947 4.931865 4.767796
 [57] 5.192605 4.571766 4.964695 4.688365 4.988027 4.610332 4.421622 5.253123
 [65] 4.336354 4.339095 4.583634 4.512310 4.740320 4.904109 5.148162 4.948719
 [73] 4.386315 5.102028 4.383329 4.215316 4.716833 4.701634 4.490838 4.427061
 [81] 4.322670 4.828525 4.610247 4.428890 4.703014 4.641656 4.965105 4.956500
 [89] 5.095534 4.696261 5.029934 5.024695 5.091339 4.819369 4.511168 5.412903
 [97] 4.909512 4.939771 5.217307 4.853442

Modelo sencillo en TMB

5. Preparar los objetos y correr el modelo:

# Data input (list)
data <- list(y = y) 
# fixed effects and random variables (list):
parameters <- list(a = 1, b = 0.5, log_sigma_proc = -1,
                   log_sigma_obs = -1, u = rep(mean(y), 100))
# Construct objective function (define random variable):
obj <- MakeADFun(data, parameters, random = "u", DLL = "Ejercicio_1")
obj$hessian <- FALSE # no calculate hessian
# Optimization:
opt <- nlminb(start = obj$par, obj = obj$fn, gr = obj$gr)
# Calculate SD of model parameters
est_pars <- sdreport(obj)

Modelo básico de evaluación de stocks

Abrir archivo Ejercicio_2.xlsx

Referencias

Aeberhard, William H., Joanna Mills Flemming, and Anders Nielsen. 2018. “Review of State-Space Models for Fisheries Science.” Annual Review of Statistics and Its Application 5 (1): 215–35. https://doi.org/10.1146/annurev-statistics-031017-100427.

Auger-Méthé, Marie, Ken Newman, Diana Cole, Fanny Empacher, Rowenna Gryba, Aaron A. King, Vianey Leos-Barajas, et al. 2021. “A Guide to Statespace Modeling of Ecological Time Series.” Ecological Monographs 91 (4). https://doi.org/10.1002/ecm.1470.

Fournier, Daid A, John Hampton, and John R Sibert. 1998. “MULTIFAN-CL: A Length-Based, Age-Structured Model for Fisheries Stock Assessment, with Application to South Pacific Albacore, Thunnus Alalunga.” Canadian Journal of Fisheries and Aquatic Sciences 55 (9): 2105–16.

Fournier, David A., Hans J. Skaug, Johnoel Ancheta, James Ianelli, Arni Magnusson, Mark N. Maunder, Anders Nielsen, and John Sibert. 2012. “AD Model Builder: Using Automatic Differentiation for Statistical Inference of Highly Parameterized Complex Nonlinear Models.” Optimization Methods and Software 27 (2): 233–49. https://doi.org/10.1080/10556788.2011.597854.

Fournier, David, and Chris P. Archibald. 1982. “A General Theory for Analyzing Catch at Age Data.” Canadian Journal of Fisheries and Aquatic Sciences 39 (8): 1195–1207. https://doi.org/10.1139/f82-157.

Gudmundsson, Gudmundur. 1994. “Time Series Analysis of Catch-at-Age Observations.” Applied Statistics 43 (1): 117. https://doi.org/10.2307/2986116.

Haddon, Malcolm. 2011. Modelling and Quantitative Methods in Fisheries. Chapman; Hall/CRC. https://doi.org/10.1201/9781439894170.

Kristensen, Kasper, Anders Nielsen, Casper W. Berg, Hans Skaug, and Bradley M. Bell. 2016. “TMB: Automatic Differentiation and Laplace Approximation.” Journal of Statistical Software 70 (5). https://doi.org/10.18637/jss.v070.i05.

Methot, Richard D., and Chantell R. Wetzel. 2013. “Stock Synthesis: A Biological and Statistical Framework for Fish Stock Assessment and Fishery Management.” Fisheries Research 142 (May): 86–99. https://doi.org/10.1016/j.fishres.2012.10.012.

Nielsen, Anders, and Casper W. Berg. 2014. “Estimation of Time-Varying Selectivity in Stock Assessments Using State-Space Models.” Fisheries Research 158 (October): 96–101. https://doi.org/10.1016/j.fishres.2014.01.014.

Stock, Brian C., and Timothy J. Miller. 2021. “The Woods Hole Assessment Model (WHAM): A General State-Space Assessment Framework That Incorporates Time- and Age-Varying Processes via Random Effects and Links to Environmental Covariates.” Fisheries Research 240 (August): 105967. https://doi.org/10.1016/j.fishres.2021.105967.

Sullivan, Patrick J. 1992. “A Kalman Filter Approach to Catch-at-Length Analysis.” Biometrics 48 (1): 237. https://doi.org/10.2307/2532752.