multimodelClosed {multimark}

Multimodel inference for 'multimark' closed population abundance models

Description

This function performs Bayesian multimodel inference for a set of 'multimark' closed population abundance models using the reversible jump Markov chain Monte Carlo (RJMCMC) algorithm proposed by Barker & Link (2013).

Usage

multimodelClosed(
  modlist,
  modprior = rep(1/length(modlist), length(modlist)),
  monparms = "N",
  miter = NULL,
  mburnin = 0,
  mthin = 1,
  M1 = NULL,
  pbetapropsd = 1,
  zppropsd = NULL,
  sigppropshape = 6,
  sigppropscale = 4,
  printlog = FALSE
)

Arguments

modlist

A list of individual model output lists returned by multimarkClosed or markClosed. The models must have the same number of chains and MCMC iterations.

modprior

Vector of length length(modlist) containing prior model probabilities. Default is modprior = rep(1/length(modlist), length(modlist)).

monparms

Parameters to monitor. Only parameters common to all models can be monitored (e.g., "pbeta[(Intercept)]", "N"), but derived capture ("p") and recapture ("c") probabilities can also be monitored. Default is monparms = "N".

miter

The number of RJMCMC iterations per chain. If NULL, then the number of MCMC iterations for each individual model chain is used.

mburnin

Number of burn-in iterations (0 <= mburnin < miter).

mthin

Thinning interval for monitored parameters.

M1

Integer vector indicating the initial model for each chain, where M1_j=i initializes the RJMCMC algorithm for chain j in the model corresponding to modlist[[i]] for i=1,..., length(modlist). If NULL, the algorithm for all chains is initialized in the most general model. Default is M1=NULL.

pbetapropsd

Scaler specifying the standard deviation of the Normal(0, pbetapropsd) proposal distribution for "pbeta" parameters. Default is pbetapropsd=1. See Barker & Link (2013) for more details.

zppropsd

Scaler specifying the standard deviation of the Normal(0, zppropsd) proposal distribution for "zp" parameters. Only applies if at least one (but not all) model(s) include individual hetergeneity in detection probability. If NULL, zppropsd = sqrt(sigma2_zp) is used. Default is zppropsd=NULL. See Barker & Link (2013) for more details.

sigppropshape

Scaler specifying the shape parameter of the invGamma(shape = sigppropshape, scale = sigppropscale) proposal distribution for sigma_zp. Only applies if at least one (but not all) model(s) include individual hetergeneity in detection probability. Default is sigppropshape=6. See Barker & Link (2013) for more details.

sigppropscale

Scaler specifying the scale parameter of the invGamma(shape = sigppropshape, scale = sigppropscale) proposal distribution for sigma_zp. Only applies if at least one (but not all) model(s) include individual hetergeneity in detection probability. Default is sigppropscale=4. See Barker & Link (2013) for more details.

printlog

Logical indicating whether to print the progress of chains and any errors to a log file in the working directory. Ignored when nchains=1. Updates are printed to log file as 1% increments of iter of each chain are completed. With >1 chains, setting printlog=TRUE is probably most useful for Windows users because progress and errors are automatically printed to the R console for "Unix-like" machines (i.e., Mac and Linux) when printlog=FALSE. Default is printlog=FALSE.

Details

Note that setting parms="all" is required when fitting individual multimarkClosed or markClosed models to be included in modlist.

Value

A list containing the following:

rjmcmc

Reversible jump Markov chain Monte Carlo object of class mcmc.list. Includes RJMCMC output for monitored parameters and the current model at each iteration ("M").

pos.prob

A list of calculated posterior model probabilities for each chain, including the overall posterior model probabilities across all chains.

Author(s)

Brett T. McClintock

References

Barker, R. J. and Link. W. A. 2013. Bayesian multimodel inference by RJMCMC: a Gibbs sampling approach. The American Statistician 67: 150-156.

See Also

multimarkClosed, markClosed, processdata

Examples

# This example is excluded from testing to reduce package check time
# Example uses unrealistically low values for nchain, iter, and burnin

#Generate object of class "multimarksetup"
setup <- processdata(bobcat)
 
#Run single chain using the default model for bobcat data. Note parms="all".
bobcat.dot <- multimarkClosed(mms=setup,parms="all",iter=1000,adapt=500,burnin=500)

#Run single chain for bobcat data with time effects. Note parms="all".
bobcat.time <- multimarkClosed(mms=setup,mod.p=~time,parms="all",iter=1000,adapt=500,burnin=500)

#Perform RJMCMC using defaults
modlist <- list(mod1=bobcat.dot,mod2=bobcat.time)
bobcat.M <- multimodelClosed(modlist=modlist,monparms=c("N","p"))

#Posterior model probabilities
bobcat.M$pos.prob
 
#multimodel posterior summary for abundance
summary(bobcat.M$rjmcmc[,"N"])