Diagonal Case

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

The two-sample capture-recapture experiment is one of the simplest possible studies with a long history. The standard Lincoln-Petersen estimator used to estimate abundance is \[ \widehat{N} = \frac{n_1 n_2}{m_2}\] where \(n_1\) is the number of animals captured, marked and released at the first capture event; \(n_2\) is the number of animals captured at the second capture event; and \(m_2\) is the number of animals from \(n_1\) that were recaptured at the second event (i.e. the number of recaptured marked animals).

A key assumption of the Lincoln-Petersen estimator is that there is no correlation in capture-probabilities at the two events. The most common way in which this can occur is if the probability of capture is equal for all animals at either the first or second event.

If capture probabilities are heterogeneous, then estimates can be biased. One way to account for heterogeneous capture-probabilities is through stratification. For example, if males and females have different catchabilities, then separate Lincoln-Petersen estimators can be computed for each sex, and the estimated abundance for the entire population is found by summing the two estimated abundances (one for each sex).

Stratification can be based on animal attributes (e.g. sex), geographic location (e.g. parts of a study area), or temporal (e.g. when captured). The \(BTSPAS\) package deals with temporal stratification.

3 Experimental Protocol

Consider an experiment to estimate the number of outgoing smolts on a small river. The run of smolts extends over several weeks. As smolts migrate, they are captured and marked with individually numbered tags and released at the first capture location using, for example, a fishwheel. The migration continues, and a second fishwheel takes a second sample several kilometers down stream. At the second fishwheel, the captures consist of a mixture of marked (from the first fishwheel) and unmarked fish.

The efficiency of the fishwheels varies over time in response to stream flow, run size passing the wheel and other uncontrollable events. So it is unlikely that the capture probabilities are equal over time at either location, i.e. are heterogeneous over time.

We suppose that we can temporally stratify the data into, for example, weeks, where the capture-probabilities are (mostly) homogeneous at each wheel in each week. Furthermore, suppose that fish captured and marked in each week tend to migrate together so that they are captured in a single subsequent stratum. For example, suppose that in each julian week \(j\), \(n1[j]\) fish are marked and released above the rotary screw trap. Of these, \(m2[j]\) are recaptured. All recaptures take place in the week of release, i.e. the matrix of releases and recoveries is diagonal. The \(n1[j]\) and \(m2[j]\) establish the capture efficiency of the second trap in julian week \(j\).

At the same time, \(u2[j]\) unmarked fish are captured at the screw trap.

This implies that the data can be structured as a diagonal array similar to:

Recovery Stratum
               tagged    rs1   rs2    rs3 ...  rsk  
Marking   ms1    n1[1]  m2[1]    0      0  ...   0  
Stratum   ms2    n1[2]   0    m2[2]     0  ...   0  
          ms3    n1[3]   0       0  m2[3]  ...   0  
          ...  
          msk    n1[k]   0       0      0  ... m2[k]  
Newly  
Untagged               u2[1]  u2[2] u2[3]  ... u2[k]  
captured  

Here the tagging and recapture events have been stratified into \(k\) temporal strata. Marked fish from one stratum tend to move at similar rates and so are recaptured together with unmarked fish. Recaptures of marked fish take place along the “diagonal.”

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   154

The 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

Empirical capture probabilities

There are several unusual features

Similarly, let us look at the pattern of unmarked fish captured at the second trap:

Observed number of unmarked recaptures

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

Estimated log(total unmarked) by julian week

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,

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:

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

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 run

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

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

Distribution of posterior samples

A plot of the \(logit(P)\) is

demo.fit$plots$logitP.plot
Estimates of logit(p)

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   500

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

5 Fixing p’s

In some cases, the second trap is not running and so there are no recaptures of tagged fish and no captures of untagged fish.

We need to set the p’s in these strata to 0 rather than letting BTSPAS impute a value.

5.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
 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,  0,     0
38,  26,  0,     0
39,  20,  0,     0
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,  0,    0')

demo2.data <- read.csv(demo2.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo2.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   0     0
#> 30    38   26   0     0
#> 31    39   20   0     0
#> 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   0     0

Notice that there are no recaptures of marked fish and no recaptures of unmarked fish in julian weeks 37, 38, 39, and 46 when the trap was not operating. Notice that this differs from the case where a small number of tagged fish released with no recaptures but captures of tagged fish occur such as in julian week 13.

5.2 Fitting the BTSPAS diagonal model and fixing p. 

Two additional arguments to the BTSPAS allow you specify the julian weeks in which the capture probability is fixed to a known (typically zero) value.

demo2.logitP.fixed <- c(37,38,39, 46)
demo2.logitP.fixed.values <- rep(-10, length(demo2.logitP.fixed))

The strata where the value of p is to be fixed is specified along with the value (on the logit scale) at which the capture probabilities are fixed. Technically, the \(logit(0.0)\) is negative infinity, but \(-10\) is ``close enough’’.

The rest of the call is basically the same – don’t forget to specify the additional arguments in the call

library("BTSPAS")  

# After which weeks is the spline allowed to jump?
demo2.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.
demo2.bad.m2     <- c(41)   # list julian weeks with bad m2 values. This is used in the Trinity Example
demo2.bad.u2     <- c(11)   # list julian weeks with bad u2 values. [This was arbitrary to demostrate the feature.]
demo2.bad.n1     <- c(38)   # list julian weeks with bad n1 values. [This was arbitrary to demonstrate the feature.]

# The prefix for the output files:
demo2.prefix <- "demo2-JC-2003-CH-TSPDE" 

# Title for the analysis
demo2.title <- "Junction City 2003 Chinook with p fixed "

cat("*** Starting ",demo2.title, "\n\n")
#> *** Starting  Junction City 2003 Chinook with p fixed

# Make the call to fit the model and generate the output files
demo2.fit <- TimeStratPetersenDiagError_fit(
                  title=demo2.title,
                  prefix=demo2.prefix,
                  time=demo2.data$jweek,
                  n1=demo2.data$n1, 
                  m2=demo2.data$m2, 
                  u2=demo2.data$u2,
                  jump.after=demo2.jump.after,
                  logitP.fixed=demo2.logitP.fixed,               # ***** NEW ****8
                  logitP.fixed.values=demo2.logitP.fixed.values, # ***** NEW *****
                  bad.n1=demo2.bad.n1,
                  bad.m2=demo2.bad.m2,
                  bad.u2=demo2.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: 100
#>    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 ***

5.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available in each recovery sample. Note how the spline interpolates across the julian weeks when no unmarked fish were captured in julian weeks 37, 38, 39, and 46 but the uncertainty is much larger.

demo2.fit$plots$fit.plot
Fitted spline curve with fixed p's

Fitted spline curve with fixed p’s

The jump in the spline when hatchery fish are released is evident.

The distribution of the posterior sample for the total number unmarked and total abundance is available as before:

demo2.fit$plots$post.UNtot.plot
Distribution of posterior samples

Distribution of posterior samples

A plot of the \(logit(P)\) is

demo2.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

Notice how the fixed values at \(-10\) (on the logit scale) occur.

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:

demo2.fit$summary[ row.names(demo2.fit$summary) %in% c("Ntot","Utot"),]
#>         mean       sd    2.5%     25%     50%     75%   97.5%     Rhat n.eff
#> Ntot 5851224 506714.6 5140894 5516261 5769217 6062302 7073362 1.007011   300
#> Utot 5800759 506714.6 5090429 5465796 5718752 6011837 7022897 1.007020   300

The estimated total abundance from \(BTSPAS\) is 5,851,224 (SD 506,715 ) fish.

6 Using covariates to model the p’s

BTSPAS also allows you to model the p’s with additional covariates, such a temperature, stream flow, etc. It is not possible to use covariates to model the total number of unmarked fish.

6.1 Reading in the data

Here is an example of some raw data that includes the covariate \(log(flow)\):

demo3.data.csv <- textConnection(
'jweek, n1, m2,    u2, logflow
 9,    0,    0,  4135, 6.617212
10, 1465,   51, 10452, 6.51217
11, 1106,  121,  2199, 7.193686
12,  229,   25,   655, 6.960754
13,   20,    0,   308, 7.008376
14,  177,   17,   719, 6.761573
15,  702,   74,   973, 6.905753
16,  633,   94,   972, 7.062314
17, 1370,   62,  2386, 7.600188
18,  283,   10,   469, 8.246509
19,  647,   32,   897, 8.110298
20,  276,   11,   426, 8.035001
21,  277,   13,   407, 7.859965
22,  333,   15,   526, 7.774255
23, 3981,  242, 39969, 7.709116
24, 3988,   55, 17580, 7.653766
25, 2889,  115,  7928, 7.622105
26, 3119,  198,  6918, 7.593734
27, 2478,   80,  3578, 7.585063
28, 1292,   71,  1713, 7.291072
29, 2326,  153,  4212, 6.55556
30, 2528,  156,  5037, 6.227665
31, 2338,  275,  3315, 6.278789
32, 1012,  101,  1300, 6.273685
33,  729,   66,   989, 6.241111
34,  333,   44,   444, 6.687999
35,  269,   33,   339, 7.222566
36,   77,    7,   107, 7.097194
37,   62,    9,    79, 6.949993
38,   26,    3,    41, 6.168714
39,   20,    1,    23, 6.113682
40, 4757,  188, 35118, 6.126557
41, 2876,    8, 34534, 6.167217
42, 3989,   81, 14960, 5.862413
43, 1755,   27,  3643, 5.696614
44, 1527,   30,  1811, 5.763847
45,  485,   14,   679, 5.987528
46,  115,    4,   154, 5.912344')

demo3.data <- read.csv(demo3.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo3.data)
#>    jweek   n1  m2    u2  logflow
#> 1      9    0   0  4135 6.617212
#> 2     10 1465  51 10452 6.512170
#> 3     11 1106 121  2199 7.193686
#> 4     12  229  25   655 6.960754
#> 5     13   20   0   308 7.008376
#> 6     14  177  17   719 6.761573
#> 7     15  702  74   973 6.905753
#> 8     16  633  94   972 7.062314
#> 9     17 1370  62  2386 7.600188
#> 10    18  283  10   469 8.246509
#> 11    19  647  32   897 8.110298
#> 12    20  276  11   426 8.035001
#> 13    21  277  13   407 7.859965
#> 14    22  333  15   526 7.774255
#> 15    23 3981 242 39969 7.709116
#> 16    24 3988  55 17580 7.653766
#> 17    25 2889 115  7928 7.622105
#> 18    26 3119 198  6918 7.593734
#> 19    27 2478  80  3578 7.585063
#> 20    28 1292  71  1713 7.291072
#> 21    29 2326 153  4212 6.555560
#> 22    30 2528 156  5037 6.227665
#> 23    31 2338 275  3315 6.278789
#> 24    32 1012 101  1300 6.273685
#> 25    33  729  66   989 6.241111
#> 26    34  333  44   444 6.687999
#> 27    35  269  33   339 7.222566
#> 28    36   77   7   107 7.097194
#> 29    37   62   9    79 6.949993
#> 30    38   26   3    41 6.168714
#> 31    39   20   1    23 6.113682
#> 32    40 4757 188 35118 6.126557
#> 33    41 2876   8 34534 6.167217
#> 34    42 3989  81 14960 5.862413
#> 35    43 1755  27  3643 5.696614
#> 36    44 1527  30  1811 5.763847
#> 37    45  485  14   679 5.987528
#> 38    46  115   4   154 5.912344

A preliminary plot of the empirical logit (excluding those weeks when the trap was not running) shows an approximate quadratic fit to \(log(flow)\), but the uncertainty in each week is enormous!

6.2 Fitting the BTSPAS diagonal model and fixing p and covariates for p. 

We need to create a matrix with the covariate values. We will need three columns - one for the intercept, the value of log(flow) and the square of its values. In practice, it is often advisable to standardize covariates to prevent numerical difficulties, but in this case, the values are small enough that standardization is not really needed.

demo3.logitP.cov <- cbind(1, demo3.data$logflow, demo3.data$logflow^2)
head(demo3.logitP.cov)
#>      [,1]     [,2]     [,3]
#> [1,]    1 6.617212 43.78749
#> [2,]    1 6.512170 42.40836
#> [3,]    1 7.193686 51.74912
#> [4,]    1 6.960754 48.45210
#> [5,]    1 7.008376 49.11733
#> [6,]    1 6.761573 45.71887

The rest of the call is basically the same – don’t forget to specify the additional arguments in the call

library("BTSPAS")  

# After which weeks is the spline allowed to jump?
demo3.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.
demo3.bad.m2     <- c(41)   # list julian weeks with bad m2 values. This is used in the Trinity Example
demo3.bad.u2     <- c(11)   # list julian weeks with bad u2 values. [This was arbitrary to demostrate the feature.]
demo3.bad.n1     <- c(38)   # list julian weeks with bad n1 values. [This was arbitrary to demonstrate the feature.]

# The prefix for the output files:
demo3.prefix <- "demo3-JC-2003-CH-TSPDE" 

# Title for the analysis
demo3.title <- "Junction City 2003 Chinook with covariates for p "

cat("*** Starting ",demo3.title, "\n\n")
#> *** Starting  Junction City 2003 Chinook with covariates for p

# Make the call to fit the model and generate the output files
demo3.fit <- TimeStratPetersenDiagError_fit(
                  title=demo3.title,
                  prefix=demo3.prefix,
                  time=demo3.data$jweek,
                  n1=demo3.data$n1, 
                  m2=demo3.data$m2, 
                  u2=demo3.data$u2,
                  jump.after=demo3.jump.after,
                  logitP.cov = demo3.logitP.cov,                 # ***** NEW *****
                  bad.n1=demo3.bad.n1,
                  bad.m2=demo3.bad.m2,
                  bad.u2=demo3.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: 106
#>    Total graph size: 1825
#> 
#> 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 ***

6.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available in each recovery sample.

demo3.fit$plots$fit.plot
Fitted spline curve with covariate for p

Fitted spline curve with covariate for p

The jump in the spline when hatchery fish are released is evident.

The distribution of the posterior sample for the total number unmarked and total abundance is available as before:

demo3.fit$plots$post.UNtot.plot
Distribution of posterior samples

Distribution of posterior samples

A plot of the \(logit(P)\) is

demo3.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

Here is a plot of the estimated \(logit(p)\)’s vs. the log(flow) and its fitted curve:

There is virtually no evidence of a relationship with flow because of the very large uncertainties in each of the estimated \(logit(p)\).

The estimated coefficients of the quadratic relationship between logit(p) and log(flow) are:

round(demo3.fit$summary[demo3.coeff.row.index,],3)
#>                  mean    sd   2.5%    25%    50%    75%  97.5%  Rhat n.eff
#> beta.logitP[1] -3.294 0.948 -5.170 -3.927 -3.284 -2.663 -1.407 1.003  1500
#> beta.logitP[2]  0.115 0.328 -0.541 -0.101  0.111  0.333  0.743 1.001  1500
#> beta.logitP[3] -0.007 0.031 -0.068 -0.028 -0.007  0.012  0.054 1.001  1500

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:

demo3.fit$summary[ row.names(demo3.fit$summary) %in% c("Ntot","Utot"),]
#>         mean       sd    2.5%     25%     50%     75%   97.5%     Rhat n.eff
#> Ntot 5849700 496356.9 5126187 5519334 5762912 6082048 7141051 1.002868   730
#> Utot 5799235 496356.9 5075722 5468869 5712447 6031583 7090586 1.002868   730

The estimated total abundance from \(BTSPAS\) is 5,849,700 (SD 496,357 ) fish.

7 Using covariates to model the p’s - prior information about beta coefficients

BTSPAS also allows you to model the p’s with additional covariates, such a temperature, stream flow, etc. It is not possible to use covariates to model the total number of unmarked fish.

In some cases, you may have additional information about the effect of the covariates that you would like to incorporate into the analysis. For example, a rotary screw trap may have run for many years and plots of the relationship between the logit(catchability) vs. log(flow) generally shows a relationship that you may not be able to capture in a single year because of lack of contrast (i.e. the flow within a year does not vary enough) or because of smallish sample sizes.

BTSPAS allows you to specify prior information on the coefficients of the relationship between logit(catchability) and covariates.

We will show an example where four models are fit to the same data

7.1 Reading in the data

Here is an example of some (fictitious) raw data that includes the covariate \(log(flow)\):

demo5.data.csv <- textConnection(
'time,   n1,     m2,       u2,       logFlow
1   ,      0   ,    0   ,      21   ,   5.415
2   ,     56   ,    8   ,    2266   ,   5.358
3   ,   1009   ,   59   ,   11314   ,   6.737
4   ,   1284   ,   25   ,    5035   ,   7.993
5   ,   1504   ,   13   ,     396   ,   8.693
6   ,      0   ,    0   ,      45   ,   8.861
7   ,      0   ,    0   ,      26   ,   8.587
8   ,   1560   ,   17   ,      12   ,   8.347
9   ,   1643   ,   14   ,      43   ,   8.260
10   ,     0   ,    0   ,      63   ,   8.606
11   ,  1487   ,    7   ,      24   ,   8.671
12   ,     0   ,    0   ,       5   ,   8.737
13   ,     0   ,    0   ,       4   ,   7.862')

demo5.data <- read.csv(demo5.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo5.data)
#>    time   n1 m2    u2 logFlow
#> 1     1    0  0    21   5.415
#> 2     2   56  8  2266   5.358
#> 3     3 1009 59 11314   6.737
#> 4     4 1284 25  5035   7.993
#> 5     5 1504 13   396   8.693
#> 6     6    0  0    45   8.861
#> 7     7    0  0    26   8.587
#> 8     8 1560 17    12   8.347
#> 9     9 1643 14    43   8.260
#> 10   10    0  0    63   8.606
#> 11   11 1487  7    24   8.671
#> 12   12    0  0     5   8.737
#> 13   13    0  0     4   7.862

A preliminary plot of the empirical logit (excluding those weeks when the trap was not running) shows an approximate linear fit to \(log(flow)\), but the uncertainty in each week is enormous!

7.2 Fitting the BTSPAS diagonal model with no covariates

We start with fitting BTSPAS with no covariates

fit1 <- BTSPAS::TimeStratPetersenDiagError_fit(
          title="no covariates"
         ,prefix="fit1"
         ,time  = demo5.data$time
         ,n1    = demo5.data$n1
         ,m2    = demo5.data$m2
         ,u2    = demo5.data$u2
         ,InitialSeed=234323
         ,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: 234323 
#> Random number seed for chain  991738 
#> Random number seed for chain  937097 
#> Random number seed for chain  622697 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 43
#>    Total graph size: 427
#> 
#> 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 ***

Estimated values of the overall mean logitP and the residual variation around the fitted line are:

select <-  grepl("beta.logitP", row.names(fit1$summary)) 
fit1$summary[select,][1]
#> [1] -4.469161

select <-  grepl("sigmaP", row.names(fit1$summary)) # sd(logitP) around mean logitP or regression curve if covariates present
fit1$summary[select,]
#>         mean           sd         2.5%          25%          50%          75%        97.5%         Rhat        n.eff 
#>    1.3816123    0.4477664    0.7319510    1.0635651    1.3070595    1.6137544    2.4928057    1.0014468 2900.0000000

The plot of the estimated logit(catchability) vs. log(flow) and the overall mean is:

plotdata1 <- data.frame( logFlow = demo5.data$logFlow,
                         logitP  = fit1$mean$logitP,
                         time  = demo5.data$time,
                         n1    = demo5.data$n1)
plotdata1$fit <- 'fit1 - nocovariates'
ggplot(data=plotdata1, aes(x=logFlow, y=logitP))+
   ggtitle("no covariates")+
   geom_point(aes(size=n1))+
   #geom_smooth(method="lm", se=FALSE)+
   geom_text(aes(label=time), vjust=-.5)+
   geom_hline(yintercept=fit1$mean$beta.logitP[1], color="red")
Fit1: logitP vs log(flow)

Fit1: logitP vs log(flow)

The appears to be a relationship with log(flow) that has not been captured with this model.

7.3 Fitting the BTSPAS diagonal model with using log(flow) and default priors

We need to create a matrix with the covariate values. We will need two columns - one for the intercept, and one for the value of log(flow). In practice, it is often advisable to standardize covariates to prevent numerical difficulties, but in this case, the values are small enough that standardization is not really needed.


fit2 <- BTSPAS::TimeStratPetersenDiagError_fit(
          title="log(Flow)- default prior"
         ,prefix="fit2"
         ,time  = demo5.data$time
         ,n1    = demo5.data$n1
         ,m2    = demo5.data$m2
         ,u2    = demo5.data$u2
         ,logitP.cov=cbind( 1, demo5.data$logFlow)
         ,InitialSeed=3542343
         ,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: 3542343 
#> Random number seed for chain  226849 
#> Random number seed for chain  411513 
#> Random number seed for chain  660797 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 44
#>    Total graph size: 471
#> 
#> 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 ***

Estimated values of the beta coefficients (the intercept and slope) and the residual variation around the fitted line are:

select <-  grepl("beta.logitP", row.names(fit2$summary)) 
fit2$summary[select,][1:2]
#> [1]  0.7610026 -0.6357489

select <-  grepl("sigmaP", row.names(fit2$summary)) # sd(logitP) around mean logitP or regression curve if covariates present
fit2$summary[select,]
#>         mean           sd         2.5%          25%          50%          75%        97.5%         Rhat        n.eff 
#>   0.57240707   0.47283283   0.04603051   0.22107098   0.43432392   0.79748650   1.76785558   1.00329914 790.00000000

The plot of the estimated logit(catchability) vs. log(flow) and the overall mean is:

plotdata2 <- data.frame(logFlow = demo5.data$logFlow,
                        logitP  = fit2$mean$logitP,
                        time  = demo5.data$time,
                        n1    = demo5.data$n1)
plotdata2$fit <- 'fit2  - log flow default priors'
ggplot(data=plotdata2, aes(x=logFlow, y=logitP))+
   ggtitle("log(flow) - default priors")+
   geom_point(aes(size=n1))+
   #geom_smooth(method="lm", se=FALSE)+
   geom_text(aes(label=time), vjust=-.5)+
   #geom_hline(yintercept=fit1$mean$beta.logitP[1], color="red")+
   geom_abline(intercept=fit2$mean$beta.logitP[1],
               slope    =fit2$mean$beta.logitP[2], color="green")
Fit2: logitP vs log(flow)

Fit2: logitP vs log(flow)

We now have a relationship with log(flow), but as you will see later, the evidence is not very strong.

7.4 Fitting the BTSPAS diagonal model with using log(flow) and weak prior

Prior information on the beta coefficients (the intercept and slope) are given using the prior.beta.logitP.mean and prior.beta.logitP.sd parameters in the call. The first specifies the values of the intercept and slope and the second specifies the uncertainty in these prior values.

Consider the fit:


fit3 <- BTSPAS::TimeStratPetersenDiagError_fit(
          title="log(Flow) - weak priors"
         ,prefix="fit3"
         ,time  = demo5.data$time
         ,n1    = demo5.data$n1
         ,m2    = demo5.data$m2
         ,u2    = demo5.data$u2
         ,logitP.cov=cbind( 1, demo5.data$logFlow)
         ,prior.beta.logitP.mean=c( .3, -.7)  # prior for intercept and slope on logit scale - mean
         ,prior.beta.logitP.sd  =c(  2,   2)  # prior for intercept and slope on logit scale - sd
         ,InitialSeed=3542343
         ,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: 3542343 
#> Random number seed for chain  226849 
#> Random number seed for chain  411513 
#> Random number seed for chain  660797 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 44
#>    Total graph size: 469
#> 
#> 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 ***

Here the prior for the intercept is set to 0.3 and the prior for the slope is set to -.7 (the prior.beta.logitP.mean).

The prior.beta.logitP.sd gives the uncertainty (like a standard error) in these values. In this fit, I used a large uncertainty, a standard deviation of 2 for both. This indicates the prior puts information on \(.3 \pm 2 \times 2\) for the intercept and \(-7 \pm 2 \times 2\) for the slope, i.e. about 95% of the information in the prior within two standard deviations of the mean.

Estimated values of the beta coefficients (the intercept and slope) and the residual variation around the fitted line are:

select <-  grepl("beta.logitP", row.names(fit3$summary)) 
fit3$summary[select,][1:2]
#> [1]  3.0685003 -0.9115914

select <-  grepl("sigmaP", row.names(fit3$summary)) # sd(logitP) around mean logitP or regression curve if covariates present
fit3$summary[select,]
#>         mean           sd         2.5%          25%          50%          75%        97.5%         Rhat        n.eff 
#> 2.416068e-01 2.067102e-01 3.054742e-02 9.338818e-02 1.847661e-01 3.249542e-01 7.776406e-01 1.000772e+00 6.000000e+03

The plot of the estimated logit(catchability) vs. log(flow) and the overall mean is:

plotdata3 <- data.frame( logFlow = demo5.data$logFlow,
                        logitP  = fit3$mean$logitP,
                        time  = demo5.data$time,
                        n1    = demo5.data$n1)

plotdata3$fit <- 'fit3 - log(flow) - weak priors'
ggplot(data=plotdata3, aes(x=logFlow, y=logitP))+
   ggtitle("log(flow) - weak priors")+
   geom_point(aes(size=n1))+
   #geom_smooth(method="lm", se=FALSE)+
   geom_text(aes(label=time), vjust=-.5)+
   #geom_hline(yintercept=fit1$mean$beta.logitP[1], color="red")+
   #geom_abline(intercept=fit2$mean$beta.logitP[1],
   #            slope    =fit2$mean$beta.logitP[2], color="green")+
   geom_abline(intercept=fit3$mean$beta.logitP[1],
               slope    =fit3$mean$beta.logitP[2], color="brown", linetype="solid")
Fit3: logitP vs log(flow)

Fit3: logitP vs log(flow)

The weak information on the slope gives a boost to the relationship between log(flow) and the the logit(catchability) especially for those weeks when the sample size is very small.

7.5 Fitting the BTSPAS diagonal model with using log(flow) and strong prior

We repeat the fit but with very strong prior information:


fit4 <- BTSPAS::TimeStratPetersenDiagError_fit(
          title="log(Flow) - strong priors"
         ,prefix="fit4"
         ,time  = demo5.data$time
         ,n1    = demo5.data$n1
         ,m2    = demo5.data$m2
         ,u2    = demo5.data$u2
         ,logitP.cov=cbind( 1, demo5.data$logFlow)
         ,prior.beta.logitP.mean=c( .3, -.7)  # prior for intercept and slope on logit scale - mean
         ,prior.beta.logitP.sd  =c(.01, .01)  # prior for intercept and slope on logit scale - sd
         ,InitialSeed=3542343
         ,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: 3542343 
#> Random number seed for chain  226849 
#> Random number seed for chain  411513 
#> Random number seed for chain  660797 
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 44
#>    Total graph size: 469
#> 
#> 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 ***

Here the prior for the intercept is set to 0.3 and the prior for the slope is set to -.7 (the prior.beta.logitP.mean).

The prior.beta.logitP.sd gives the uncertainty (like a standard error) in these values. In this fit, I used a very small uncertainty, a standard deviation of .01 for both. This indicates the prior puts information on \(.3 \pm 2 \times .01\) for the intercept and \(-7 \pm 2 \times .01\) for the slope, i.e. about 95% of the information in the prior within two standard deviations of the mean.

Estimated values of the beta coefficients (the intercept and slope) and the residual variation around the fitted line are:

select <-  grepl("beta.logitP", row.names(fit4$summary)) 
fit4$summary[select,][1:2]
#> [1]  0.3002301 -0.6965861

select <-  grepl("sigmaP", row.names(fit4$summary)) # sd(logitP) around mean logitP or regression curve if covariates present
fit4$summary[select,]
#>         mean           sd         2.5%          25%          50%          75%        97.5%         Rhat        n.eff 
#>    1.2702315    0.4088080    0.7081187    0.9870669    1.1955885    1.4680574    2.2657949    1.0021155 1500.0000000

The plot of the estimated logit(catchability) vs. log(flow) and the overall mean is:

plotdata4 <- data.frame( logFlow = demo5.data$logFlow,
                         logitP  = fit4$mean$logitP,
                         time    = demo5.data$time,
                         n1      = demo5.data$n1)
plotdata4$fit <- 'fit4 - log(flow) - strong priors'
ggplot(data=plotdata4, aes(x=logFlow, y=logitP))+
   ggtitle("log(flow) - strong prior")+
   geom_point(aes(size=n1))+
   #geom_smooth(method="lm", se=FALSE)+
   geom_text(aes(label=time), vjust=-.5)+
   #geom_hline(yintercept=fit1$mean$beta.logitP[1], color="red")+
   #geom_abline(intercept=fit2$mean$beta.logitP[1],
   #            slope    =fit2$mean$beta.logitP[2], color="green")+
   #geom_abline(intercept=fit3$mean$beta.logitP[1],
   #            slope    =fit3$mean$beta.logitP[2], color="brown", linetype="solid")+
   geom_abline(intercept=fit4$mean$beta.logitP[1],
               slope    =fit4$mean$beta.logitP[2], color="brown", linetype="dashed")
Fit4: logitP vs log(flow)

Fit4: logitP vs log(flow)

Notice that the (strong) prior is now in conflict with the data. The model now believes that the variation around the line must be large, which allows it to move estimates of logit(catchability) below the line.

7.6 Comparing the results of the fits

7.6.1 Estimates of abundance

Comparing estimates of abundance among fit
Source Mean SD
No covariates 610227 103992
Default prior 601688 74420
Weak prior 612352 65280
Strong prior 636210 94309

There appears to be little impact on the estimate of abundance, but notice that the uncertainty declines as you add information from fit1 to *fit3(), but when you add a strong conflicting prior (see below), the uncertainty now increases.

7.6.2 Estimates of coefficients

Comparing estimates of beta coefficients among fit
Source Estimate Mean SD
No covariates Intercept -4.469 0.514
Default prior Intercept 0.761 2.490
Default prior Slope -0.636 0.297
Weak prior Intercept 3.069 1.130
Weak prior Slope -0.912 0.140
Strong prior Intercept 0.300 0.010
Strong prior Slope -0.697 0.010

With the strong prior, the data plays essentially no role in determining the slope and intercept. With a weak prior, the estimated slope has lower precision than with the even weaker (default) prior.

7.6.3 Comparing the residual standard deviation around the line of fit

Comparing estimates of residual variation in logitP among fit
Source Mean SD
No covariates 1.382 0.448
Default prior 0.572 0.473
Weak prior 0.242 0.207
Strong prior 1.270 0.409

The residual variation declines as more information is added via the prior for the first 3 fits, but then increases when a strong (conflicting) prior is added (last fit).

7.6.4 Comparing the fits

All fits: logitP vs log(flow)

All fits: logitP vs log(flow)

In the plots of logitP vs. log(flow), the red horizontal line is from fit1 (no covariates) and represents the mean. The green line is from the fit of logitP vs logFlow using the default prior. The solid brown line is from the fit of logitP vs logFlow for the weak prior and the dashed brown line is for the strong prior. The size of the dots represents n1 the number of fish released.

So for fit1, there is very little data for time 1 and time 2, and and so their logitP are basically free to vary. For fit2, the default prior gives some (but not much information) about the relationship between logitP and logFlow, so the estimates of logitP for periods 1 and 2 are moved “closer” to the green line. For fit3, the weak prior gives more information and so the points move very close to the line. As well, the periods with lots of data can be well fit with very little noise and so the model says the the noise of periods in logitP must also be small, and so all the points are forced to lie on the line. For fit4, the very strong prior is placed on a WRONG line (the dashed brown). Now the model is in a quandary and says that the only way the logitP for weeks 3, 4, etc. with lots of data are consistent with the dashed brown line is if the among period variance is very large and so the logitP for weeks with very poor data are allowed to move away from the brown dashed line.

We can also compare the spline fit to the number of unmarked fish:

Allfts: spline fit to U

Allfts: spline fit to U

The spline is only affected slightly.

We can also compare the trend of logitP over time:

All fits: logitP vs time

All fits: logitP vs time

The trend over time is mostly unchanged for period with lots of data, and for periods with very small amount of data, the points are allowed to vary as needed.

8 Adding structure to the p’s with a spline.

The previous example still showed some temporal structure in the \(p\)’s. This additional structure can be imposed by using a spline for the \(p\)’s.

8.1 Reading in the data

Here is an example of some raw data:

demo4.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')

demo4.data <- read.csv(demo4.data.csv, header=TRUE, as.is=TRUE, strip.white=TRUE)

print(demo4.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   154

A preliminary plot of the empirical logit (excluding those weeks when the trap was not running) shows some temporal structure:

8.2 Fitting the BTSPAS diagonal model and a spline for p. 

We can create a set of covariates that serve as the basis for the spline over time using the \(bs()\) function from the splines package:

library(splines)
demo4.logitP.cov <- bs(1:length(demo4.data$n1), df=floor(length(demo4.data$n1)/4), intercept=TRUE)
head(demo4.logitP.cov)
#>                1         2         3           4 5 6 7 8 9
#> [1,] 1.000000000 0.0000000 0.0000000 0.000000000 0 0 0 0 0
#> [2,] 0.588138906 0.3756145 0.0355359 0.000710718 0 0 0 0 0
#> [3,] 0.308471364 0.5593351 0.1265078 0.005685744 0 0 0 0 0
#> [4,] 0.135411525 0.5959371 0.2494620 0.019189387 0 0 0 0 0
#> [5,] 0.043373542 0.5301956 0.3809449 0.045485953 0 0 0 0 0
#> [6,] 0.006771563 0.4068861 0.4975026 0.088839753 0 0 0 0 0

The rest of the call is basically the same – don’t forget to specify the additional arguments in the call

library("BTSPAS")  

# After which weeks is the spline allowed to jump?
demo4.jump.after <- c(22,39)  # julian weeks after which jump occurs

# Which julian weeks have "bad" recapture values. These will be set to 0 or missing prior to the run.
demo4.bad.m2     <- c(41)   # list julian weeks with bad m2 values. This is used in the Trinity Example
demo4.bad.u2     <- c(11)   # list julian weeks with bad u2 values. [This was arbitrary to demostrate the feature.]
demo4.bad.n1     <- c(38)   # list julian weeks with bad n1 values. [This was arbitrary to demonstrate the feature.]

# The prefix for the output files:
demo4.prefix <- "demo4-JC-2003-CH-TSPDE" 

# Title for the analysis
demo4.title <- "Junction City 2003 Chinook with spline for p "

cat("*** Starting ",demo4.title, "\n\n")
#> *** Starting  Junction City 2003 Chinook with spline for p

# Make the call to fit the model and generate the output files
demo4.fit <- TimeStratPetersenDiagError_fit(
                  title=demo4.title,
                  prefix=demo4.prefix,
                  time=demo4.data$jweek,
                  n1=demo4.data$n1, 
                  m2=demo4.data$m2, 
                  u2=demo4.data$u2,
                  jump.after=demo4.jump.after,
                  logitP.cov = demo4.logitP.cov,                 # ***** NEW *****
                  bad.n1=demo4.bad.n1,
                  bad.m2=demo4.bad.m2,
                  bad.u2=demo4.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: 112
#>    Total graph size: 2071
#> 
#> 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 ***

8.3 The output from the fit

Here is the fitted spline curve to the number of unmarked fish available in each recovery sample.

demo4.fit$plots$fit.plot
Fitted spline curve with covariate for p

Fitted spline curve with covariate for p

The jump in the spline when hatchery fish are released is evident.

The distribution of the posterior sample for the total number unmarked and total abundance is available as before:

demo4.fit$plots$post.UNtot.plot
Distribution of posterior samples

Distribution of posterior samples

A plot of the \(logit(P)\) is

demo4.fit$plots$logitP.plot
Estimates of logit(p)

Estimates of logit(p)

Here is a plot of the estimated \(logit(p)\)’s with the fitted spline for the \(p\)’s:

The underlying spline smooths the p’s somewhat, especially when the credible intervals are very wide (e.g. around julian weeks 37-40).

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:

demo4.fit$summary[ row.names(demo4.fit$summary) %in% c("Ntot","Utot"),]
#>         mean       sd    2.5%     25%     50%     75%   97.5%     Rhat n.eff
#> Ntot 6245617 526927.4 5428102 5881117 6181515 6545797 7439444 1.003715   830
#> Utot 6195152 526927.4 5377637 5830652 6131050 6495332 7388979 1.003711   830

The estimated total abundance from \(BTSPAS\) is 6,245,617 (SD 526,927 ) fish.

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