Diagonal Case - Multiple Stocks/Ages

Carl James Schwarz

2021-10-24

1 Location of vignette source and code.

Because of the length of time needed to run the vignettes, only static vignettes have been included with this package.

The original of the vignettes and the code can be obtained from the GitHub site at https://github.com/cschwarz-stat-sfu-ca/BTSPAS

2 Introduction

In some cases, the population of fish consist of a mixture of ages (young of year, and juvenile) and stocks (wild and hatchery). BTSPAS has a number of routines to handle three common occurrences as explained below.

In all cases, only diagonal recoveries are allowed.

2.1 Fixing values of \(p\) or using covariates.

Refer to the vignette on the Diagonal Case for information about fixing values of \(p\) or modelling \(p\) using covariates such a stream flow or smoothing \(p\) using a temporal spline.

3 Wild and Hatchery Chinook

In each stratum \(j\), \(n1[j]\) fish are marked and released above a rotary screw trap. Of these, \(m2[j]\) are recaptured in the stratum of release, i.e. the matrix of releases and recoveries is diagonal. The \(n1[j]\) and \(m2[j]\) establish the capture efficiency of the trap.

At the same time, \(u2[j]\) unmarked fish are captured at the screw trap. These fish are a mixture of wild and hatchery raised Chinook Salmon. A portion (clip.rate) of the hatchery raised fish are adipose fin clipped and can be recognized as hatchery raised. The unclipped fish are a mixture of wild and hatchery fish which must be separated. Hence the \(u2[j]\) are separated into:

3.1 Reading in the data

Here is an example of some raw data that is read in:

demo.data.csv <- textConnection(
"     jweek,      n1,      m2,      u2.A,      u2.N
          9,       0,       0,         0,      4135
         10,    1465,      51,         0,     10452
         11,    1106,     121,         0,      2199
         12,     229,      25,         0,       655
         13,      20,       0,         0,       308
         14,     177,      17,         0,       719
         15,     702,      74,         0,       973
         16,     633,      94,         0,       972
         17,    1370,      62,         0,      2386
         18,     283,      10,         0,       469
         19,     647,      32,         0,       897
         20,     276,      11,         0,       426
         21,     277,      13,         0,       407
         22,     333,      15,         0,       526
         23,    3981,     242,      9427,     30542
         24,    3988,      55,      4243,     13337
         25,    2889,     115,      1646,      6282
         26,    3119,     198,      1366,      5552
         27,    2478,      80,       619,      2959
         28,    1292,      71,       258,      1455
         29,    2326,     153,       637,      3575
         30,    2528,     156,       753,      4284
         31,    2338,     275,       412,      2903
         32,    1012,     101,       173,      1127
         33,     729,      66,        91,       898
         34,     333,      44,        38,       406
         35,     269,      33,        22,       317
         36,      77,       7,         8,       99
         37,      62,       9,         2,       77
         38,      26,       3,         4,        37
         39,      20,       1,         1,        22")

demo.data <- read.csv(demo.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo.data)
#>    jweek   n1  m2 u2.A  u2.N
#> 1      9    0   0    0  4135
#> 2     10 1465  51    0 10452
#> 3     11 1106 121    0  2199
#> 4     12  229  25    0   655
#> 5     13   20   0    0   308
#> 6     14  177  17    0   719
#> 7     15  702  74    0   973
#> 8     16  633  94    0   972
#> 9     17 1370  62    0  2386
#> 10    18  283  10    0   469
#> 11    19  647  32    0   897
#> 12    20  276  11    0   426
#> 13    21  277  13    0   407
#> 14    22  333  15    0   526
#> 15    23 3981 242 9427 30542
#> 16    24 3988  55 4243 13337
#> 17    25 2889 115 1646  6282
#> 18    26 3119 198 1366  5552
#> 19    27 2478  80  619  2959
#> 20    28 1292  71  258  1455
#> 21    29 2326 153  637  3575
#> 22    30 2528 156  753  4284
#> 23    31 2338 275  412  2903
#> 24    32 1012 101  173  1127
#> 25    33  729  66   91   898
#> 26    34  333  44   38   406
#> 27    35  269  33   22   317
#> 28    36   77   7    8    99
#> 29    37   62   9    2    77
#> 30    38   26   3    4    37
#> 31    39   20   1    1    22

There are several unusual features of the data set:

3.2 Fitting the BTSPAS diagonal model

We already read in the data above. Here we set the rest of the parameters.

Don’t forget to set the working directory as appropriate:

library("BTSPAS")  

# After which weeks do the hatchery fish start to arrive. Prior to this point, all fish are wild and it is not
# necessary to separate out the wild vs hatchery
demo.hatch.after <- c(22)  # julian weeks after which hatchery fish arrive.

# Which julian weeks have "bad"  values. These will be set to 0 (releases or recaptures) or missing (unmarked) and estimated.
demo.bad.m2     <- c()   # list of julian weeks with bad m2 values
demo.bad.u2.A   <- c()     # list of julian weeks with bad u2.A values
demo.bad.u2.N   <- c()     # list of julian weeks with bad u2.N values

# The clipping fraction
demo.clip.frac.H <- .25    # what fraction of the hatchery fish are adipose fin clipped


# The prefix for the output files:
demo.prefix <- "JC-2003-CH"

# Title for the analysis
demo.title <- "Junction City 2003 Chinook - Separation of Wild and Hatchery YoY Chinook"


cat("*** Starting ",demo.title, "\n\n")
#> *** Starting  Junction City 2003 Chinook - Separation of Wild and Hatchery YoY Chinook

# Make the call to fit the model and generate the output files
demo.fit <- TimeStratPetersenDiagErrorWHChinook_fit(
                  title      =demo.title,
                  prefix     =demo.prefix,
                  time       =demo.data$jweek ,
                  n1         =demo.data$n1,
                  m2         =demo.data$m2,
                  u2.A       =demo.data$u2.A,
                  u2.N       =demo.data$u2.N ,
                  clip.frac.H= demo.clip.frac.H,
                  hatch.after=demo.hatch.after,
                  bad.m2     =demo.bad.m2, 
                  bad.u2.A   =demo.bad.u2.A, 
                  bad.u2.N   =demo.bad.u2.N,
                  debug=TRUE,  # this generates only 10,000 iterations of the MCMC chain for checking.
                  save.output.to.files=FALSE)
#> 
#> 
#> *** Start of call to JAGS 
#> Working directory:  /Users/cschwarz/Dropbox/SPAS-Bayesian/BTSPAS/vignettes 
#> Initial seed for JAGS set to: 786588 
#> Random number seed for chain  309430 
#> Random number seed for chain  41861 
#> Random number seed for chain  1207 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 78
#>    Unobserved stochastic nodes: 177
#>    Total graph size: 1630
#> 
#> Initializing model
#> 
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |++                                                |   4%
  |                                                        
  |++++                                              |   8%
  |                                                        
  |++++++                                            |  12%
  |                                                        
  |++++++++                                          |  16%
  |                                                        
  |++++++++++                                        |  20%
  |                                                        
  |++++++++++++                                      |  24%
  |                                                        
  |++++++++++++++                                    |  28%
  |                                                        
  |++++++++++++++++                                  |  32%
  |                                                        
  |++++++++++++++++++                                |  36%
  |                                                        
  |++++++++++++++++++++                              |  40%
  |                                                        
  |++++++++++++++++++++++                            |  44%
  |                                                        
  |++++++++++++++++++++++++                          |  48%
  |                                                        
  |++++++++++++++++++++++++++                        |  52%
  |                                                        
  |++++++++++++++++++++++++++++                      |  56%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  60%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  64%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  68%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  72%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  76%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  84%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  88%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  92%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  96%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |**                                                |   4%
  |                                                        
  |****                                              |   8%
  |                                                        
  |******                                            |  12%
  |                                                        
  |********                                          |  16%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |************                                      |  24%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |******************                                |  36%
  |                                                        
  |********************                              |  40%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |************************                          |  48%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |************************************              |  72%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |**************************************************| 100%
#> 
#> 
#> *** Finished JAGS ***
#> Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

3.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available in each recovery sample

demo.fit$plots$fit.plot
Fitted spline curve

Fitted spline curve

The separation of wild and hatchery fish is evident as well has the start of the spline for the hatchery fish.

A plot of the \(logit(P)\) is

demo.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

In cases where there is no information (such as the first julian week), \(BTSPAS\) has interpolated based on the distribution of catchability in the other strata and so the credible interval is very wide.

A summary of the posterior for each parameter is also available. In particular, here are the summary statistics on the posterior sample for the total number unmarked separated by wild and hatchery origin:

demo.fit$summary[ row.names(demo.fit$summary) %in% c("Ntot","Utot","Utot.H","Utot.W"),]
#>           mean      sd      2.5%     25%       50%     75%   97.5%      Rhat n.eff
#> Utot   4886878 2201403 2876388.7 3190186 3550698.7 7225287 8138957  3.741584     3
#> Utot.H 3672542 2007050 2006948.2 2205647 2392307.5 6300015 6924714 17.726043     3
#> Utot.W 1214336 1199008  722701.6  825935  926705.4 1169012 3487257  1.072534    50

4 Wild and Hatchery Chinook with YoY and Age 1 fish.

In this example, BTSPAS allows for separating wild from hatchery Chinook salmon when Age-1 Chinook Salmon are present (residualized) from last year.

In each stratum \(j\), \(n1[j]\) fish are marked and released above a rotary screw trap. Of these, \(m2[j]\) are recaptured in the stratum of release, i.e. the matrix of releases and recoveries is diagonal. The \(n1[j]\) and \(m2[j]\) establish the capture efficiency of the trap.

At the same time, \(u2[j]\) unmarked fish are captured at the screw trap. These fish are a mixture of YoY and Age-1 wild and hatchery raised Chinook Salmon. A portion (clip.rate.H.YoY, clip.rate.H.1) of the YoY and Age1 hatchery raised fish are adipose fin clipped and can be recognized as hatchery raised. The unclipped fish are a mixture of wild and hatchery fish which must be separated. Hence the \(u2[j]\) are separated into

4.1 Reading in the data

Here is an example of some raw data that is read in:

demo2.data.csv <- textConnection(
"     jweek,      n1,      m2,      u2.A.YoY,      u2.N.YoY,      u2.A.1,      u2.N.1
          2,       0,       0,             0,            15,           0,          10
          3,       0,       0,             0,            94,           1,          90
          4,     833,      52,             0,           385,           2,         112
          5,     852,      67,             0,          1162,          12,         140
          6,    1495,      77,             0,           592,          10,         103
          7,    1356,     182,             0,          1151,           4,          94
          8,    1889,     145,             0,          2258,           7,         121
          9,    2934,      89,             0,          1123,           2,          80
         10,    1546,      53,             0,          2277,           5,          57
         11,    4001,     232,             0,          2492,           4,          27
         12,    2955,     158,             0,          1579,          14,          88
         13,     529,      14,             0,          1046,           5,          45
         14,    1172,      49,             0,           766,           3,          13
         15,    3204,     232,             0,          2702,           1,           9
         16,    1328,      57,             0,         10408,           2,          18
         17,    3540,     114,             0,         12145,           3,          15
         18,    4791,      45,             0,           186,           0,           1
         19,    4808,      11,             0,           407,           0,           2
         20,    5952,      44,             0,           862,           0,           0
         21,    3852,      55,             0,           465,           0,           0
         22,    2621,      17,             0,           724,           0,          27
         23,    2131,      37,           854,          4860,           0,           0
         24,    5002,     152,           794,          3539,           0,           1
         25,    3706,     120,           904,          4597,           0,           0
         26,    1225,      44,           708,          3819,           0,           0
         27,     723,      45,           762,          3300,           0,           0
         28,    2895,     167,          1356,          5460,           0,           0
         29,    1395,     117,           614,          2918,           0,           0
         30,     479,      77,           420,          2252,           0,           0
         31,     964,      74,           289,          1240,           0,           0
         32,    2803,     288,            87,           428,           0,           0
         33,     952,      51,           114,           464,           0,           0
         34,     880,     126,            53,           515,           0,           0
         35,       0,       0,            21,            93,           0,           0")

demo2.data <- read.csv(demo2.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo2.data)
#>    jweek   n1  m2 u2.A.YoY u2.N.YoY u2.A.1 u2.N.1
#> 1      2    0   0        0       15      0     10
#> 2      3    0   0        0       94      1     90
#> 3      4  833  52        0      385      2    112
#> 4      5  852  67        0     1162     12    140
#> 5      6 1495  77        0      592     10    103
#> 6      7 1356 182        0     1151      4     94
#> 7      8 1889 145        0     2258      7    121
#> 8      9 2934  89        0     1123      2     80
#> 9     10 1546  53        0     2277      5     57
#> 10    11 4001 232        0     2492      4     27
#> 11    12 2955 158        0     1579     14     88
#> 12    13  529  14        0     1046      5     45
#> 13    14 1172  49        0      766      3     13
#> 14    15 3204 232        0     2702      1      9
#> 15    16 1328  57        0    10408      2     18
#> 16    17 3540 114        0    12145      3     15
#> 17    18 4791  45        0      186      0      1
#> 18    19 4808  11        0      407      0      2
#> 19    20 5952  44        0      862      0      0
#> 20    21 3852  55        0      465      0      0
#> 21    22 2621  17        0      724      0     27
#> 22    23 2131  37      854     4860      0      0
#> 23    24 5002 152      794     3539      0      1
#> 24    25 3706 120      904     4597      0      0
#> 25    26 1225  44      708     3819      0      0
#> 26    27  723  45      762     3300      0      0
#> 27    28 2895 167     1356     5460      0      0
#> 28    29 1395 117      614     2918      0      0
#> 29    30  479  77      420     2252      0      0
#> 30    31  964  74      289     1240      0      0
#> 31    32 2803 288       87      428      0      0
#> 32    33  952  51      114      464      0      0
#> 33    34  880 126       53      515      0      0
#> 34    35    0   0       21       93      0      0

There are several unusual features of the data set:

4.2 Fitting the BTSPAS diagonal model

We already read in the data above. Here we set the rest of the parameters.

Don’t forget to set the working directory as appropriate:

library("BTSPAS")  

# After which weeks do the YoY hatchery fish start to arrive. 
# Prior to this point, all YoY fish are wild and it is not
# necessary to separate out the YoY wild vs hatchery
demo2.hatch.after.YoY <- c(22)  # julian weeks after which YoY hatchery fish arrive.

# Which julian weeks have "bad"  values. These will be set 0 (releases or recaptures) or  missing (unmarked) and estimated.
demo2.bad.m2       <- c()     # list of julian weeks with bad m2 values
demo2.bad.u2.A.YoY <- c()     # list of julian weeks with bad u2.A.YoY values
demo2.bad.u2.N.YoY <- c()     # list of julian weeks with bad u2.N.YoY values
demo2.bad.u2.A.1   <- c()     # list of julian weeks with bad u2.A.YoY values
demo2.bad.u2.N.1   <- c()     # list of julian weeks with bad u2.N.YoY values

# The clipping fraction for the current YoY and last year's YoY (which are now Age 1 fish)
demo2.clip.frac.H.YoY <- .25    # what fraction of the YoY  hatchery fish are adipose fin clipped
demo2.clip.frac.H.1   <- .25    # what fraction of the Age1 hatchery fish are adipose fin clipped

# The prefix for the output files:
demo2.prefix <- "NF-2009-CH-WH-YoY-Age1" 

# Title for the analysis
demo2.title <- "North Fork 2009 Chinook - Separation of YoY and Age 1 Wild and Hatchery Chinook"


cat("*** Starting ",demo2.title, "\n\n")
#> *** Starting  North Fork 2009 Chinook - Separation of YoY and Age 1 Wild and Hatchery Chinook

# Make the call to fit the model and generate the output files
demo2.fit  <- TimeStratPetersenDiagErrorWHChinook2_fit(
                  title    = demo2.title,
                  prefix   = demo2.prefix,
                  time     = demo2.data$jweek,
                  n1       = demo2.data$n1,          
                  m2       = demo2.data$m2,         
                  u2.A.YoY = demo2.data$u2.A.YoY,
                  u2.N.YoY = demo2.data$u2.N.YoY,
                  u2.A.1   = demo2.data$u2.A.1,
                  u2.N.1   = demo2.data$u2.N.1,
                  clip.frac.H.YoY= demo2.clip.frac.H.YoY,
                  clip.frac.H.1  = demo2.clip.frac.H.1,
                  hatch.after.YoY= demo2.hatch.after.YoY,
                  bad.m2       = demo2.bad.m2, 
                  bad.u2.A.YoY = demo2.bad.u2.A.YoY, 
                  bad.u2.N.YoY = demo2.bad.u2.N.YoY,
                  bad.u2.A.1   = demo2.bad.u2.A.1  , 
                  bad.u2.N.1   = demo2.bad.u2.N.1,
                  debug=TRUE,  # this generates only 10,000 iterations of the MCMC chain for checking.
                  save.output.to.files=FALSE)
#> 
#> 
#> *** Start of call to JAGS 
#> Working directory:  /Users/cschwarz/Dropbox/SPAS-Bayesian/BTSPAS/vignettes 
#> Initial seed for JAGS set to: 370578 
#> Random number seed for chain  669181 
#> Random number seed for chain  956873 
#> Random number seed for chain  878382 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 146
#>    Unobserved stochastic nodes: 331
#>    Total graph size: 2933
#> 
#> Initializing model
#> 
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |++                                                |   4%
  |                                                        
  |++++                                              |   8%
  |                                                        
  |++++++                                            |  12%
  |                                                        
  |++++++++                                          |  16%
  |                                                        
  |++++++++++                                        |  20%
  |                                                        
  |++++++++++++                                      |  24%
  |                                                        
  |++++++++++++++                                    |  28%
  |                                                        
  |++++++++++++++++                                  |  32%
  |                                                        
  |++++++++++++++++++                                |  36%
  |                                                        
  |++++++++++++++++++++                              |  40%
  |                                                        
  |++++++++++++++++++++++                            |  44%
  |                                                        
  |++++++++++++++++++++++++                          |  48%
  |                                                        
  |++++++++++++++++++++++++++                        |  52%
  |                                                        
  |++++++++++++++++++++++++++++                      |  56%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  60%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  64%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  68%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  72%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  76%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  84%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  88%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  92%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  96%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |**                                                |   4%
  |                                                        
  |****                                              |   8%
  |                                                        
  |******                                            |  12%
  |                                                        
  |********                                          |  16%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |************                                      |  24%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |******************                                |  36%
  |                                                        
  |********************                              |  40%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |************************                          |  48%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |************************************              |  72%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |**************************************************| 100%
#> 
#> 
#> *** Finished JAGS ***

4.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available of both stocks and ages.

demo2.fit$plots$fit.plot
Fitted spline curve

Fitted spline curve

The separation of wild and hatchery fish is evident as well has the start of the spline for the hatchery fish.

A plot of the \(logit(P)\) is

demo2.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

In cases where there is no information (such as the first julian week), \(BTSPAS\) has interpolated based on the distribution of catchability in the other strata and so the credible interval is very wide.

A summary of the posterior for each parameter is also available. In particular, here are the summary statistics on the posterior sample for the total number unmarked separated by wild and hatchery origin and the different ages:

round(demo2.fit$summary[ grepl("Utot", row.names(demo2.fit$summary)),],1)
#>                 mean      sd      2.5%       25%       50%       75%     97.5% Rhat n.eff
#> Utot       2386940.8 84166.8 2237244.0 2327203.5 2384582.0 2437539.0 2558058.1  1.1    50
#> Utot.1       25942.0  2626.9   22200.7   24138.7   25525.5   27284.5   32443.7  1.0  1500
#> Utot.H.1      6152.0   839.7    4687.4    5553.0    6132.0    6646.5    8014.0  1.0   240
#> Utot.H.YoY  691196.5 36095.7  635368.9  667992.2  685298.0  706758.0  774673.6  1.2    24
#> Utot.W.1     19790.1  2576.2   15872.2   18034.5   19437.0   21026.2   25996.6  1.0   890
#> Utot.W.YoY 1669802.3 66913.1 1547504.9 1624302.5 1667991.5 1713981.7 1806870.2  1.0   220
#> Utot.YoY   2360998.7 83896.9 2213493.1 2301055.7 2358232.0 2411267.0 2531455.4  1.1    50

5 Wild and Hatchery Steelhead with YoY and Age 1 fish.

In this analysis we fit a diagonal time-stratified Petersen estimator separating wind from hatchery Steelhead salmon..

This differs from the Wild vs Hatchery Chinook salmon in previous sections in that all hatchery raised steelhead are marked, so there is complete separation by age and (wild/hatchery). There are 3 population of interest, Wild.YoY, Hatchery.Age1+, and Wild.Age1+.

This analysis is based on the analysis of California Junction City 2003 Steelhead data and is the example used in the Trinity River Project.

In each stratum \(j\), \(n1[j]\) fish are marked and released above a rotary screw trap. Of these, \(m2[j]\) are recaptured in the stratum of release, i.e. the matrix of releases and recoveries is diagonal. The \(n1[j]\) and \(m2[j]\) establish the capture efficiency of the trap.

At the same time, \(u2[j]\) unmarked fish are captured at the screw trap. These fish are a mixture of wild and hatchery raised steelhead salmon. The \(u2[j]\) are separated into

5.1 Reading in the data

Here is an example of some raw data that is read in:

demo3.data.csv <- textConnection(
"     jweek,      n1,      m2,      u2.W.YoY,      u2.W.1,      u2.H.1
          9,       0,       0,             0,          58,           0
         10,       0,       0,             0,         357,           2
         11,       0,       0,             0,         720,           0
         12,     999,       5,             0,         850,        4643
         13,    1707,      13,            11,         585,        5758
         14,    1947,      39,             0,         532,        4220
         15,    2109,       7,             0,         873,        2328
         16,     972,       1,             0,         303,        1474
         17,     687,       0,             1,         291,         875
         18,       0,       0,            33,          12,          39
         19,       0,       0,            31,         101,          15
         20,       0,       0,            11,          47,          13
         21,       0,       0,            78,          49,          26
         22,       0,       0,            46,          44,          22
         23,       0,       0,            35,          50,          59
         24,       0,       0,            30,          38,          15
         25,       0,       0,           309,          58,           8
         26,       3,       0,           278,          36,           4
         27,       0,       0,           207,          13,           2
         28,       0,       0,           196,           5,           0
         29,       0 ,      0,           613,          12,           0
         30,       0,       0,           764,          15,           0
         31,       0,       0,           556,          11,           0
         32,       0,       0,           250,          12,           0
         33,       0,       0,           106,          13,           0
         34,       0,       0,           413,          12,           0
         35,       0,       0,           995,          28,           1
         36,       0,       0,           357,          10 ,          0
         37,       0,       0,           181,           8,          27
         38,       0,       0,            53,           3,           2
         39,       0,       0,            29,           2,           0
         40,       0,       0,             3,           0,           0
         41,       0,       0,             5,           0,           0
         42,       0,       0,            14,           4,           0
         43,       0,       0,             8,          10,           0
         44,       0,       0,            19,           7,           0
         45,       0,       0,            46,           4,           0
         46,       0,       0,           229,           7,          0")

demo3.data <- read.csv(demo3.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo3.data)
#>    jweek   n1 m2 u2.W.YoY u2.W.1 u2.H.1
#> 1      9    0  0        0     58      0
#> 2     10    0  0        0    357      2
#> 3     11    0  0        0    720      0
#> 4     12  999  5        0    850   4643
#> 5     13 1707 13       11    585   5758
#> 6     14 1947 39        0    532   4220
#> 7     15 2109  7        0    873   2328
#> 8     16  972  1        0    303   1474
#> 9     17  687  0        1    291    875
#> 10    18    0  0       33     12     39
#> 11    19    0  0       31    101     15
#> 12    20    0  0       11     47     13
#> 13    21    0  0       78     49     26
#> 14    22    0  0       46     44     22
#> 15    23    0  0       35     50     59
#> 16    24    0  0       30     38     15
#> 17    25    0  0      309     58      8
#> 18    26    3  0      278     36      4
#> 19    27    0  0      207     13      2
#> 20    28    0  0      196      5      0
#> 21    29    0  0      613     12      0
#> 22    30    0  0      764     15      0
#> 23    31    0  0      556     11      0
#> 24    32    0  0      250     12      0
#> 25    33    0  0      106     13      0
#> 26    34    0  0      413     12      0
#> 27    35    0  0      995     28      1
#> 28    36    0  0      357     10      0
#> 29    37    0  0      181      8     27
#> 30    38    0  0       53      3      2
#> 31    39    0  0       29      2      0
#> 32    40    0  0        3      0      0
#> 33    41    0  0        5      0      0
#> 34    42    0  0       14      4      0
#> 35    43    0  0        8     10      0
#> 36    44    0  0       19      7      0
#> 37    45    0  0       46      4      0
#> 38    46    0  0      229      7      0

There are several unusual features of the data set:

5.2 Fitting the BTSPAS diagonal model

We already read in the data above. Here we set the rest of the parameters.

Don’t forget to set the working directory as appropriate:

library("BTSPAS")  

# After which weeks do the hatchery fish start to arrive. Prior to this point, all fish are wild and it is not
# necessary to separate out the wild vs hatchery
demo3.hatch.after <- c(11)  # julian weeks after which hatchery fish arrive.

# Which julian weeks have "bad"  values. These will be set to 0 (releases and recapture) or missing (unmarked captured) and estimated.
demo3.bad.m2       <- c()     # list of julian weeks with bad m2 values
demo3.bad.u2.W.YoY <- c()     # list of julian weeks with bad u2.W.YoY values
demo3.bad.u2.W.1   <- c()     # list of julian weeks with bad u2.W.1   values
demo3.bad.u2.H.1   <- c()     # list of julian weeks with bad u2.H.1   values



# The prefix for the output files:
demo3.prefix <- "demo-JC-2003-ST-TSPDE-WH" 

# Title for the analysis
demo3.title <- "Junction City 2003 Steelhead - Separation of Wild and Hatchery YoY and Age 1+ Steelhead"



cat("*** Starting ",demo3.title, "\n\n")
#> *** Starting  Junction City 2003 Steelhead - Separation of Wild and Hatchery YoY and Age 1+ Steelhead

# Make the call to fit the model and generate the output files
demo3.fit <- TimeStratPetersenDiagErrorWHSteel_fit(
                  title     = demo3.title,
                  prefix    = demo3.prefix,
                  time      = demo3.data$jweek,
                  n1        = demo3.data$n1,
                  m2        = demo3.data$m2,
                  u2.W.YoY  = demo3.data$u2.W.YoY, 
                  u2.W.1    = demo3.data$u2.W.1, 
                  u2.H.1    = demo3.data$u2.H.1,
                  hatch.after=demo3.hatch.after,
                  bad.m2     = demo3.bad.m2,
                  bad.u2.W.YoY= demo3.bad.u2.W.YoY,
                  bad.u2.W.1  = demo3.bad.u2.W.1,
                  bad.u2.H.1  = demo3.bad.u2.H.1,
                  debug=TRUE,  # this generates only 10,000 iterations of the MCMC chain for checking.
                  save.output.to.files=FALSE)
#> 
#> 
#> *** Start of call to JAGS 
#> Working directory:  /Users/cschwarz/Dropbox/SPAS-Bayesian/BTSPAS/vignettes 
#> Initial seed for JAGS set to: 830315 
#> Random number seed for chain  665415 
#> Random number seed for chain  477144 
#> Random number seed for chain  618123 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 118
#>    Unobserved stochastic nodes: 235
#>    Total graph size: 3030
#> 
#> Initializing model
#> 
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |++                                                |   4%
  |                                                        
  |++++                                              |   8%
  |                                                        
  |++++++                                            |  12%
  |                                                        
  |++++++++                                          |  16%
  |                                                        
  |++++++++++                                        |  20%
  |                                                        
  |++++++++++++                                      |  24%
  |                                                        
  |++++++++++++++                                    |  28%
  |                                                        
  |++++++++++++++++                                  |  32%
  |                                                        
  |++++++++++++++++++                                |  36%
  |                                                        
  |++++++++++++++++++++                              |  40%
  |                                                        
  |++++++++++++++++++++++                            |  44%
  |                                                        
  |++++++++++++++++++++++++                          |  48%
  |                                                        
  |++++++++++++++++++++++++++                        |  52%
  |                                                        
  |++++++++++++++++++++++++++++                      |  56%
  |                                                        
  |++++++++++++++++++++++++++++++                    |  60%
  |                                                        
  |++++++++++++++++++++++++++++++++                  |  64%
  |                                                        
  |++++++++++++++++++++++++++++++++++                |  68%
  |                                                        
  |++++++++++++++++++++++++++++++++++++              |  72%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++            |  76%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++          |  80%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++        |  84%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++      |  88%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++    |  92%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++  |  96%
  |                                                        
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
#> 
  |                                                        
  |                                                  |   0%
  |                                                        
  |**                                                |   4%
  |                                                        
  |****                                              |   8%
  |                                                        
  |******                                            |  12%
  |                                                        
  |********                                          |  16%
  |                                                        
  |**********                                        |  20%
  |                                                        
  |************                                      |  24%
  |                                                        
  |**************                                    |  28%
  |                                                        
  |****************                                  |  32%
  |                                                        
  |******************                                |  36%
  |                                                        
  |********************                              |  40%
  |                                                        
  |**********************                            |  44%
  |                                                        
  |************************                          |  48%
  |                                                        
  |**************************                        |  52%
  |                                                        
  |****************************                      |  56%
  |                                                        
  |******************************                    |  60%
  |                                                        
  |********************************                  |  64%
  |                                                        
  |**********************************                |  68%
  |                                                        
  |************************************              |  72%
  |                                                        
  |**************************************            |  76%
  |                                                        
  |****************************************          |  80%
  |                                                        
  |******************************************        |  84%
  |                                                        
  |********************************************      |  88%
  |                                                        
  |**********************************************    |  92%
  |                                                        
  |************************************************  |  96%
  |                                                        
  |**************************************************| 100%
#> 
#> 
#> *** Finished JAGS ***

The final parameter (save.output.to.files) can be set to automatically to save plots and reports in files with the appropriate prefix in the working directory.

5.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available of both stocks and ages.

demo3.fit$plots$fit.plot
Fitted spline curve

Fitted spline curve

The separation of wild and hatchery fish is evident as well has the start of the spline for the hatchery fish.

A plot of the \(logit(P)\) is

demo3.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

In cases where there is no information (such as the first julian week), \(BTSPAS\) has interpolated based on the distribution of catchability in the other strata and so the credible interval is very wide.

A summary of the posterior for each parameter is also available. In particular, here are the summary statistics on the posterior sample for the total number unmarked separated by wild and hatchery origin and the different ages:

demo3.fit$summary[ grepl("Utot", row.names(demo3.fit$summary)),]
#>               mean        sd      2.5%     25%     50%     75%    97.5%     Rhat n.eff
#> Utot       6898903 1832665.1 4375291.6 5650086 6524989 7644093 12032087 1.078609    34
#> Utot.H.1   3875433  829188.6 2589298.5 3276220 3757773 4341926  5844094 1.150846    20
#> Utot.W.1   1348184  420103.9  800107.3 1073633 1260605 1514400  2424269 1.057564    46
#> Utot.W.YoY 1675286  891964.8  772846.3 1145625 1423931 1861434  4504528 1.014674   140

6 References

Bonner, S. J., & Schwarz, C. J. (2011). Smoothing population size estimates for Time-Stratified Mark–Recapture experiments Using Bayesian P-Splines. Biometrics, 67, 1498–1507. https://doi.org/10.1111/j.1541-0420.2011.01599.x

Schwarz, C. J., & Dempson, J. B. (1994). Mark-recapture estimation of a salmon smolt population. Biometrics, 50, 98–108.