bma.cr: Bayesiam Model Averaging for Capture-Recapture

In HRDAG/DGA: Capture-Recapture Estimation using Bayesian Model Averaging

Description

This function averages over all decomposable graphical models for p lists.

Usage

bma.cr(
  Y,
  Nmissing,
  delta,
  graphs,
  logprior = NULL,
  log.prior.model.weights = NULL
)

Arguments

Y	a 2^p array of list intersection counts. See details.
Nmissing	A vector of all possible values for the number of individuals that appear on no list.
delta	The hyper-parameter for the hyper-Dirichlet prior distribution on list intersection probabilities. A smaller value indicates fewer prior observations per cell. A suggested default is 2^-p
graphs	A pre-computed list of all decomposable graphical models for p lists. These should be loaded using data(graphsp); see example. Currently, this package includes a list of graphs for three, four, or five lists.
logprior	The log of the prior probability of each value in Nmissing. If left blank, this will default to the -log(Nmissing).
log.prior.model.weights	Prior weights on the graphs. This should be a vector of length length(graphs).

Details

This is the main function in this package. It performs capture-recapture (or multiple systems estimation) using Bayesian model averaging as outlined in Madigan and York (1997).

Y can be created by the array() command from a vector that is ordered lexigraphically by the cell names, e.g., c(x000, x001, x010, x011, x100, x101, x110, x111).

Value

This function returns a matrix of weights, where rows correspond to models and columns correspond to values of Nmissing. Thus, the ijth entry of the matrix is the posterior probability of the ith model and the jth entry of Nmissing. Row sums return posterior probabilities by model.Column sums return posterior probabilities by value of Nmissing.

Note

This function is pretty robust relative to the more common log-linear model approach to capture-recapture. It will not fail (or issue a numerical warning) even if there are no overlaps among the lists. The user should take care that there is adequate list overlap and that there are sufficient cases in the stratum.

Author(s)

James Johndrow [e-mail] and Kristian Lum ([e-mail])

References

Madigan, David, and Jeremy C. York. "Bayesian methods for estimation of the size of a closed population." Biometrika 84.1 (1997): 19-31.

Examples

#### 5 list example from M & Y ##########
delta <- .5
Y <- c(0, 27, 37, 19, 4, 4, 1, 1, 97, 22, 37, 25, 2, 1, 3, 5,
       83, 36, 34, 18, 3, 5, 0, 2, 30, 5, 23, 8, 0, 3, 0, 2)
Y <- array(Y, dim = c(2, 2, 2, 2, 2))
Nmissing <- 1:300
N <- Nmissing + sum(Y)
data(graphs5)
weights <- bma.cr(Y, Nmissing, delta, graphs5)
plotPosteriorN(weights, N)

##### 3 list example from M & Y #######
Y <- c(0, 60, 49, 4, 247, 112, 142, 12)
Y <- array(Y, dim = c(2, 2, 2))

delta <- 1
a <- 13.14
b <- 55.17


Nmissing <- 1:300
N <- Nmissing + sum(Y)

logprior <- N * log(b) - (N + a) * log(1 + b) + lgamma(N + a) - lgamma(N + 1) - lgamma(a)

data(graphs3)
weights <- bma.cr(Y, Nmissing, delta, graphs3, logprior)
plotPosteriorN(weights, N)

HRDAG/DGA documentation built on Dec. 17, 2021, 10:27 p.m.