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

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

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

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

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

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

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