R: Wrapper (*_fit) to call the Time Stratified Petersen...

LP_BTSPAS_fit_NonDiag {Petersen}

R Documentation

Wrapper (*_fit) to call the Time Stratified Petersen Estimator with NON-Diagonal Entries function in BTSPAS.

Description

Takes the data structure as described below, and uses Bayesian methods to fit a fit a spline through the population numbers and a hierarchical model for the trap efficiency over time. An MCMC object is also created with samples from the posterior.

Usage

LP_BTSPAS_fit_NonDiag(
  data,
  p_model = ~1,
  p_model_cov = NULL,
  jump.after = NULL,
  logitP.fixed = NULL,
  logitP.fixed.values = NULL,
  InitialSeed = ceiling(stats::runif(1, min = 0, max = 1e+06)),
  n.chains = 3,
  n.iter = 2e+05,
  n.burnin = 1e+05,
  n.sims = 2000,
  trace = FALSE,
  remove_MCMC_files = TRUE,
  quietly = FALSE
)

Arguments

`data`	Data frame containing the variables: cap_hist Capture history (see details below) freq Number of times this capture history was observed plus any other covariates (e.g. discrete strata and/or continuous covariates) to be used in the model fitting.
`p_model`	Model for the captured probabilities. This can reference other variables in the data frame, plus a special reserved term `..time` to indicate a time dependence in the capture probabilities. For example, `p_model=~1` would indicate that the capture probabilities are equal across the sampling events; `p_model=~..time` would indicate that the capture probabilities vary by sampling events; `p_model=~sex*..time` would indicate that the capture probabilities vary across all combination of sampling events (`..time`) and a stratification variable (`sex`). The `sex` variable also needs to be in the data frame. For some models (e.g., tag loss models), the `..time` variable cannot be used because the models only have one capture probability (e.g., only for event 1).
`p_model_cov`	Data frame with covariates for the model for prob capture at second sampling event. If this data frame is given, it requires one line for each of the temporal strata at the second sampling event (even if missing in the `data` that has the capture histories) with one variable being `..time` to represent the second temporal stratum.
`jump.after`	A numeric vector with elements belonging to `time`. In some cases, the spline fitting the population numbers should be allowed to jump. For example, the population size could take a jump when hatchery released fish suddenly arrive at the trap. The jumps occur AFTER the strata listed in this argument.
`logitP.fixed`	A numeric vector (could be null) of the time strata where the logit(P) would be fixed. Typically, this is used when the capture rates for some strata are 0 and logit(P) is set to -10 for these strata. The fixed values are given in `logitP.fixed.values`
`logitP.fixed.values`	A numerical vector (could be null) of the fixed values for logit(P) at strata given by logitP.fixed. Typically this is used when certain strata have a 0 capture rate and the fixed value is set to -10 which on the logit scale gives $p_i$ essentially 0. Don't specify values such as -50 because numerical problems could occur when this is converted to the 0-1 scale.
`InitialSeed`	Numeric value used to initialize the random numbers used in the MCMC iterations.
`n.chains`	Number of chains to fit in the MCMC
`n.iter`	Total number of iterations
`n.burnin`	Number of burnin iterations
`n.sims`	Total number of simulations to keep in output (implies a thinning)
`trace`	Internal tracing flag.
`remove_MCMC_files`	Should the temporary MCMC files (init.txt, data.text, model.txt, CODA*txt) removed after the fit.
`quietly`	Suppress all console messages that occur during the fit. This includes the progress bar when a model that requires MCMC is fit (LP_BTSPAS_fit_Diag and LP_BTSPAS_fit_NonDiag), or a trace of the likelihood during the fit (LP_SPAS_fit).

Details

Use the Petersen::LP_BTSPAS_fit_Diag function for cases where recaptures take place in a single stratum (diagonal case).

The frequency variable (freq in the data argument) is the number of animals with the corresponding capture history.

Capture histories (cap_hist in the data argument) are character values of the format xx..yy is a capture_history where xx and yy are the temporal stratum (e.g., julian week) and '..' separates the two temporal strata. If a fish is released in temporal stratum and never captured again, then yy is set to 0; if a fish is newly captured in temporal stratum yy, then xx is set to zero. For example, a capture history of 23..23 indicates animals released in temporal stratum 23 and recaptured in temporal stratum 23; a capture history of 23..00 indicates animals released in temporal stratum 23 and never seen again; a capture history of 00..23 indicates animals newly captured in temporal stratum 23 at the second sampling event.

In the non-diagonal case, fish are allowed to move among temporal strata.

It is not necessary to label the temporal strata starting at 1; BTSPAS will treat the smallest value of the temporal strata seen as the first stratum and will interpolate for temporal strata without any data. Temporal strata labels should be numeric, i.e., do NOT use A, B, C etc.

Value

An list object of class LP_BTSPAS_fit_Diag with the following elements

summary A data frame with the information on the number of observations in the fit
data Data used in the fit
p_model, p_model_cov Information on modelling the capture probabilities at the second occasion
fit n MCMC object with samples from the posterior distribution. A series of graphs and text file are also created with summary information. Refer to the BTSPAS package for more details.
datetime Date and time the fit was done

References

Bonner, S. J. and Schwarz, C. J. (2021). BTSPAS: Bayesian Time Stratified Petersen Analysis System.R package version 2021.11.2.

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. doi:10.1111/j.1541-0420.2011.01599.x

Examples


# NOTE. To keep execution time to a small value as required by CRAN
# I've made a very small example.
# Additionally, I've set the number of MCMC chains, iterations, burning, simulation to save to
# small values. Proper mixing may not have occurred yet.
# When using this routine, you likely want to the use the default values
# for these MCMC parameters.
data(data_btspas_nondiag1)
temp<- cbind(data_btspas_nondiag1,
             split_cap_hist( data_btspas_nondiag1$cap_hist,
                             sep="..", make.numeric=TRUE))
xtabs(~t1, data=temp)

# only use data up to week 10 to keep example small
temp <- temp[ temp$t1 %in% c(0, 27:32) & temp$t2 %in% c(0, 27:32),]

fit <- Petersen::LP_BTSPAS_fit_NonDiag(
  temp,
  p_model=~1,
  InitialSeed=23943242,
  # the number of chains and iterations are too small to be useful
  # they are set to a small number to pare execution time to <5 seconds for an example
  n.chains=2, n.iter=20000, n.burnin=1000, n.sims=100,
  quietly=TRUE
)
fit$summary

# now get the estimates of abundance
est <-  Petersen::LP_BTSPAS_est (fit)
est$summary

[Package Petersen version 2023.12.1 Index]