Skip to contents
library(hbm4ecology)
library(rjags)
library(posterior)
library(tidyverse)
library(bayesplot)
plot_theme <-  theme_classic(base_size = 15L) + 
  theme(axis.text.x =  element_text(angle = 90),
        strip.background = element_rect(fill = "#618288"),
        strip.text = element_text(size = 15))

Motivation and context

We illustrate the flexibility of the Bayesian state-space modeling approach for stage-structured population dynamics models fitted to series of sequential observations of different nature. The example is inspired from Rivot et al. (2004).

The method is applied to a fully stage-structured model for the Atlantic salmon life cycle with multiple life histories. The model describes the dynamics of the numbers of individuals at various life stages, with a discrete annual time step. It includes nonlinear regulation and has a probabilistic structure accommodating for both environmental and demographic stochasticity. The model is fitted to a dataset resulting from the comprehensive survey of the salmon population of the Oir River (Lower Normandy, France) between 1984 and 2001. Observation models are constructed to relate the field data to the hidden states at the various life stages. The observation process corresponds essentially to capture-mark-recapture (CMR) experiments for the evaluation of migrating juvenile and spawner runs and random sampling for demographic features.

We take into consideration a detailed description of A. salmon life histories: fish from the two smolt age classes (1+ Smolts and 2+ Smolts) can either spend one or two winters at sea (1SW and 2SW in the following) before returning as spawners. One and two sea winter spawners resulting from 1+ Smolts are denoted \(Sp11\) and \(Sp12\), respectively, and those issued from 2+ Smolts are denoted \(Sp21\) and \(Sp22\), respectively. Both \(Sp11\) and \(Sp21\) are 1SW fish but with different smolt-ages, and \(Sp12\) and \(Sp22\) are 2SW fish with different smolt-ages.

Figure 1 – Atlantic salmon life cycle model including the four main life histories encountered in French populations. Eggs spawned at year t give rise to spawners Sp11 returning year t+3. Sp12 and Sp21 returning year t+4 and Sp22 returning year t+5.

Modeling

Here we adopt a slightly simplified model than the one described in the corresponding chapter of the book (Chap. 11.3, Parent and Prévost (2013)).

Process equations describing the hidden population dynamics

  1. Spawner –> Eggs
    • The number of eggs laid by female in year \(t\) is a deterministic function of the number of spawners, the proportion of females and the fecondity rate: \[ W_t = (Sp11_t + Sp21t - y^{Sp}_{2,t}) \times pf_{1,t} \times fec_1 + (Sp12_t + Sp22t - y^{Sp}_{3,t}) \times pf_{2,t} \times fec_2\] with the following distributions:
  • \(y^2_{Sp}\) and \(y^3_{Sp}\) are known number of fish removed.
  • \(Sp21_{t=1} \sim Uniform(1, 50)\) spawners issued from 1+ smolts
  • \(Sp12_{t=1} \sim Uniform(1, 50)\) spawners issued from 2+ smolts
  • \(fec_1\) and \(fec_2\) are the mean fecondity of female which is known and constant over time
  • \(pf_{1,t}\) and \(pf_{2,t}\) are the proportion of female among 1SW and 2SW fish.
  1. Eggs –> 0+ juveniles

We model this with a density dependent Ricker relationship, with lognormal error to represent the environmental variability: \[ 0+_{t+1}^* \sim Normal(\log(W_t^* \times \alpha \times e^{-\beta\times W_{t}^*}), \sigma^2) \] \(0+_{t}^* = 0+_t/h\) and \(W_{t}^* = W_t/h\) are standardized by the surface area of habitat available for juveniles production \(h = 25229m^2\). And we give the following prior distributions:

  • \(\log(\alpha) \sim Uniform(-10,10)\)
  • \(\beta \sim Normal (0,100)\)
  • \(\log (\sigma^2) \sim Uniform(-20,20)\)
  1. 0+ juveniles –> Smolts

The juveniles \(0+_{t+1}\) will survive to next spring as pre-smolts \(PSm_{t+2}\) with probability \(\gamma_{0+}\). Here, we use an approximation of \(PSm_{t+2} \sim Binomial(0+_{t+1}, \gamma_{0+}\) by saying that a fixed rate \(\gamma_{0+}\) of juveniles will survive. We will make a similar approximation for all the survival and migration steps.

  • \(PSm_{t+2} = \gamma_{0+}\times 0+_{t+1}\)
  • \(\gamma_{0+} \sim Beta(15, 15)\)

Among pre-smolts a proportion \(\theta_{Sm1, t+2}\) will migrate as \(1+\) Smolts while the remaining will stay an additional year as \(1+\) Parrs, which will survive and migrate next year with probability \(\gamma_{parr1}\). Again we use a deterministic approximation for the binomial:

  • \(Sm1_{t+2} = PSm_{t+2} \times \theta_{Sm1, t+2}\)
  • \(Parr1_{t+2} = PSm_{t+2} - Sm1_{t+2}\)
  • \(Sm2_{t+3} = Parr1_{t+2} \times \gamma_{parr1}\)
  1. Smolts –> Returning spawners

A proportion \(\gamma_{Sm1, t+2}\) of 1+ Smolts will survive the first winter at sea and will become post-smolts. Then among these post-smolts, a proportion \(\theta_{m1}\) will mature as sea-winter fish \(Sp11\) the next summer, while the rest will not mature this year (\(Res1\)).

  • \(PostSm1_{t+3} = Sm1_{t+2} \times \gamma_{Sm1, t+2}\)
  • \(Sp11_{t+3} = PostSm1_{t+3} \times \theta_{m1}\)
  • \(Res1_{t+3} = PostSm1_{t+3} - Sp11_{t+3}\)

Also for year \(t +4\):

  • \(Sp12_{t+4} = Res1_{t+3} \times \gamma_{res}\)
  • \(PostSm2_{t+4} = Sm2_{t+2} \times \gamma_{Sm2, t+3}\)
  • \(Sp21_{t+4} = PostSm2_{t+4} \times \theta_{m1}\)
  • \(Res2_{t+4} = PostSm2_{t+4} - Sp21_{t+4}\)

and year \(t + 5\):

  • \(SP22_{t+5} = Res2_{t+4} \times \gamma_{res}\)

Finally the number of spawners in year \(t\) is obtained by the mass balance equation: \[Sp_t = Sp11_t + Sp12_t + Sp21_t + Sp22_t\]

Initial state variable

We give the following prior distribution at \(t=1\) for the fish population at different stages of the life cycle:

  • \(0+_{t=1}^* \sim Uniform(0,1)\)
  • \(PSm_{t=1} \sim Uniform(1,10000)\)
  • \(Sm2_{t=1} \sim Uniform(0,300)\)
  • \(PostSm1_{t=1} \sim Uniform(0,1000)\)
  • \(PostSm2_{t=1} \sim Uniform(0,1000)\)
  • \(Sp12_{t=1}^* \sim Uniform(0,100)\)
  • \(Sp21_{t=1}^* \sim Uniform(0,50)\)

Modeling the priors on vital rates trough an exchangeable hierarchical structures

The model encompasses both demographic and environmental stochasticity by allowing an interannual random variability of the proportion of smolts migrating as 1+ Smolts and the marine survival rate of post-smolts. Interannual random variations are modeled through an exchangeable hierarchical structure (gelman+2004?), by assuming that the parameters of all years \(t\) are randomly drawn from the same probability distribution conditioned by unknown parameters common to all years. Hierarchical structures allow to probabilistically share information among the different years and can improve the estimation of key parameters ((McAllister+2004?); (MichielsenMcAllister04?); (RivotPrevost02?)). Here, an exchangeable hierarchical structure is used to mimic a between-year variability that is assumed to be random, without any particular time trends or covariates (i.e., climate) to explain the variations.

An exchangeable hierarchical structure is used to capture the between-year variability of the probability for a pre-smolt to smoltify as 1+ smolt. The \(logit\)-transform of the \(\theta_{Sm1,t}\)’s are a priori drawn in a Normal distribution with mean and variance parameters \(\phi_{1}=(\mu_{\theta_{Sm1}},\sigma_{\theta_{Sm1}})\) shared between years. The \(\theta_{Sm1,t}\)’s are then calculated by the inverse transformation \(logit^{-1}\):

  • \(logit(\theta_{Sm1,t}) \sim Normal(\mu_{\theta_{Sm1}},\sigma_{\theta_{Sm1}})\)

The same procedure was applied for the survival rate of one year river-aged post-smolts, \(\gamma_{Sm1,t}\). The survival rates of two year river-aged post-smolts, \(\gamma_{Sm2,t}\), are defined based on the reasonable a priori constraint that the survival rate of two year river-aged post-smolts is always greater than for one year river-aged post-smolts, but with the same time variations. This constraint is specified in the \(logit\) scale.

  • \(logit(\gamma_{Sm1,t}) \; \sim \; Normal(\mu_{\gamma_{Sm1}},\sigma_{\gamma_{Sm1}})\)
  • \(logit(\gamma_{Sm2,t}) \; = \; logit(\gamma_{Sm1,t}) + \delta_{\gamma}\), \(\delta_{\gamma} \geq 0\)

Finally, the \(logit\)-transform of the proportions of females in 1SW and 2SW fish, \(pf1_{t}\) and \(pf2_{t}\) are a priori independently drawn in Normal distributions with mean and variance parameters shared between years:

  • \(logit(pf1_{t}) \; \sim \; Normal(\mu_{pf1},\sigma_{pf1})\)
  • \(logit(pf2_{t}) \; \sim \; Normal(\mu_{pf2},\sigma_{pf2})\)

Informative beta distributions based on prior ecological expertise are set for \(\gamma_{0+}\) (survival rates of 0+ Juveniles) and \(\gamma_{Parr1}\) (survival rate of 1+ Parrs), with mean \(0.5\) and \(0.66\), respectively. All the remaining parameters of the population dynamics are considered constant over time with rather diffuse prior distribution.

Where we give the following prior distributions:

We give the following prior distributions for the survival and migration rates:

  • \(\theta_{m1} \sim Beta(3,2)\)
  • \(\gamma_{res} \sim Beta(3,2)\)
  • \(\gamma_{parr1} \sim Beta(20,10)\)
  • \(\mu_{\theta_{Sm1}} \sim Uniform(0, 100)\) and \(\sigma_{\theta_{Sm1}} \sim Uniform(0, 5)\)
  • \(\mu_{\theta_{Sm1}} \sim Uniform(0, 100)\) and \(\sigma_{\theta_{Sm1}} \sim Uniform(0, 5)\)
  • \(\mu_{pf1} \sim Normal(0,100)\) and \(\sigma_{pf1} \sim Uniform(0,5)\)
  • \(\mu_{pf1} \sim Normal(0,100)\) and \(\sigma_{pf2} \sim Uniform(0,5)\)

Figure 2 – DAG of the hidden dynamics for the Atlantic salmon life cycle including the observation process on both the smolt and spaner runs.

Data and observation equations

Let first have a glimpse of the data:

glimpse(SalmonLifeCycle)
#> List of 25
#>  $ n             : num 18
#>  $ surf          : num 25229
#>  $ fec1          : num 4635
#>  $ fec2          : num 7965
#>  $ Year          : int [1:18] 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 ...
#>  $ c_Sp          : num [1:18] 167 264 130 16 226 235 15 44 31 100 ...
#>  $ x_Sp1         : num [1:18] 10 37 28 3 35 31 4 0 10 17 ...
#>  $ x_Sp2         : num [1:18] 3 11 9 1 8 5 4 0 1 2 ...
#>  $ mad           : num [1:18] 154 216 93 12 183 199 7 44 20 81 ...
#>  $ rmad          : num [1:18] 12 21 5 2 12 56 2 23 4 4 ...
#>  $ r_unm         : num [1:18] 10 4 4 22 0 0 15 1 5 3 ...
#>  $ sample_Sp_age : num [1:18] 159 211 111 16 197 220 9 41 28 98 ...
#>  $ sample_Sp11   : num [1:18] 113 116 61 13 85 129 3 38 22 85 ...
#>  $ sample_Sp21   : num [1:18] 20 50 19 2 74 54 2 1 4 6 ...
#>  $ sample_Sp12   : num [1:18] 23 43 24 1 36 31 3 1 2 7 ...
#>  $ sample_Sp22   : num [1:18] 3 2 7 0 2 6 1 1 0 0 ...
#>  $ sample_Sp1_sex: num [1:18] 141 203 93 15 182 197 9 42 28 93 ...
#>  $ sample_Sp2_sex: num [1:18] 26 61 37 1 44 38 6 2 3 7 ...
#>  $ sample_Sp1f   : num [1:18] 40 69 31 1 63 64 3 13 11 22 ...
#>  $ sample_Sp2f   : num [1:18] 21 42 26 1 31 21 4 1 2 3 ...
#>  $ c_Sm          : num [1:18] NA 439 887 283 307 553 746 151 580 209 ...
#>  $ sample_Sm_age : num [1:18] 1 439 887 283 307 553 746 151 580 209 ...
#>  $ sample_Sm1    : num [1:18] NA 232 848 146 282 495 708 101 571 171 ...
#>  $ m_Sm          : num [1:18] 86 86 135 31 59 65 38 35 50 26 ...
#>  $ r_Sm          : num [1:18] NA NA 91 24 43 43 35 27 43 24 ...

Capture-mark-recapture and sampling data for smolts at the downstream trap by migration year (\(y_i^{Sm}\), \(i = 1, \dots, 5)\):

knitr::kable(bind_cols(SalmonLifeCycle[c(5, 21:25)]))
Year c_Sm sample_Sm_age sample_Sm1 m_Sm r_Sm
1984 NA 1 NA 86 NA
1985 439 439 232 86 NA
1986 887 887 848 135 91
1987 283 283 146 31 24
1988 307 307 282 59 43
1989 553 553 495 65 43
1990 746 746 708 38 35
1991 151 151 101 35 27
1992 580 580 571 50 43
1993 209 209 171 26 24
1994 329 329 323 17 10
1995 618 618 541 63 53
1996 767 767 684 76 58
1997 205 205 186 63 31
1998 511 511 438 91 44
1999 195 195 43 59 45
2000 1849 1849 1835 300 232
2001 688 688 636 264 123

Capture-mark-recapture and sampling data for spawners at the upstream trap by migration year (\(y_i^{S}\), \(i = 1, \dots, 6)\):

knitr::kable(bind_cols(SalmonLifeCycle[c(5, 6:11)]))
Year c_Sp x_Sp1 x_Sp2 mad rmad r_unm
1984 167 10 3 154 12 10
1985 264 37 11 216 21 4
1986 130 28 9 93 5 4
1987 16 3 1 12 2 22
1988 226 35 8 183 12 0
1989 235 31 5 199 56 0
1990 15 4 4 7 2 15
1991 44 0 0 44 23 1
1992 31 10 1 20 4 5
1993 100 17 2 81 4 3
1994 32 12 2 18 1 4
1995 109 6 1 102 39 7
1996 70 13 2 55 25 57
1997 56 19 3 34 12 3
1998 34 3 1 30 6 30
1999 154 5 1 148 13 22
2000 53 0 0 53 4 33
2001 160 1 0 159 31 4

Capture-mark-recapture and sampling data for spawners at the upstream trap by migration year (\(y_i^{S}\), \(i = 7, \dots, 15)\) (continued):

knitr::kable(bind_cols(SalmonLifeCycle[c(5, 12:20)]))
Year sample_Sp_age sample_Sp11 sample_Sp21 sample_Sp12 sample_Sp22 sample_Sp1_sex sample_Sp2_sex sample_Sp1f sample_Sp2f
1984 159 113 20 23 3 141 26 40 21
1985 211 116 50 43 2 203 61 69 42
1986 111 61 19 24 7 93 37 31 26
1987 16 13 2 1 0 15 1 1 1
1988 197 85 74 36 2 182 44 63 31
1989 220 129 54 31 6 197 38 64 21
1990 9 3 2 3 1 9 6 3 4
1991 41 38 1 1 1 42 2 13 1
1992 28 22 4 2 0 28 3 11 2
1993 98 85 6 7 0 93 7 22 3
1994 29 24 3 2 0 30 2 18 2
1995 108 88 17 3 0 106 3 45 3
1996 60 48 9 2 1 67 3 27 3
1997 52 47 4 1 0 55 1 21 1
1998 28 22 5 1 0 33 1 10 1
1999 140 105 18 12 5 136 18 43 13
2000 51 45 2 2 2 49 4 26 2
2001 140 120 14 5 1 151 9 51 3

The observations (1984 - 2001) have been gathered under an homogeneous experimental design and few data are missing. Each annual survey provides two complementary sources of information, the number of smolts and spawners in each annual run, and the demographic structure (age classes and sex-ratio) in smolts and adult runs.

Updating the population size at the various life stages

Available information to estimate the size of smolt and spawner runs does not distinguish between the age classes. For both the smolts and the spawners, a proportion of the migrating population is captured at a partial counting fence and marked for further recapture. CMR experiments provide data to update the total number of smolts migrating in year \(t\), \(Sm_{t} = Sm1_{t}+Sm2_{t}\), and the total number of spawners returning in year \(t\), \(Sp_{t}=Sp11_{t}+Sp12_{t}+Sp21_{t}+Sp22_{t}\).

Smolt run

Figure 3 – Schame of the Cerisel trapping facility (Oir river, Normandy, France) working as a double downstream (for smolt runs) and upstream (for spawner runs) partial counting fence. The black circles at the top of the figure indicates the release site of tagged smolts and spawners. The release site and the downstream trap are 1 km away from each other.

CMR experiments for smolts are analogous to the two-stage Petersen experiment (). For each year \(t\) between 1986 and 2001, let us denote \(y_{1,t}^{Sm}\) the number of smolts caught in the downstream trapping facility during the migration time. Among the \(y_{1,t}^{Sm}\) smolts captured, a number \(y_{2,t}^{Sm}\) have been tagged and released upstream from the trapping facility used for capture. Some of these tagged and released smolts, denoted \(y_{3,t}^{Sm}\) (\(y_{3,t}^{Sm}<y_{2,t}% ^{Sm}\)), will be recaptured at the same downstream trap. The capture and recapture data for each year \(t\) can be modeled by two successive Binomial distributions with the same capture probabilities \(\pi_{t}^{Sm}\) (interpreted as the trapping efficiency that may vary between years):

  • \(y^{Sm}_{1,t} \; \sim \; Binomial(Sm_{t},\pi^{Sm}_{t})\)
  • \(y^{Sm}_{3,t} \; \sim \; Binomial(y^{Sm}_{2,t},\pi^{Sm}_{t})\)

Intuitively, the second Binomial equation will contribute to specify the trapping efficiency \(\pi_{Sm,t}\), while the first one brings information to update the total number of downstream migrating smolts \(Sm_{t}\).

Spawner run

Conversely to smolts, CMR experiments for spawners cannot be assimilated to a simple Petersen experiment. For each year \(t\) from \(1984\) to \(2001\), \(y_{1,t}^{Sp}\) denotes the number of fish trapped at the counting fence. Among these \(y_{1,t}^{Sp}\) fish, \(y_{2,t}^{Sp}\) one sea-winter and \(y_{3,t}^{Sp}\) two sea-winter fish are removed from the population. A number \(y_{4,t}% ^{Sp}=y_{1,t}^{Sp}-(y_{2,t}^{Sp}+y_{3,t}^{Sp})\) of fish are tagged and released upstream from the trap. The recapture sample is gathered during and after spawning by three methods (electrofishing on the spawning grounds, collection of dead fish after spawning, and trapping of spent fish at the downstream trap). Let us denote as \(y_{5,t}^{Sp}\) and \(y_{6,t}^{Sp}\) the number of marked and unmarked fish among recaptured fish, respectively.

  • \(y^{Sp}_{1,t} \; \sim \; Binomial(Sp_{t},\pi^{Sp,1}_{t})\)
  • \(y^{Sp}_{5,t} \; \sim \; Binomial(y^{Sp}_{4,t},\pi^{Sp,2}_{t})\)
  • \(y^{Sp}_{6,t} \; \sim \; Binomial(Sp_{t}-y^{Sp}_{1,t},\pi^{Sp,2}_{t})\)

The reasoning underlying how the information contained in the data is used in the three lines is as follows: the second line essentially conveys information to estimate the recapture probability \(\pi _{t}^{Sp,2}\). Both the third and the first line carry information on the total number of spawners \(Sp_{t}\).

Hierarchical prior on the trapping efficiencies

The between-year variability of the probabilities of capture \(\pi_{t}^{Sm}\), \(\pi_{t}^{Sp,1}\) and \(\pi_{t}^{Sp,2}\) is modeled through exchangeable hierarchical structures with hyper-parameters denoted \(\phi_{Sm}\), \(\phi_{Sp,1}\) and \(\phi_{Sp,2}\) respectively. Instead of drawing the \(\pi\)’s in beta distributions, the \(logit\)-transform of the \(\pi\)’s are drawn from Normal distributions with mean and variance parameters shared between years. The \(\pi\)’s are then calculated by the inverse transformation \(logit^{-1}\).

  • \(logit(\pi^{Sm}_{t}) \; \sim \; Normal(\mu_{\pi_{Sm}},\sigma_{\pi_{Sm}})\)
  • \(logit(\pi^{Sp,1}_{t}) \; \sim \; Normal(\mu_{\pi_{Sp,1}},\sigma_{\pi_{Sp,1}})\)
  • \(logit(\pi^{Sp,2}_{t}) \; \sim \; Normal(\mu_{\pi_{Sp,2}},\sigma_{\pi_{Sp,2}})\)

We set diffuse prior distribution on the parameters of the Normal distributions:

  • \(\mu_{\pi_{Sm}} \sim Normal(0,100)\) and \(\sigma_{\pi_{Sm}} \sim Uniform(0,5)\)
  • \(\mu_{\pi_{Sp,1}} \sim Normal(0,100)\) and \(\sigma_{\pi_{Sp,1}} \sim Uniform(0,5)\)
  • \(\mu_{\pi_{Sp,2}} \sim Normal(0,100)\) and \(\sigma_{\pi_{Sp,2}} \sim Uniform(0,5)\)

Updating the demographic structure

In addition to the CMR data, sampling among the fish caught at the trap enables us to update the yearly proportions of the age classes in the smolt and spawner runs \(Sm_{t}\) and \(Sp_{t}\), and the mean sex ratio in the spawner run. For each demographic feature, sampling is modeled by independent Binomial or multinomial processes.

River age in smolt run

\(y^{Sm}_{4,t}\) denotes the sample size of smolts examined for river-age in the run of year \(t\) and \(y^{Sm}_{5}\) is the number of 1+ Smolts among them. The proportion of 1+ Smolts in year \(t\), denoted \(\rho_{Sm1,t}\) is updated by modeling \(y^{Sm}_{5}\) as the result of a Binomial sampling with a sample size \(y^{Sm}_{4,t}\) and a probability \(\rho_{Sm1,t}\):

  • \(\rho_{Sm1,t} = \frac{Sm1_{t}}{Sm1_{t} + Sm2_{t}}\)
  • \(y^{Sm}_{5,t} \; \sim \; Binomial(y^{Sm}_{4,t},\rho_{Sm1,t})\)

Sea and River-age in spawner run

\(y_{7,t}^{Sp}\) denotes the sample size of spawners examined for sea and river-age in the run of year \(t\). \(y_{8,t}^{Sp}\) and \(y_{9,t}^{Sp}\) are the number of 1SW fish resulting from 1+ and 2+ Smolts, respectively among \(y_{7,t}^{Sp}\), and \(y_{10,t}^{Sp}\) and \(y_{11,t}^{Sp}\) are the number of 2SW fish issued from 1+ and 2+ Smolts, respectively among \(y_{7,t}^{Sp}\). The proportions of the four different age classes in each spawning run (summing to \(1\)) are updated through a multinomial model:

  • $ _{Sp,t} = ( , , , )$
  • \((y^{Sp}_{8,t},y^{Sp}_{9,t},y^{Sp}_{10,t},y^{Sp}_{11,t}) \; \sim \; Multinomial(y^{Sp}_{7,t},\rho_{Sp,t})\)

Sex ratio in spawner run

Sex-ratio in spawner runs are updated each year \(t\) through Binomial sampling experiments considered independent for one and two sea-winter fish. Each year \(t\), \(y^{Sp}_{12,t}\) 1SW fish (resp. \(y^{Sp}_{14,t}\) 2SW fish) are examined for sex identification and \(y^{Sp}_{13,t}\) (resp. \(y^{Sp}_{15,t}\)) are identified as females among them. Proportion of females in the two sea-age classes \(pf_{1,t}\) and \(pf_{2,t}\) are updated through Binomial sampling:

  • $ y^{Sp}{13,t} Binomial(y^{Sp}{12,t},pf_{1,t})$
  • $ y^{Sp}{14,t} Binomial(y^{Sp}{15,t},pf_{2,t})$

Model estimation

We will remove the data from 1984 and 1985 as there are missing data (NA).

data <- SalmonLifeCycle
data$n <-  16
data[c(5:25)] <- bind_cols(SalmonLifeCycle[c(5:25)]) %>% filter(Year >= 1986) %>% as.list()
data$sample_Sp <- bind_cols(data[13:16])
model_str <- "
model {
# ---------------------------------------------------------------------------
#                               PRIORS  
# ---------------------------------------------------------------------------

# Priors for Ricker stock-recruitmnt parameters
# R = alpha*exp(-beta*S)
# alpha constrained <1 becasue it is a survival rate eggs --> 0+

    log.alpha ~ dunif(-10,0)        
    alpha <- exp(log.alpha)
    beta ~ dnorm(0,0.01)
    exp.beta <- exp(beta)
    
# Prior for recruitment process error variance
# (Corresponds to approximative 1/V between 1E-6 and 1E+6)

    log.sigma_2 ~ dunif(-13.8,13.8)
    sigma_2 <- exp(log.sigma_2)
    tau <- 1/sigma_2

# Survival rate Juveniles 0+ --> Pre-smolts1 (constant between years)

    s01 ~ dbeta(15,15)

# Life history choice - probability to migrate as Smolt 1
# (exchangeable hierarchical structure between years)

    mu_theta_Sm1 ~ dnorm(0,0.01)
    sd_theta_Sm1  ~ dunif(0,5)
    tau_theta_Sm1 <- 1/(sd_theta_Sm1*sd_theta_Sm1)
    # Indice [n+1] is the predictive
    for( i in 1:(n+1) )     {
    logit_theta_Sm1[i] ~ dnorm(mu_theta_Sm1,tau_theta_Sm1)
    logit(theta_Sm1[i]) <- logit_theta_Sm1[i]           }

# Survival rate Parr 1+ resident --> Smolts 2 (constant between years)

    s12 ~ dbeta(20,10)

# Transition Smolt1 --> post-smolts (survival during the first year at sea)
# (exchangeable hierarchical structure between years)

    mu_ss11 ~ dnorm(0,0.01)
    sd_ss11 ~ dunif(0,5)
    tau_ss11 <- 1/(sd_ss11*sd_ss11)
    # Indice [n+1] is the predictive
    for( i in 1 : (n+1) )   {
    logit_ss11[i] ~ dnorm(mu_ss11,tau_ss11)
    logit(ss11[i]) <- logit_ss11[i]
    logit(ss21[i]) <- logit_ss11[i] + delta_ss1 }

# Differential in survival rate between Sm1 and Sm2 (constant between years)
# logit(ss21) = logit(ss11) + delta_ss1

    delta_ss1 ~ dunif(0,10)

# Proportion of maturing adults as 1SW (constant between years)
    
    theta_m1 ~ dbeta(3,2)

# Marine survival rate for non maturing 1SW post-smolts
# Survival probability is the same for non maturing post.smolts1 and post.smolts2 
# (ss12 = ss22) and constant accross year

    ss2 ~ dbeta(3,2)

# Proportion of females pf1 and pf2 (exchangeable hierarchical structure between years)
    
    mu_pf1 ~ dnorm(0,0.01)
    sd_pf1 ~ dunif(0,5)
    tau_pf1 <- 1/(sd_pf1*sd_pf1)

    mu_pf2 ~ dnorm(0,0.01)
    sd_pf2 ~ dunif(0,5)
    tau_pf2 <- 1/(sd_pf2*sd_pf2)

    # Indice [n+1] is the predictive
    for( i in 1 : (n+1) )   {
    logit_pf1[i] ~ dnorm(mu_pf1,tau_pf1)
    logit(pf1[i]) <- logit_pf1[i]
    logit_pf2[i] ~ dnorm(mu_pf2,tau_pf2)
    logit(pf2[i]) <- logit_pf2[i]       }

#  Trapping efficiencies (exchangeable hierarchical structure between years)

    # Smolts
    mu_pi_sm ~ dnorm(0,0.01)
    sd_pi_sm ~ dunif(0,5)
    tau_pi_sm <- 1/(sd_pi_sm*sd_pi_sm)
    
    # Adults
    mu_pi_sp1 ~ dnorm(0,0.01)
    sd_pi_sp1 ~ dunif(0,5)
    tau_pi_sp1 <- 1/(sd_pi_sp1*sd_pi_sp1)
    
    # Adults recapture efficiency
    mu_pi_sp2 ~ dnorm(0,0.01)
    sd_pi_sp2 ~ dunif(0,5)
    tau_pi_sp2 <- 1/(sd_pi_sp2*sd_pi_sp2)

    # Indice [n+1] is the predictive
    for( i in 1 : (n+1) )   {
    # Smolts
    logit_pi_sm[i] ~ dnorm(mu_pi_sm,tau_pi_sm)
    logit(pi_sm[i]) <- logit_pi_sm[i]
    
    # Adults
    logit_pi_sp1[i] ~ dnorm(mu_pi_sp1,tau_pi_sp1)
    logit(pi_sp1[i]) <- logit_pi_sp1[i]
    
    # Adults recapture efficiency
    logit_pi_sp2[i] ~ dnorm(mu_pi_sp2,tau_pi_sp2)
    logit(pi_sp2[i]) <- logit_pi_sp2[i]     }


# ---------------------------------------------------------------------------
#           PROCESS EQUATIOnS (hidden population dynamics)
# ---------------------------------------------------------------------------

# Prior on states for the first year
    dSp12 <- rep(1, 100)
    dSp22 <- rep(1, 50)
    
    Sp12[1] ~ dcat(dSp12)
    Sp22[1] ~ dcat(dSp22)
    
    dSm2 <- rep(1, 300)
    
    Sm2[1] ~ dcat(dSm2)
    
    Jint[1] ~ dunif(0,1)

  dPSm <- rep(1, 10000)
    PSm[1] ~ dcat(dPSm)

    dPostSm1 <- rep(1, 1000)
    dPostSm2 <- rep(1, 1000)    
    post.smolt1[1] ~ dcat(dPostSm1)
    post.smolt2[1] ~ dcat(dPostSm2)

#  Loop on time series (n = 18, i = 1 is year 1984, i = n is year 2001)

for( i in 1 : n ){
    
    Sp1[i] <- Sp11[i] + Sp21[i]
    Sp2[i] <- Sp12[i] + Sp22[i]
    Sp[i] <- Sp1[i] + Sp2[i]
    
# Stock (eggs) : number of female = marked and released + escape from the trap
# Standardized per m² of wetted production area

    Eggs[i] <- max( ( (Sp1[i]-x_Sp1[i]) * pf1[i] * fec1  + (Sp2[i]-x_Sp2[i]) * pf2[i] * fec2 )  /  surf, 1)

    p_Eggs1[i] <- max(  (Sp1[i]-x_Sp1[i]) * pf1[i] * fec1 / surf ,1) / Eggs[i] 
    
# Ricker * Lognormal Process errors for recruitment

    LogRmean[i] <- log(Eggs[i]*alpha*exp(-beta*Eggs[i]))
    Jint[i+1] ~ dlnorm(LogRmean[i],tau)
    J[i] <- round(Jint[i]*surf)
        
    egg_juv_surv[i] <- Jint[i+1]/Eggs[i]
    
# Transition Juveniles 0+ --> pre-smolts PSm with survival s01

    PSm[i+1] <-  max(round(s01 * J[i]), 1)
    
# Transition to smolts 1+ and 2+ with life history choice theta_Sm1
# Survival of resident Parr1 = s12 gamma_Parr1
    
    Sm1[i]   <- max(round(theta_Sm1[i] * PSm[i]), 1)
    Parr1[i] <- max(PSm[i] - Sm1[i],1)

    Sm2[i+1] <- max(round(s12 * Parr1[i]), 1)

# Transition Smolt1 --> post-smolts (survival during the first year at sea)
    
    post.smolt1[i+1] <- max(round(ss11[i] * Sm1[i]), 1)
    post.smolt2[i+1] <- max(round(ss21[i] * Sm2[i]), 1)

# Maturation of post-smolts (1 and 2) as 1SW fish (with probability theta_m1)
# Maturing probability (theta_m1) is the same for post.smolts1 and post.smolts2 
# and constant accross year

    Sp11[i] <- max(round(theta_m1 * post.smolt1[i]), 1)
    Sp21[i] <- max(round(theta_m1 * post.smolt2[i]), 1)

# Survival during the second year at sea
    
    Res1[i]   <- max(post.smolt1[i] - Sp11[i],1)
    Sp12[i+1] <- max(round(ss2 * Res1[i]),1)

    Res2[i]   <- max(post.smolt2[i] - Sp21[i],1)
    Sp22[i+1] <- max(round(ss2 * Res2[i]),1)

}  # end of the loop on i

mean_p_Eggs1 <- mean(p_Eggs1[])
mean_egg_juv_surv <- mean(egg_juv_surv[])


# ---------------------------------------------------------------------------
#                       OBSERVATIOn EQUATIONS
# ---------------------------------------------------------------------------

for( i in 1 : n )
{

# Smolts
# --------------------------------

# Catches

    Sm[i] <- Sm1[i]+Sm2[i]
    c_Sm[i] ~ dbin(pi_sm[i],Sm[i])

# Recaptures

    r_Sm[i] ~ dbin(pi_sm[i],m_Sm[i])

# Updating of river-age proportions in smolts runs

    p_Sm1[i] <- Sm1[i]/(Sm1[i]+Sm2[i])
    sample_Sm1[i] ~ dbin(p_Sm1[i],sample_Sm_age[i])


# Adults
# --------------------------------

# Updating total numbers
    
# Catches       

    c_Sp[i] ~ dbin(pi_sp1[i],Sp[i])

# Recaptured marked

    rmad[i] ~ dbin(pi_sp2[i],mad[i])

# Recaptured unmarked (r_unm[i] among a total of escnet[i] unmarked fish)
# where escnet[i] are adults that escaped to the trap = Sp - c_Sp

    escnet[i] <- Sp[i] - c_Sp[i]
    r_unm[i] ~ dbin(pi_sp2[i], escnet[i])

# Updating the demographic structure

    p_Sp[i,1] <- Sp11[i]/Sp[i]
    p_Sp[i,2] <- Sp12[i]/Sp[i]
    p_Sp[i,3] <- Sp21[i]/Sp[i]
    p_Sp[i,4] <- Sp22[i]/Sp[i]

    p_Sp1[i] <- Sp1[i] / Sp[i] 
    p_Sp11[i] <- Sp11[i] / Sp1[i]
    p_Sp12[i] <- Sp12[i] / Sp2[i]
    
    #x[i,1] <-  sample_Sp11[i]
    #x[i,2] <-  sample_Sp12[i]
    #x[i,3] <-  sample_Sp21[i]
    #x[i,4] <-  sample_Sp22[i]
    
    sample_Sp[i,1:4] ~ dmulti(p_Sp[i,1:4], sample_Sp_age[i])
    
# Updating female proportions in 1SW and 2SW

    sample_Sp1f[i] ~ dbin(pf1[i],sample_Sp1_sex[i])
    sample_Sp2f[i] ~ dbin(pf2[i],sample_Sp2_sex[i])

}  # end of the loop on i

    mean_p_Sm1 <- mean(p_Sm1[])

    mean_p_Sp1 <- mean(p_Sp1[])
    mean_p_Sp11 <- mean(p_Sp11[])
    mean_p_Sp12 <- mean(p_Sp12[])

}  # End of the model"
var_prod <- c("alpha", "beta", "sigma_2")
var_survival <- c("s01", "s12", "ss2", "ss11", "ss21")
var_migration <- c("theta_Sm1", "delta_ss1", "theta_m1")
var_female <- c("pf1", "pf2")
var_trapping <- c("pi_sm", "pi_sp1", "pi_sp2")
var_egg <- c("Eggs", "p_Eggs1`", "J", "LogRmean", "egg_juv_surv")
var_pop <- c("Jint", "PSm", "post.smolt1", "post.smolt2", # Juveniles and smolts
                     "Sm1", "Parr1", "Sm2", "Res1", "Sm",
                     "Sp11", "Sp12", "Sp21", "Sp22", "Sp1", "Sp2", "Sp", # Spawners
                        "mean_p_Sm1", "mean_p_Sp1", "mean_p_Sp11", "mean_p_Sp12")
inits1 <- source(file = "model_init/SalmonLifeCycle/inits1_SalmonLifeCycle.txt")$value
inits <- list(inits1, inits1, inits1) 
model <- jags.model(file = textConnection(model_str), 
                    data = data, inits = inits,
                    n.chains = 3)
# Inferences
update(model, n.iter = 20000)
posterior_sample <- coda.samples(
  model = model,
  variable.names = c("alpha", "beta", "sigma_2", # Production
                     "s01", "s12", "ss2", "ss11", "ss21", #Survival rate gamma
                     "theta_Sm1", "delta_ss1", "theta_m1", # Migration rate
                     "pf1", "pf2", # Proportion of female
                     "pi_sm", "pi_sp1", "pi_sp2", # Trapping efficiency
                     "Eggs", "p_Eggs1", "J", "LogRmean", "egg_juv_surv",
                     "Jint", "PSm", "post.smolt1", "post.smolt2", # Juveniles and smolts
                     "Sm1", "Parr1", "Sm2", "Res1", "Sm",
                     "Sp11", "Sp12", "Sp21", "Sp22", "Sp1", "Sp2", "Sp", # Spawners
                        "mean_p_Sm1", "mean_p_Sp1", "mean_p_Sp11", "mean_p_Sp12"
                     ),
                       n.iter = 30000,
                       thin = 10)
summarise_draws(as_draws_df(posterior_sample),
                default_convergence_measures()) %>%
  pivot_longer(-variable) %>% 
  ggplot(aes(x = value)) +
  geom_histogram() +
  facet_wrap(~ name, scales = "free") +
  plot_theme
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 6 rows containing non-finite values (stat_bin).

draws <- as_draws_df(posterior_sample)
mcmc_pairs(draws, 
           pars = c("s01", "s12", "theta_m1", "delta_ss1", "sigma_2"),
           diag_fun = "dens")

mcmc_pairs(draws, pars = c("alpha", "beta"),
           diag_fun = "dens")

mcmc_hist(draws, pars = c("pi_sp1[1]", "pi_sp1[2]", "pi_sp1[3]",
                          "pi_sp1[4]", "pi_sp1[5]", "pi_sp1[6]"))
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# mcmc_pairs(draws, pars = c("mu_pf1", "mu_pf2"),
#            diag_fun = "dens")
dfw <- extract_wider(posterior_sample) %>% 
  rename()
dfl <- extract_longer(posterior_sample)
dfw %>%
  select(starts_with("Sm.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Smolts (Sm)") +
  ylim(0,3000) +
  # scale_x_discrete(labels = as.character(seq(1986, 1986+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  plot_theme

dfw %>%
  select(starts_with("Sp.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Spawners (Sp)") +
  # scale_x_discrete(labels = as.character(seq(1986, 1986+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  plot_theme

dfw %>%
  select(starts_with("pi_sm.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Trapping efficiency for Smolts (pi_Sm)") +
  ylim(c(0,1)) +
  # scale_x_discrete(labels = as.character(seq(1984, 1984+17)),
                   # guide = guide_axis(angle = 90),
                   # breaks = seq(1984, 1984+17)) +
  plot_theme

dfw %>%
  select(starts_with("pi_Sp1.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab (bquote(paste("Trapping efficiency for Spawners : ", pi[Sp1]))) +
  ylim(c(0,1)) +
  plot_theme

  # scale_x_discrete(labels = as.character(seq(1984, 1984+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  # theme_classic(base_size = 15L)
dfw %>%
  select(starts_with("pi_Sp2.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab (bquote(paste("Trapping efficiency for Spawners : ", pi[Sp2]))) +
  ylim(c(0,1)) +
  plot_theme

  # scale_x_discrete(labels = as.character(seq(1984, 1984+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  # theme_classic(base_size = 15L)
name_years <- as.character(seq(1986, 1986+15))
dfw %>%
  select(starts_with("theta_Sm1.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Probability to smolitify") +
  ylim(c(0,1)) +
  plot_theme

  # scale_x_discrete(labels = name_years,
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  # theme_classic(base_size = 15L)
dfw %>%
  select(starts_with("ss11.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("1+ Post-smolts survival rate") +
  ylim(c(0,1)) +
  plot_theme

  # scale_x_discrete(labels = as.character(seq(1984, 1984+17)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1984, 1984+17)) +
  # theme_classic(base_size = 15L)
dfw %>%
  select(starts_with("Sp11.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Spawners (Sp11)") +
  plot_theme

  # scale_x_discrete(labels = as.character(seq(1986, 1986+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  # theme_classic(base_size = 15L)
dfw %>%
  select(starts_with("Sm1.")) %>%
  pivot_longer(cols = everything()) %>%
  mutate(name = as_factor(str_sub(name, 4))) %>%
  ggplot(aes(group = name, y = value, x = as.integer(name) + 1985)) +
  geom_boxplot(outlier.shape = NA, width = .5, fill = "gray70") +
  xlab("Years") + ylab ("Smolts (Sm1)") +
  plot_theme

  # scale_x_discrete(labels = as.character(seq(1986, 1986+15)),
  #                  guide = guide_axis(angle = 90),
  #                  breaks = seq(1986, 1986+15)) +
  # theme_classic(base_size = 15L)

References

Parent, E., and E. Prévost. 2013. Introduction to Hierarchical Bayesian Modeling for Ecological Data. CRC press Boca Raton, Florida.
Rivot, E., E. Prévost, E. Parent, and J. L. Baglinière. 2004. “A Bayesian State-Space Modelling Framework for Fitting a Salmon Stage-Structured Population Dynamic Model to Multiple Time Series of Field Data.” Ecological Modelling 179: 463–85.