4 Example of basic BTSPAS fit.
4.1 Reading in the data
Because the matrix is diagonal, and because the \(u2\) vector is the same length as the \(n1\) and \(m2\) vectors, the data can be entered as several columns. Here is an example of some raw data:
demo.data.csv <- textConnection(
'jweek, n1, m2, u2
9 ,0, 0, 4135
10,1465, 51, 10452
11,1106,121, 2199
12, 229, 25, 655
13, 20, 0, 308
14, 177, 17, 719
15, 702, 74, 973
16, 633, 94, 972
17,1370, 62, 2386
18, 283, 10, 469
19, 647, 32, 897
20, 276, 11, 426
21, 277, 13, 407
22, 333, 15, 526
23,3981,242, 39969
24,3988, 55, 17580
25,2889,115, 7928
26,3119,198, 6918
27,2478, 80, 3578
28,1292, 71, 1713
29,2326,153, 4212
30,2528,156, 5037
31,2338,275, 3315
32,1012,101, 1300
33, 729, 66, 989
34, 333, 44, 444
35, 269, 33, 339
36, 77, 7, 107
37, 62, 9, 79
38, 26, 3, 41
39, 20, 1, 23
40,4757,188, 35118
41,2876, 8, 34534
42,3989, 81, 14960
43,1755, 27, 3643
44,1527, 30, 1811
45, 485, 14, 679
46, 115, 4, 154')
demo.data <- read.csv(demo.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)
print(demo.data)
#> jweek n1 m2 u2
#> 1 9 0 0 4135
#> 2 10 1465 51 10452
#> 3 11 1106 121 2199
#> 4 12 229 25 655
#> 5 13 20 0 308
#> 6 14 177 17 719
#> 7 15 702 74 973
#> 8 16 633 94 972
#> 9 17 1370 62 2386
#> 10 18 283 10 469
#> 11 19 647 32 897
#> 12 20 276 11 426
#> 13 21 277 13 407
#> 14 22 333 15 526
#> 15 23 3981 242 39969
#> 16 24 3988 55 17580
#> 17 25 2889 115 7928
#> 18 26 3119 198 6918
#> 19 27 2478 80 3578
#> 20 28 1292 71 1713
#> 21 29 2326 153 4212
#> 22 30 2528 156 5037
#> 23 31 2338 275 3315
#> 24 32 1012 101 1300
#> 25 33 729 66 989
#> 26 34 333 44 444
#> 27 35 269 33 339
#> 28 36 77 7 107
#> 29 37 62 9 79
#> 30 38 26 3 41
#> 31 39 20 1 23
#> 32 40 4757 188 35118
#> 33 41 2876 8 34534
#> 34 42 3989 81 14960
#> 35 43 1755 27 3643
#> 36 44 1527 30 1811
#> 37 45 485 14 679
#> 38 46 115 4 154The first stratum was defined as julian week 9. At this time, 0 fish were captured, tagged, and released, but 4135 unmarked fish were recovered in the first recovery stratum. In the next week, 1465 fish were captured, tagged, and released, with 51 fish recaptured, and 10452 unmarked fish recovered.
4.2 Preliminary screening of the data
A pooled-Petersen estimator would add all of the marked, recaptured and unmarked fish to give an estimate of 10,602,698,039 which seems unbelievable!
What went wrong? Let us first examine a plot of the estimated capture efficiency at the second trap for each (recovery) julian week.
Empirical capture probabilities
There are several unusual features
- 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 at least two 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 appears to be heterogeneity in the capture probabilities.
- In some julian weeks, the number of marked fish released and recaptured is very small which lead to estimates with poor precision.
Similarly, let us look at the pattern of unmarked fish captured at the second trap:
Observed number of unmarked recaptures
There are two julian weeks where the number of unmarked fish captured suddenly jumps by several orders of magnitude (remember the above plot is on the log() scale). These jumps correspond to releases of hatchery fish into the system.
Finally, let us look at the individual estimates for each stratum found by computing a Petersen estimator for the total number of unmarked fish for each individual stratum:
Estimated log(total unmarked) by julian week
- The sudden jumps in abundance due to the hatchery releases is apparent, but the estimate in julian week 41 is far too large(!) with an estimate of almost 13 million fish from the simple stratified Petersen.
- There is a fairly regular pattern in abundance with a slow increase until the first hatchery release, followed by a steady decline, followed by a jump for the second hatchery release, followed by a steady decline.
4.3 Fitting the basic BTSPAS diagonal model
The \(BTSPAS\) package attempts to strike a balance between the completely pooled Petersen estimator and the completely stratified Petersen estimator. In the former, capture probabilities are assumed to equal for all fish in all strata, while in the latter, capture probabilities are allowed to vary among strata in no structured way. Furthermore, fish populations often have a general structure to the run, rather than arbitrarily jumping around from stratum to stratum.
Bonner and Schwarz (2011) developed a suite of models that add structure. In the basic model,
- A spline is used to smooth the total number of unmarked fish presenting themselves at the second trap over the strata
- A hierarchical model for the capture-probabilities is assumed where individual stratum capture probabilities are assumed to vary around a common mean.
The model also allows the user to use covariates to explain some of the variation in the capture probabilities. Bonner and Schwarz (2011) also developed models where fish are recaptured in more than stratum.
The \(BTSPAS\) package also has additional features and options:
- if \(u2\) is missing for any stratum, the program will use the spline to interpolate the number of unmarked fish in the population for the missing stratum.
- if \(n1\) and \(m2\) are 0, then these strata provide no information towards recapture probabilities. This is useful when no release take place in a stratum (e.g. trap did not run) and so you need ‘dummy’ values as placeholders. Of course is \(n1>0\) and \(m2=0\), this provides information that the capture probability may be small. If \(n1=0\) and \(m2>0\), this is an error (recoveries from no releases).
- the program allows you specify break points in the underlying spline to account for external events. We was in the above example, that hatchery fish were released at in julian weeks 23 and 40 resulting in sudden jump in abundance. The \(jump.after\) parameter gives the julian weeks just BEFORE the sudden jump, i.e. the spline is allowed to jump AFTER the julian weeks in jump.after.
- sometimes bad thing happen. The vector \(bad.m2\) indicates which julian weeks something went wrong. In the above example, the number of recoveries in julian week 41 is far below expectations and leads to an impossible Petersen estimate for julian week 41. Similarly, the vector \(bad.u2\) indicates in which julian weeks, the number of unmarked fish is suspect. In both cases, the suspect values of \(m2\) and \(u2\) are set to missing. Alternatively, the user can set the \(m2\) and \(u2\) values to missing in the data input directly. I arbitrarily chose the third julian week to demonstrate this feature.
I find it easiest if bad recaptures (a value of \(m2\)) result in zeroing out both \(n1\) and \(m2\) for that stratum.
The \(BTSPAS\) function also allows you specify
- The prefix is used to identify the output files for this run.
- The title is used to title the output.
- Various parameters to control the Bayesian MCMC phase of model fitting. Please contact us for help in setting these if problem arise.
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 is the spline allowed to jump?
demo.jump.after <- c(22,39) # julian weeks after which jump occurs
# Which julian weeks have "bad" recapture values. These will be set to missing and estimated.
demo.bad.m2 <- c(41) # list julian weeks with bad m2 values. This is used in the Trinity Example
demo.bad.u2 <- c(11) # list julian weeks with bad u2 values. [This was arbitrary to demostrate the feature.]
demo.bad.n1 <- c(38) # list julian weeks with bad n1 values. [This was arbitrary to demonstrate the feature.]
# The prefix for the output files:
demo.prefix <- "demo-JC-2003-CH-TSPDE"
# Title for the analysis
demo.title <- "Junction City 2003 Chinook "
cat("*** Starting ",demo.title, "\n\n")
#> *** Starting Junction City 2003 Chinook
# Make the call to fit the model and generate the output files
demo.fit <- TimeStratPetersenDiagError_fit(
title=demo.title,
prefix=demo.prefix,
time=demo.data$jweek,
n1=demo.data$n1,
m2=demo.data$m2,
u2=demo.data$u2,
jump.after=demo.jump.after,
bad.n1=demo.bad.n1,
bad.m2=demo.bad.m2,
bad.u2=demo.bad.u2,
InitialSeed=890110,
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: 890110
#> Random number seed for chain 127579
#> Random number seed for chain 693284
#> Random number seed for chain 289944
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 72
#> Unobserved stochastic nodes: 104
#> Total graph size: 1665
#>
#> 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.
The output object contains all of the results and can be saved for later interrogations. This is useful if the run takes considerable time (e.g. overnight) and you want to save the results for later processing. Notice that I didn’t save the results below as part of this vignette.
# save the results in a data dump that can be read in later using the load() command.
#save(list=c("demo.fit"), file="demo-fit-saved.Rdata") # save the results from this run4.4 The output from the basic fit
The final object has many components
names(demo.fit)#> [1] "n.chains" "n.iter" "n.burnin"
#> [4] "n.thin" "n.keep" "n.sims"
#> [7] "sims.array" "sims.list" "sims.matrix"
#> [10] "summary" "mean" "sd"
#> [13] "median" "root.short" "long.short"
#> [16] "dimension.short" "indexes.short" "last.values"
#> [19] "program" "model.file" "isDIC"
#> [22] "DICbyR" "pD" "DIC"
#> [25] "model" "parameters.to.save" "plots"
#> [28] "runTime" "report" "data"The plots sub-object contains many plots:
names(demo.fit$plots)#> [1] "init.plot" "fit.plot" "logitP.plot"
#> [4] "acf.Utot.plot" "post.UNtot.plot" "gof.plot"
#> [7] "trace.logitP.plot" "trace.logU.plot"In particular, it contains plots of the initial spline fit (init.plot), the final fitted spline (fit.plot), the estimated capture probabilities (on the logit scale) (logitP.plot), plots of the distribution of the posterior sample for the total unmarked and marked fish (post.UNtot.plot) and model diagnostic plots (goodness of fit (gof.plot), trace (trace…plot), and autocorrelation plots (act.Utot.plot).
These plots are all created using the \(ggplot2\) packages, so the user can modify the plot (e.g. change titles etc).
The \(BTSPAS\) program also creates a report, which includes information about the data used in the fitting, the pooled- and stratified-Petersen estimates, a test for pooling, and summaries of the posterior. Only the first few lines are shown below:
head(demo.fit$report)
#> [1] "Time Stratified Petersen with Diagonal recaptures and error in smoothed U - Sun Oct 24 15:37:06 2021"
#> [2] "Version: 2021-11-02 "
#> [3] ""
#> [4] ""
#> [5] ""
#> [6] " Junction City 2003 Chinook Results "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 jump in the spline when hatchery fish are released is evident. The actual number of unmarked fish is allowed to vary around the spline as shown below.
The distribution of the posterior sample for the total number unmarked and total abundance is available:
demo.fit$plots$post.UNtot.plot
Distribution of posterior samples
A plot of the \(logit(P)\) is
demo.fit$plots$logitP.plot
Estimates of logit(p)
In cases where there is no information, \(BTSPAS\) has interpolated based on the distribution of catchability in the other strata and so the credible interval is very wide (e.g. julian weeks 7, 13, and 41).
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 and total abundance:
demo.fit$summary[ row.names(demo.fit$summary) %in% c("Ntot","Utot"),]
#> mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
#> Ntot 5834690 520650.8 5091090 5498342 5765732 6066297 6949601 1.004049 500
#> Utot 5784225 520650.8 5040625 5447877 5715267 6015832 6899136 1.004053 500This also includes the Rubin-Brooks-Gelman statistic (\(Rhat\)) on mixing of the chains and the effective sample size of the posterior (after accounting for autocorrelation).
The estimated total abundance from \(BTSPAS\) is 5,834,690 (SD 520,651 ) fish.
Samples from the posterior are also included in the sims.matrix, sims.array and sims.list elements of the results object.
It is always important to do model assessment before accepting the results from the model fit. Please contact me for details on how to interpret the goodness of fit, trace, and autocorrelation plots.
