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.

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:

\(u2.A[j]\) representing the number of adipose clipped fish known to be hatchery fish, and
\(u2.N[j]\) representing the number of unclipped fish which are mixture of hatchery and wild fish.

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:

No fish were tagged and released in the first stratum. So no information is available to estimate the capture efficiency at the second trap in the first week.
In one of the recovery strata, there were no recaptures and so the estimated recapture probability will be zero. But the data shows that some unmarked fish were captured in these strata, so the actual efficiency must have been non-zero.
There one julian week where the number of unmarked fish captured suddenly jumps by several orders of magnitude. This jump correspond to releases of hatchery fish into the system.
Before the hatchery release, we know that all of the unmarked fish without an adipose clip are wild fish.

3.2 Fitting the BTSPAS diagonal model

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

the hatch.after variable indicates the stratum after which hatchery fish are released
the bad.m2 variable identifies the stratum where the \(m2\) value is unusual.
the clip.frac.H variable identify what fraction of the hatchery fish have adipose fin clips.

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 missing 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: 495608 
#> Random number seed for chain  268831 
#> Random number seed for chain  635035 
#> Random number seed for chain  141795 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 78
#>    Unobserved stochastic nodes: 177
#>    Total graph size: 1628
#> 
#> 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 ***

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

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)

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   4840557 2208939 2823738.3 3120941.4 3530534.4 7063046 8276518  3.589497     4
#> Utot.H 3580254 1992414 1992683.8 2130237.0 2269267.5 6129085 6796352 24.221222     3
#> Utot.W 1260303 1226999  731979.1  833369.2  945090.5 1238914 3785597  1.048881    70

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 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

\(u2.A.YoY[j]\) representing the number of YoY adipose clipped fish known to be hatchery fish, and
\(u2.N.YoY[j]\) representing the number YoY unclipped fish which are mixture of hatchery and wild fish
\(u2.A.1 [j]\) representing the number of Age1 adipose clipped fish known to be hatchery fish, and
\(u2.N.1 [j]\) representing the number of Age1 unclipped fish) which are mixture of hatchery and wild fish.

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:

No fish were tagged and released in the first two strata nor in the last stratum. So no information is available to estimate the capture efficiency at the second trap in these strata.
There one julian week where the number of unmarked fish captured suddenly jumps by several orders of magnitude. This jump correspond to releases of hatchery fish into the system.
Before the hatchery release, we know that all of the unmarked fish without an adipose clip are wild fish.

4.2 Fitting the BTSPAS diagonal model

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

the hatch.after variable indicates the stratum after which hatchery fish are released
the bad.m2 variable identifies the stratum where the \(m2\) value is unusual.
the clip.frac.H variables identify what fraction of the hatchery fish have adipose fin clips in the YoY and Age 1 fish.

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 to missing 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: 951019 
#> Random number seed for chain  339637 
#> Random number seed for chain  93791 
#> Random number seed for chain  563994 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 146
#>    Unobserved stochastic nodes: 331
#>    Total graph size: 2931
#> 
#> 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

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)

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       2409097.0 78830.5 2254767.8 2353611.5 2412489.5 2464537.7 2564559.8  1.3    12
#> Utot.1       25982.0  2801.6   21972.3   24168.0   25494.0   27165.5   33032.5  1.0   110
#> Utot.H.1      7389.3  1822.6    4959.9    5950.0    6909.5    8676.7   11383.1  2.6     4
#> Utot.H.YoY  693935.0 33621.5  624460.0  670662.7  692258.5  721546.7  749831.4  2.0     5
#> Utot.W.1     18592.7  2983.6   13499.2   16672.7   18448.0   20213.2   25440.2  1.3    10
#> Utot.W.YoY 1689180.0 64919.3 1571547.0 1642933.7 1687950.0 1731618.2 1820170.9  1.0    73
#> Utot.YoY   2383115.0 78564.3 2229938.8 2328795.5 2385771.5 2438823.2 2538630.4  1.3    12

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

\(u2.W.YoY[j]\) representing wild, YoY steelhead,
\(u2.W.1 [j]\) representing wild, age 1+ steelhead, and
\(u2.H.1 [j]\) representing hatchery, age 1+ steelhead.

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:

Marking and releasing of steelhead took place in only a few weeks! The hierarchical model will extrapolate outside these weeks to estimate the capture rate. This is likely a very dangerous thing to do.
There one julian week where the number of unmarked fish captured suddenly jumps by several orders of magnitude. This jump correspond to releases of hatchery fish into the system.
Before the hatchery release, we know that all of the unmarked fish without an adipose clip are wild fish.

5.2 Fitting the BTSPAS diagonal model

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

the hatch.after variable indicates the stratum after which hatchery fish are released
the bad.m2 variable identifies the stratum where the \(m2\) value is unusual.
the clip.frac.H variables identify what fraction of the hatchery fish have adipose fin clips in the YoY and Age 1 fish.

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 missing 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: 975987 
#> Random number seed for chain  431631 
#> Random number seed for chain  256644 
#> Random number seed for chain  265350 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 118
#>    Unobserved stochastic nodes: 235
#>    Total graph size: 3028
#> 
#> 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

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)

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       6172172 1058443.2 4525794.0 5385168 6026926 6769758 8671418 1.067168    39
#> Utot.H.1   3637034  618269.2 2665230.7 3197699 3539495 3961527 5177890 1.055898    46
#> Utot.W.1   1199063  266006.0  826481.2 1014370 1152217 1338084 1858644 1.064779    41
#> Utot.W.YoY 1336076  421846.6  709618.6 1045576 1255249 1554357 2353400 1.024553    99

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.

Diagonal Case - Multiple Stocks/Ages

Carl James Schwarz

2019-12-03

1 Location of vignette source and code.

2 Introduction

3 Wild and Hatchery Chinook

3.1 Reading in the data

3.2 Fitting the BTSPAS diagonal model

3.3 The output from the fit

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

4.1 Reading in the data

4.2 Fitting the BTSPAS diagonal model

4.3 The output from the fit

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

5.1 Reading in the data

5.2 Fitting the BTSPAS diagonal model

5.3 The output from the fit

6 References