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