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

Dr. Giancarlo M. Correa

Cousteau Consultant Group

WHAM: Abundancia a la edad

Abundancia a la edad

Transiciones de la abundancia entre edades (\(N_{a,y}\)) (Stock and Miller 2021):

\[\small{ log(N_{a,y}) = \begin{cases} log(f(SSB_{y-1})) + \varepsilon_{1,y},\hspace{5cm}a=1\\ log(N_{a-1,y-1})-Z_{a-1,y-1}+\varepsilon_{a,y},\hspace{1.4cm}1<a<A\\ log(N_{A-1,y-1}exp(-Z_{A-1,y-1}) + \\ \hspace{1cm}N_{A,y-1}exp(-Z_{A,y-1}))+\varepsilon_{A,y},\hspace{2cm}a=A \end{cases} }\]

Donde \(a\) es la edad, \(y\) es año, \(A\) es el grupo plus, \(Z\) es la mortalidad total. \(\varepsilon_{a,y}\) son desviaciones estocásticas (e.g., immigración, emmigración, misreported catch).

Stock reclutamiento

¿Qué opciones tenemos para estimar reclutamiento? Vamos a omitir las desviaciones \(\varepsilon_{1,y}\):

Estimar \(N_{1,y}\) como efectos fijos (1 parámetro por año).
Estima reclutamiento medio a lo largo de los años (\(N_{1,y}=\bar{R}\), 1 parámetro).

Stock reclutamiento

Relación Beverton-Holt (2 parámetros)

\[N_{1,y} = \frac{\alpha SSB_{y-1}}{1+\beta SSB_{y-1}}\]

Relación Ricker (2 parámetros)

\[N_{1,y} = \alpha SSB_{y-1}exp(-\beta SSB_{y-1})\]

Stock reclutamiento

Para el caso de Beverton-Holt (B-H) y Ricker, es posible también usar \(h\) (steepness) y \(R_0\) (reclutamiento sin pesca).

Para B-H, \(\alpha\) y \(\beta\) es calculado como:

\[\alpha = \frac{4h}{SPR_0(1-h)}\hspace{2cm}\beta=\frac{5h-1}{(1-h)R_0SPR_0}\]

Donde \(SPR_0\) es el spawning stock biomass per recruit en condiciones sin pesca.

Stock reclutamiento

Para Ricker, \(\alpha\) y \(\beta\) es calculado como:

\[\alpha = exp(1.25log(5h)-log(SPR_0))\\ \beta=\frac{1.25log(5h)}{R_0SPR_0}\]

Biomasa desovante

\(SSB\) se calcula como:

\[SSB_y = \sum_aN_{a,y}*W_{a,y}*mat_{a,y}*exp(-Z_{a,y}*frSSB_y)\]

Donde \(W\) es peso a la edad, \(mat\) es madurez, y \(frSSB\) es la fracción del año donde ocurre el desove.

Desviaciones

\(\varepsilon_{a,y}\) son tratados como efectos aleatorios, y pueden tener diferentes estructuras.

Independientes:

\[\varepsilon_{1,y} \sim N(-\frac{\sigma_R^2}{2},\sigma_R^2)\]

donde \(\sigma_R^2\) es la varianza de reclutamiento y es estimado como efecto fijo.

Desviaciones

Correlación temporal:

\[\varepsilon_{1,y+1} \sim N(\rho_{year}\varepsilon_{1,y}-\frac{\sigma_R^2}{2(1-\rho_{year}^2)},\sigma_R^2)\]

donde \(-1< \rho_{year}<1\) es el parámetro de correlación a lo largo de los años.

Desviaciones

Separado por edades:

\[\varepsilon_{a,y} \sim \begin{cases} N(-\frac{\sigma_R^2}{2}, \sigma_R^2),\hspace{2cm}a=1\\ N(-\frac{\sigma_a^2}{2}, \sigma_a^2),\hspace{2cm}a>1 \end{cases}\]

Donde \(\sigma_a^2\) es la varianza que aplica a edades diferentes a 1.

Desviaciones

Correlacionado por edades y año

\[E \sim MVN(0, \Sigma)\]

\(E = (\varepsilon_{1,1},...,\varepsilon_{1,Y-1},\varepsilon_{2,1},...,\varepsilon_{1,Y-1},...,\varepsilon_{A,1},...,\varepsilon_{A,Y-1})'\), \(Y\) es el número de años, y \(\Sigma\) es la matriz de covarianza:

\[Cov(\varepsilon_{a,y},\tilde{\varepsilon}_{a,y}) = \frac{\sigma_a \sigma_\tilde{a}\rho_{a}^{\vert a-\tilde{a} \vert}\rho_{y}^{\vert y-\tilde{y} \vert}}{(1-\rho^2_a)(1-\rho_y^2)}\]

donde \(-1<\rho_a <1\) y \(-1<\rho_y <1\) son los coeficientes de autocorrelación (efectos fijos).

Población inicial

Recuerda que también necesitamos la estructura poblacional para \(y=1\). Tenemos dos opciones:

Estimamos una abundancia a la edad como efectos fijos (número de parámetros igual al número de edades).
Estimamos un reclutamiento inicial \(N_{1,1}\) y una mortalidad por pesca \(F_{0}\) (2 parámetros).

Definición en R

my_input = wham::prepare_wham_input(model_name = "Example_1",
                         basic_info = input_data, 
                         NAA_re = list(), # NAA parameters
                         M = list(), # M parameter
                         selectivity = list(), # Selectivity parameter
                         catchability = list(), # Catchability parameter
                         ecov = list(), # Environmental information
                         # Crecimiento somatico:
                         growth = list(), LAA = list(), # Mean length-at-age
                         LW = list(), # Length-weight relationship
                         WAA = list(), # Mean weight-at-age
                         age_comp = "multinomial", # Age composition model
                         len_comp = "multinomial" # Length composition model
                         )

Definición en R

Enfoquemos solo la parte de abundancia a la edad. Se pueden incluir los siguientes argumentos:

my_input = wham::prepare_wham_input(...,
                         NAA_re = list(N1_model = ..., N1_pars = ..., # Pob inicial
                                       recruit_model = ..., # recruitment model
                                       use_steepness = ..., # use steepness ?
                                       recruit_pars = ...,  # S-R parameters
                                       sigma = ..., cor = ...),
                         ...
                         )

donde ... representa objetos que vamos a especificar.

Población inicial

N1_model: (integer) como modelar la abundancia a la edad inicial (primer año)
- 0: estima abundancia a la edad inicial como efectos fijos
- 1: estima reclutamiento inicial y mortalidad por pesca inicial (2 parámetros como efectos fijos)

Población inicial

N1_pars: (vector) valores iniciales para la abundancia a la edad inicial. Puede ser un vector con longitud igual al número de edades (si N1_model = 0), o un vector con dos valores: \(N_{1,1,}\) y \(F_0\) (si N1_model = 1).

Modelo de reclutamiento

recruit_model: (integer) modelo para reclutamiento
- 1: reclutamientos estimados como efectos fijos
- 2: estima reclutamiento medio
- 3: relación Beverton-Holt
- 4: relación Ricker

Modelo de reclutamiento

use_steepness: (logical) TRUE or FALSE si se desea usar steepness formulation. Solo se usa si recruit_model>2.

Modelo de reclutamiento

recruit_pars: (vector) valores iniciales para el modelo de reclutamiento.
- Si recruit_model=2, entonces viene a ser reclutamiento medio.
- Si use_steepness = FALSE y recruit_model>2, entonces viene a ser \(\alpha\) y \(\beta\).
- Si use_steepness = TRUE y recruit_model>2, entonces viene a ser \(h\) and \(R_0\).

Desviaciones

sigma: (caracter) permite estimar diferentes \(\sigma_{a}\) para diferentes grupos de edades.
- rec: solo una varianza reclutamiento (edad 1, \(\sigma^2_R\))
- rec+1: dos varianzas: \(\sigma^2_R\) y \(\sigma^2_a\) (compartida entre edades diferentes a 1).

Desviaciones

cor: (caracter) estructura de desviaciones
- iid: desviaciones independientes
- ar1_a: desviaciones correlacionadas por edad
- ar1_y: desviaciones correlacionadas por año
- 2dar1: desviaciones correlacionadas por edad y año

WHAM: Mortalidad natural

Mortalidad natural

En WHAM, tenemos 3 opciones para estimar \(M\) como efecto fijo:

\(M\) único para todas las edades, por lo tanto: \(Z_{a,y} = F_yS_{a,y} + M\).
\(M_a\) por edad, por lo tanto: \(Z_{a,y} = F_yS_{a,y} + M_{a}\).
\(M\) como función de peso a la edad (Lorenzen 1996): \(M_{y,a}=MW_{a,y}^b\). Por lo tanto: \(Z_{a,y} = F_yS_{a,y} + M_{a,y}\).

Desviaciones

Al igual que para abundancia a la edad, podemos estimar desviaciones del parámetro \(M\):

\[log(M_{a,y}) = \mu_{M_a} + \delta_{a,y}\]

\(\mu_{M_a}\) es el efecto fijo en log-scale.

Desviaciones

\(\delta_{a,y}\) puede tener diferentes estructuras asumiendo una distribución normal:

Independientes (por año y edad)
Correlacionado por edades (constante en tiempo)
Correlacionado por años (constante en edades)
Correlacionado por edades y años

Definición en R

my_input = wham::prepare_wham_input(model_name = "Example_1",
                         basic_info = input_data, 
                         NAA_re = list(), # NAA parameters
                         M = list(), # M parameter
                         selectivity = list(), # Selectivity parameter
                         catchability = list(), # Catchability parameter
                         ecov = list(), # Environmental information
                         # Crecimiento somatico:
                         growth = list(), LAA = list(), # Mean length-at-age
                         LW = list(), # Length-weight relationship
                         WAA = list(), # Mean weight-at-age
                         age_comp = "multinomial", # Age composition model
                         len_comp = "multinomial" # Length composition model
                         )

Definición en R

Enfoquemos solo la parte de abundancia a la edad. Se pueden incluir los siguientes argumentos:

my_input = wham::prepare_wham_input(...,
                         M = list(model = ..., # Modelo M 
                                  re = ..., # Estructura desviaciones
                                  initial_means = ..., # Valor inicial
                                  est_ages = ..., # Edades a estimar
                                  logb_prior = ...), # Prior b (link to W)
                         ...
                         )

donde ... representa objetos que vamos a especificar.

Modelo para M

model: (caracter) modelo M a utilizar
- constant: único M para todas las edades
- age-specific: M para cada edad
- weight-at-age: relación con peso a la edad

Estructura de desviaciones

re: (caracter) estructura a utilizar
- none: no incluir desviaciones
- iid: desviaciones por año y edad y son independientes.
- ar1_a: desviaciones correlacionadas por edad (constante en tiempo)
- ar1_y: desviaciones correlacionadas por año (constante en edad)
- 2dar1: desviaciones correlacionadas por edad y año

Valores iniciales

initial_means: (vector) valores iniciales de \(M\). Longitud va a depender del model que hayamos elegido. Para constant y weight-at-age, será solo un valor de \(M\) (\(log_b\) es asumido ser 0.305).

est_ages: (vector) indexación de edades a estimar \(M\). Si model es constant, entonces el valor 1. Si model es age-specific y si queremos estimar edad 1, 2, y 4 (sin importar el grupo plus), entonces será c(1,2,4).

logb_prior: (logical) Se debe usar prior para la relación con peso a la edad? El prior para \(b\) es \(N(0.305,0.08)\).

WHAM: Capturabilidad

Capturabilidad

Definición: proporción de peces de la población pescada por una unidad de esfuerzo pesquero. Mide la interacción entre la abundancia del recurso y esfuerzo pesquero (Arreguin-Sanchez 1996).

En WHAM, especificamos un parámetro de capturabilidad por índice de abundancia. Además, podemos incluir desviaciones:

\[logit(q_{i,y}) = logit(q_i)+\delta_{i,y}\]

Donde \(\delta_{i,y}\) sigue una distribución normal y puede ser independientes en el tiempo o correlacionadas en el tiempo.

Definición en R

my_input = wham::prepare_wham_input(model_name = "Example_1",
                         basic_info = input_data, 
                         NAA_re = list(), # NAA parameters
                         M = list(), # M parameter
                         selectivity = list(), # Selectivity parameter
                         catchability = list(), # Catchability parameter
                         ecov = list(), # Environmental information
                         # Crecimiento somatico:
                         growth = list(), LAA = list(), # Mean length-at-age
                         LW = list(), # Length-weight relationship
                         WAA = list(), # Mean weight-at-age
                         age_comp = "multinomial", # Age composition model
                         len_comp = "multinomial" # Length composition model
                         )

Definición en R

Enfoquemos solo la parte de abundancia a la edad. Se pueden incluir los siguientes argumentos:

my_input = wham::prepare_wham_input(...,
                         catchability = list(re = ..., # Estructura de desviaciones
                                             initial_q = ..., # Valores iniciales para Q
                                             q_lower = ..., # Lower bound
                                             q_upper = ..., # Lower bound
                                             prior_sd = ...), # SD para prior (normal)
                         ...
                         )

donde ... representa objetos que vamos a especificar.

Estructura de desviaciones

re: (character vector) estructura de desviaciones. Longitud debe ser igual al número de índices.
- none: no desviaciones
- iid: independientes
- ar1: correlacionados en el tiempo

Parámetros de capturabilidad

initial_q: (numeric vector) valores iniciales de q. Longitud debe ser igual al número de índices.
q_lower: (numeric vector) lower bound de q. Longitud debe ser igual al número de índices.
q_upper: (numeric vector) upper bound de q. Longitud debe ser igual al número de índices.
prior_sd: (numeric vector) valores de SD si prior quiere ser usado. Especificar NA si no se quiere usar prior. Longitud debe ser igual al número de índices.

Referencias

Arreguin-Sanchez, Francisco. 1996. “Catchability: A Key Parameter for Fish Stock Assessment.” Reviews in Fish Biology and Fisheries 6 (2). https://doi.org/10.1007/bf00182344.

Lorenzen, K. 1996. “The Relationship Between Body Weight and Natural Mortality in Juvenile and Adult Fish: A Comparison of Natural Ecosystems and Aquaculture.” Journal of Fish Biology 49 (4): 627–42. https://doi.org/10.1111/j.1095-8649.1996.tb00060.x.

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.