multimodelCJS {multimark}

Multimodel inference for 'multimark' open population survival models

Description

This function performs Bayesian multimodel inference for a set of 'multimark' open population survival (i.e., Cormack-Jolly-Seber) models using the reversible jump Markov chain Monte Carlo (RJMCMC) algorithm proposed by Barker & Link (2013).

Usage

multimodelCJS(
  modlist,
  modprior = rep(1/length(modlist), length(modlist)),
  monparms = "phi",
  miter = NULL,
  mburnin = 0,
  mthin = 1,
  M1 = NULL,
  pbetapropsd = 1,
  zppropsd = NULL,
  phibetapropsd = 1,
  zphipropsd = NULL,
  sigppropshape = 1,
  sigppropscale = 0.01,
  sigphipropshape = 1,
  sigphipropscale = 0.01,
  printlog = FALSE
)

Arguments

modlist

A list of individual model output lists returned by multimarkCJS. 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)]", "phibeta[(Intercept)]", "psi"), but derived survival ("phi") and capture ("p") probabilities can also be monitored. Default is monparms = "phi".

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.

phibetapropsd

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

zphipropsd

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

sigppropshape

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

sigppropscale

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

sigphipropshape

Scaler specifying the shape parameter of the invGamma(shape = sigphipropshape, scale = sigphipropscale) proposal distribution for "sigma2_zphi". Only applies if at least one (but not all) model(s) include individual hetergeneity in survival probability. Default is sigphipropshape=1. See Barker & Link (2013) for more details.

sigphipropscale

Scaler specifying the scale parameter of the invGamma(shape = sigphipropshape, scale = sigphipropscale) proposal distribution for "sigma_zphi". Only applies if at least one (but not all) model(s) include individual hetergeneity in survival probability. Default is sigphipropscale=0.01. 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 multimarkCJS 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

multimarkCJS, 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" from simulated data
data_type = "always"
noccas <- 5
phibetaTime <- seq(2,0,length=noccas-1) # declining trend in survival
data <- simdataCJS(noccas=5,phibeta=phibetaTime,data.type=data_type)
setup <- processdata(data$Enc.Mat,data.type=data_type)

#Run single chain using the default model. Note parms="all".
sim.pdot.phidot <- multimarkCJS(mms=setup,parms="all",iter=1000,adapt=500,burnin=500)

#Run single chain with temporal trend for phi. Note parms="all".
sim.pdot.phiTime <- multimarkCJS(mms=setup,mod.phi=~Time,parms="all",iter=1000,adapt=500,burnin=500)

#Perform RJMCMC using defaults
modlist <- list(mod1=sim.pdot.phidot,mod2=sim.pdot.phiTime)
sim.M <- multimodelCJS(modlist=modlist)

#Posterior model probabilities
sim.M$pos.prob

#multimodel posterior summary for survival (display first cohort only)
summary(sim.M$rjmcmc[,paste0("phi[1,",1:(noccas-1),"]")])