knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(binom) library(BTSPAS) library(ggplot2) max.width=70
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
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.
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:
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)
There are several unusual features of the data set:
We already read in the data above. Here we set the rest of the parameters.
Don't forget to set the working directory as appropriate:
library("BTSPAS") # After which weeks do the hatchery fish start to arrive. Prior to this point, all fish are wild and it is not # necessary to separate out the wild vs hatchery demo.hatch.after <- c(22) # julian weeks after which hatchery fish arrive. # Which julian weeks have "bad" values. These will be set to 0 (releases or recaptures) or missing (unmarked) and estimated. demo.bad.m2 <- c() # list of julian weeks with bad m2 values demo.bad.u2.A <- c() # list of julian weeks with bad u2.A values demo.bad.u2.N <- c() # list of julian weeks with bad u2.N values # The clipping fraction demo.clip.frac.H <- .25 # what fraction of the hatchery fish are adipose fin clipped # The prefix for the output files: demo.prefix <- "JC-2003-CH" # Title for the analysis demo.title <- "Junction City 2003 Chinook - Separation of Wild and Hatchery YoY Chinook" cat("*** Starting ",demo.title, "\n\n") # 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)
# delete extra files that were created file.remove("data.txt" ) file.remove("CODAindex.txt" ) file.remove("CODAchain1.txt" ) file.remove("CODAchain2.txt" ) file.remove("CODAchain3.txt" ) file.remove("inits1.txt" ) file.remove("inits2.txt" ) file.remove("inits3.txt" ) file.remove("model.txt" )
Here is the fitted spline curve to the number of unmarked fish available in each recovery sample
demo.fit$plots$fit.plot
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
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"),]
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
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)
There are several unusual features of the data set:
We already read in the data above. Here we set the rest of the parameters.
Don't forget to set the working directory as appropriate:
library("BTSPAS") # After which weeks do the YoY hatchery fish start to arrive. # Prior to this point, all YoY fish are wild and it is not # necessary to separate out the YoY wild vs hatchery demo2.hatch.after.YoY <- c(22) # julian weeks after which YoY hatchery fish arrive. # Which julian weeks have "bad" values. These will be set 0 (releases or recaptures) or missing (unmarked) and estimated. demo2.bad.m2 <- c() # list of julian weeks with bad m2 values demo2.bad.u2.A.YoY <- c() # list of julian weeks with bad u2.A.YoY values demo2.bad.u2.N.YoY <- c() # list of julian weeks with bad u2.N.YoY values demo2.bad.u2.A.1 <- c() # list of julian weeks with bad u2.A.YoY values demo2.bad.u2.N.1 <- c() # list of julian weeks with bad u2.N.YoY values # The clipping fraction for the current YoY and last year's YoY (which are now Age 1 fish) demo2.clip.frac.H.YoY <- .25 # what fraction of the YoY hatchery fish are adipose fin clipped demo2.clip.frac.H.1 <- .25 # what fraction of the Age1 hatchery fish are adipose fin clipped # The prefix for the output files: demo2.prefix <- "NF-2009-CH-WH-YoY-Age1" # Title for the analysis demo2.title <- "North Fork 2009 Chinook - Separation of YoY and Age 1 Wild and Hatchery Chinook" cat("*** Starting ",demo2.title, "\n\n") # 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)
# delete extra files that were created file.remove("data.txt" ) file.remove("CODAindex.txt" ) file.remove("CODAchain1.txt" ) file.remove("CODAchain2.txt" ) file.remove("CODAchain3.txt" ) file.remove("inits1.txt" ) file.remove("inits2.txt" ) file.remove("inits3.txt" ) file.remove("model.txt" )
Here is the fitted spline curve to the number of unmarked fish available of both stocks and ages.
demo2.fit$plots$fit.plot
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
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)
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
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)
There are several unusual features of the data set:
We already read in the data above. Here we set the rest of the parameters.
Don't forget to set the working directory as appropriate:
library("BTSPAS") # After which weeks do the hatchery fish start to arrive. Prior to this point, all fish are wild and it is not # necessary to separate out the wild vs hatchery demo3.hatch.after <- c(11) # julian weeks after which hatchery fish arrive. # Which julian weeks have "bad" values. These will be set to 0 (releases and recapture) or missing (unmarked captured) and estimated. demo3.bad.m2 <- c() # list of julian weeks with bad m2 values demo3.bad.u2.W.YoY <- c() # list of julian weeks with bad u2.W.YoY values demo3.bad.u2.W.1 <- c() # list of julian weeks with bad u2.W.1 values demo3.bad.u2.H.1 <- c() # list of julian weeks with bad u2.H.1 values # The prefix for the output files: demo3.prefix <- "demo-JC-2003-ST-TSPDE-WH" # Title for the analysis demo3.title <- "Junction City 2003 Steelhead - Separation of Wild and Hatchery YoY and Age 1+ Steelhead" cat("*** Starting ",demo3.title, "\n\n") # 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)
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.
# delete extra files that were created file.remove("data.txt" ) file.remove("CODAindex.txt" ) file.remove("CODAchain1.txt" ) file.remove("CODAchain2.txt" ) file.remove("CODAchain3.txt" ) file.remove("inits1.txt" ) file.remove("inits2.txt" ) file.remove("inits3.txt" ) file.remove("model.txt" )
Here is the fitted spline curve to the number of unmarked fish available of both stocks and ages.
demo3.fit$plots$fit.plot
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
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)),]
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.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.