Learn R Programming

marked (version 1.2.8)

probitCJS: Perform MCMC analysis of a CJS model

Description

Takes design data list created with the function make.design.data for model "probitCJS" and draws a sample from the posterior distribution using a Gibbs sampler.

Usage

probitCJS(
  ddl,
  dml,
  parameters,
  design.parameters,
  burnin,
  iter,
  initial = NULL,
  imat = NULL
)

Value

A list with MCMC iterations and summarized output:

beta.mcmc

list with elements Phi and p which contain MCMC iterations for each beta parameter

real.mcmc

list with elements Phi and p which contain MCMC iterations for each real parameter

beta

list with elements Phi and p which contain summary of MCMC iterations for each beta parameter

reals

list with elements Phi and p which contain summary of MCMC iterations for each real parameter

Arguments

ddl

list of dataframes for design data; created by call to make.design.data

dml

list of design matrices created by create.dm from formula and design data

parameters

A model specification list with a list for Phi and p containing a formula and optionally a prior specification which is a named list. See 'Priors' section below.

design.parameters

Specification of any grouping variables for design data for each parameter

burnin

number of iteration to initially discard for MCMC burnin

iter

number of iteration to run the Gibbs sampler for following burnin

initial

A named list (Phi,p). If null and imat is not null, uses cjs.initial to create initial values; otherwise assigns 0

imat

A list of vectors and matrices constructed by process.ch from the capture history data

Prior distribution specification

The prior distributions used in probitCJS are multivariate normal with mean mu a and variance tau*(X'X)^(-1) on the probit scale for fixed effects. The matrix X is the design matrix based on the model specification (located in parameters$Phi$formula and parameters$p$formula respectively). Priors for random effect variance components are inverse gamma with shape parameter 'a' and rate parameter 'b'. Currently, the default values are mu=0 and tau=0.01 for phi and p parameters. For all randome effects deault values are a=2 and b=1.0E-4. In addition to the variance component each random effect can be specified with a known unscaled covariance matrix, Q, if random effects are correlated. For example, to obtain a random walk model for a serially correlated effect see Examples section below. To specify different hyper-parameters for the prior distributions, it must be done in the parameters list. See the Examples section for changing the prior distributions. Note that a and b can be vectors and the Qs are specified via a list in order of the random effects specified in the model statements.

Author

Devin Johnson

Examples

Run this code
# \donttest{
# This example is excluded from testing to reduce package check time
# Analysis of the female dipper data
data(dipper)
dipper=dipper[dipper$sex=="Female",]
# following example uses unrealistically low values for burnin and 
# iteration to reduce package testing time
fit1 = crm(dipper,model="probitCJS",model.parameters=list(Phi=list(formula=~time),
 p=list(formula=~time)), burnin=100, iter=1000)
fit1
# Real parameter summary
fit1$results$reals

# Changing prior distributions:
fit2 = crm(dipper,model="probitCJS",
   model.parameters=list(
     Phi=list(formula=~time, prior=list(mu=rep(0,6), tau=1/2.85^2)),
     p=list(formula=~time, prior=list(mu=rep(0,6), tau=1/2.85^2))
   ), burnin=100, iter=1000)
fit2
# Real parameter summary
fit2$results$reals

# Creating a Q matrix for random walk effect for 6 recapture occasions
A=1.0*(as.matrix(dist(1:6))==1)
Q = diag(apply(A,1,sum))-A

# Fit a RW survival process
fit3 = crm(dipper,model="probitCJS",
   model.parameters=list(
     Phi=list(
       formula=~(1|time), 
       prior=list(mu=0, tau=1/2.85^2, re=list(a=2, b=1.0E-4, Q=list(Q)))
       ),
     p=list(formula=~time, prior=list(mu=rep(0,6), tau=1/2.85^2))
   ), burnin=100, iter=1000)
fit3
# Real parameter summary (Not calculated for random effects yet)
fit3$results$reals

# }

Run the code above in your browser using DataLab