Capture-Recapture Estimation using Bayesian Model Averaging

Description

Performs Bayesian model averaging over all decomposable graphical models, each representing a different list dependence structure, with the goal of estimating the number of uncounted elements from the universe of elements that could have been recorded or “captured". The code here is an implementation of the model described in Madigan and York (1997).

Details

Author(s)

James Johndrow james.johndrow@gmail.com, Kristian Lum kl@hrdag.org, and Patrick Ball pball@hrdag.org

Maintainer: Kristian Lum kl@hrdag.org

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

### This simulated example goes through the whole process for 3 lists.###
library(chron)
# Simulate some data
N <- 1000 # true population size
num.lists <- 3 # number of lists

# simulate 3 lists independently 'capturing'
l1 <- rbinom(N, 1, .2)
l2 <- rbinom(N, 1, .25)
l3 <- rbinom(N, 1, .3)

overlaps <- data.frame(l1, l2, l3)

# simulate dates of recording
dates <- paste(
  rep(2015, N), "-", sample(1:4, N, replace = TRUE), "-",
  sample(1:28, N, replace = TRUE)
)
dates <- chron(dates, format = c(dates = "y-m-d"))

# save true number in each stratum for comparison later
truth <- table((months(dates)))[1:4]

# remove dates of unrecorded elements
dates <- dates[apply(overlaps, 1, sum) > 0]

# remove individuals not recorded on any list
overlaps <- overlaps[apply(overlaps, 1, sum) > 0, ]

# stratify by date
strata <- make.strata(overlaps, dates = dates, date.defs = "monthly")

# check to make sure that all strata are OK
check <- check.strata(strata)

# look at strata, just to make sure
par(mfrow = c(2, 2), mar = rep(1, 4))
for (i in 1:nrow(strata$overlap.counts)) {
  venn3(strata$overlap.counts[i, ],
    main = rownames(strata$overlap.counts)[i],
    cex.main = 1
  )
}


# load the graphs to make the estimates
data(graphs3)

# select expansion factor defining the largest number of unrecorded elements.
# this makes Nmissing <- 0:(sum(Y)*fac)
fac <- 5

# set prior
delta <- 1 / 2^num.lists

# loop over strata to calculate posterior distributions of
# the total population size for each stratum
par(mfrow = c(2, 2), mar = rep(1.75, 4))
# if using Rstudio, make sure your plot window is pretty big here!
for (i in 1:nrow(strata$overlap.counts)) {
  Nmissing <- 0:(sum(strata$overlap.counts[i, ]) * fac)
  Y <- array(strata$overlap.counts[i, ], dim = rep(2, num.lists))
  weights <- bma.cr(Y, Nmissing, delta, graphs3)
  plotPosteriorN(weights, Nmissing + sum(strata$overlap.counts[i, ]),
    main = rownames(strata$overlap.counts)[i]
  )
  points(truth[i], .5 * max(weights), col = "red", pch = 16, cex = 2)
}

Computes Marginal Likelihoods for Each Clique and Value of Nmissing

Description

Assembles all of the pieces of the marginal likelihoods to be used to calculate the posterior probability of each model/value of Nmissing.

Usage

CompLogML(D, Nmissing, delta)

Arguments

Value

The log marginal likelihood of the marginal table.

Author(s)

James Johndrow james.johndrow@gmail.com and Kristian Lum kl@hrdag.org

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

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))

# Compute marginal array over lists 1 and 3
D <- apply(Y, c(1, 3), sum)

dga:::CompLogML(D, 1:300, 0.5)

Component-wise Matrix of Log Marginal Likelihoods

Description

Calls CompLogML to create a matrix of number of possible components by length(Nmissing) log marginal likelihoods. Calculates the log marginal likehood of each possible marginal table for every value of Nmissing.

Usage

MakeCompMatrix(p, delta, Y, Nmissing)

Arguments

Value

A matrix of log marginal likelihoods.

Author(s)

James Johndrow james.johndrow@gmail.com and Kristian Lum kl@hrdag.org

Bayesiam Model Averaging for Capture-Recapture

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

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 james.johndrow@gmail.com and Kristian Lum (kl@hrdag.org)

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)

A Helper Function for make.strata

Description

A helper function used in make.strata to make list overlap counts.

Usage

cfunction(x, nlist)

Arguments

Value

a table of the number of records with each capture history

Author(s)

Kristian Lum kl@hrdag.org

Examples

## The function is currently defined as
cfunction <- function(x, nlist) {
  out <- table(c(x, 0:(2^nlist - 1))) - 1
}

Checks Each Stratum for Suitability for Capture-Recapture

Description

Takes in list intersection counts and source list totals as produced by make.strata. It then checks whether there are between three and five lists, whether all of the lists are non-empty, and whether all of the lists overlap with some other list.

Usage

check.strata(strata)

Arguments

Value

A boolean indicating whether any serious problems have been found with the strata.

Note

This does not issue a warning for cases where some subset of lists is not connected to the others, e.g. Lists A and B have overlap with each other, lists C and D have overlap with each other, but no records from A or B overlap with lists C or D. We suggest that you examine the list intersection counts manually as well.

Author(s)

Kristian Lum kl@hrdag.org

Examples

library(chron)

N <- 1000
overlaps <- data.frame(l1 = rbinom(N, 1, .5), l2 = rbinom(N, 1, .5), l3 = rbinom(N, 1, .5))
dates <- paste(
  rep(2015, N), "-", sample(1:12, N, replace = TRUE), "-",
  sample(1:28, N, replace = TRUE)
)
dates <- chron(dates, format = c(dates = "y-m-d"))
locations <- sample(c("A", "B", "C", "D"), N, replace = TRUE)

# Aggregate only by week:
strata <- make.strata(overlaps, dates, date.def = "weekly")
check <- check.strata(strata)

A Helper Function Used in venn3

Description

Takes the parameters of a circle and returns points on its perimeter to be plotted to make circles for a venn diagram.

Usage

circle(x, y, r)

Arguments

Value

Author(s)

Kristian Lum kl@hrdag.org

Examples

plot(dga:::circle(0, 0, 1), type = "l")

A Helper Function for venn3

Description

Used in venn3 to tell whether proposed points are inside of the given circle.

Usage

circle.ind(ps, x, y, r)

Arguments

Value

a length n vector telling whether each row of ps is inside the given circle.

Author(s)

Kristian Lum kl@hrdag.org

Examples

ps <- cbind(runif(100), runif(100))
plot(dga:::circle(0, 0, 1), type = "l")
inds <- dga:::circle.ind(ps, 0, 0, 1)
points(inds)

A Helper Function Used by Venn4 to Define the Perimeter of an Ellipse

Description

Draws the ellipses used in venn4.

Usage

ellipse(x, y, a, b, alpha)

Arguments

Value

points that define the perimeter of an ellipse.

Author(s)

Kristian Lum kl@hrdag.org

Examples

plot(dga:::ellipse(0, 0, .5, .2, 1))

A Helper Function Used by Venn4

Description

Takes potential points to be plotted in the venn diagrams and returns whether the point is inside or outside of the ellipse described by x, y, a, b, and alpha.

Usage

ellipse.ind(ps, x, y, a, b, alpha)

Arguments

Value

a length n vector indicating whether each point is inside the ellipse.

Author(s)

Kristian Lum kl@hrdag.org

Examples

## The function is currently defined as
ps <- cbind(runif(100), runif(100))
plot(dga:::ellipse(0, 0, .5, .3, 0), type = "l")
inds <- dga:::ellipse.ind(ps, 0, 0, .5, .3, 0)
points(inds)

All Decomposable Graphical Models on Three Lists

Description

This dataset contains all of the cliques and separators for each of the decomposable graphical models on three lists. On three lists, this is all of the models.

Usage

data(graphs3)

Format

A list of lists. graphs3[[i]] is the ith model under consideration. This consists of graphs3[[i]]$C, all of the cliques in that model, and graphs3[[i]]$S, the separators.

All Decomposable Graphical Models on Four Lists

Description

This dataset contains all of the cliques and separators for each of the decomposable graphical models on four lists.

Usage

data(graphs4)

Format

A list of lists. graphs4[[i]] is the ith model under consideration. This consists of graphs4[[i]]$C, all of the cliques in that model, and graphs4[[i]]$S, the separators.

All Decomposable Graphical Models on Five Lists

Description

This dataset contains all of the cliques and separators for each of the decomposable graphical models on five lists.

Usage

data(graphs5)

Format

A list of lists. graphs5[[i]] is the ith model under consideration. This consists of graphs5[[i]]$C, all of the cliques in that model, and graphs5[[i]]$S, the separators.

Base Converter

Description

Takes a decimal number and converts it to base b.

Usage

integer.base.b(x, b = 2)

Arguments

Details

This was harvested from the internet here: https://stat.ethz.ch/pipermail/r-help/2003-September/038978.html. Posted by Spencer Graves.

Value

A number in base b.

Author(s)

Spencer Graves

References

https://stat.ethz.ch/pipermail/r-help/2003-September/038978.html

Transforms Records to List Intersection Counts by Stratum

Description

Helps you to create list overlaps in the correct order to be used in bma.cr. This function also does some of the heavy lifting to stratify records by time (date, etc.) and other variables.

Usage

make.strata(
  overlaps,
  dates = NULL,
  locations = NULL,
  demographics = NULL,
  date.defs = "monthly",
  loc.defs = NULL,
  demog.defs = NULL,
  start.date = NULL,
  end.date = NULL
)

Arguments

Value

Author(s)

Kristian Lum kl@hrdag.org

Examples

require(chron)
N <- 10000
overlaps <- data.frame(l1 = rbinom(N, 1, .5), l2 = rbinom(N, 1, .5), l3 = rbinom(N, 1, .5))
dates <- paste(
  rep(2015, N), "-", sample(1:12, N, replace = TRUE), "-",
  sample(1:28, N, replace = TRUE)
)
dates <- chron(dates, format = c(dates = "y-m-d"))
locations <- sample(c("A", "B", "C", "D"), N, replace = TRUE)

# Aggregate only by week:
make.strata(overlaps, dates, date.def = "weekly")

# Aggregate by year and location, where locations are not grouped:
make.strata(overlaps, dates, date.def = "yearly", locations)

# Aggregate by 2 day increments and location, where there are unique location levels
#       A, B, C, and D and locations A and B are in group 1
#       and locations C and D are in group 2.
loc.defs <- list("g1" = c("A", "B"), "g2" = c("C", "D"))
make.strata(overlaps, dates, date.def = 2, locations, loc.defs = loc.defs)
# Aggregate by demographic (sex) only, where sex takes values M, F, A, NA, and U
#       and we would like to group these as M, F, and other.
sex <- sample(c("M", "F", "A", NA, "U"),
  prob = c(.4, .4, .1, .05, .05),
  N, replace = TRUE
)
demog.defs <- list("M" = "M", "F" = "F", "Other" = c("A", NA, "U"))
make.strata(overlaps, demographics = sex, demog.defs = demog.defs)

Plots Posterior Distribution of Nmissing

Description

Plots the model averaged posterior distribution of the total number of elements (the solid line) and the contribution to the posterior of each of the models (dotted lines)

Usage

plotPosteriorN(weights, N, main = NULL)

Arguments

Value

A plot.

Author(s)

Kristian Lum kl@hrdag.org

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)

A Helper Function to Tell Which Points Are Near the Boundary of a Circle

Description

Used in venn3 to tell which of the potential points to be plotted are near the boundary of the circle defned by x, y, and r.

Usage

remove.close(ps, x, y, r)

Arguments

Value

Author(s)

Kristian Lum kl@hrdag.org

Examples

ps <- cbind(runif(100), runif(100))
inds <- dga:::remove.close(ps, .5, .5, .1)

A Helper Function to Tell Which Points are Near the Boundary of the Ellipse

Description

A helper function.

Usage

remove.close.ellipse(ps, x, y, a, b, alpha)

Arguments

Value

Author(s)

Kristian Lum kl@hrdag.org

Examples

## The function is currently defined as
ps <- cbind(runif(100), runif(100))
inds <- dga:::remove.close.ellipse(ps, .5, .5, .1, .3, 1)

A Helper Function for make.strata.

Description

This is the simplest function ever. It's just an apply to sum across columns.

Usage

sfunction(x)

Arguments

Value

apply(x, 2, sum)

Author(s)

Kristian Lum kl@hrdag.org

Examples

## The function is currently defined as
sfunction <- function(x) {
  out <- apply(x, 2, sum)
}

Three List Venn Diagram

Description

A function that plots a venn diagram of 3 lists. One point is plotted in each region for each record that falls into the corresponding list overlap.

Usage

venn3(
  overlap.counts,
  main = NULL,
  num.test.points = 1e+05,
  p.cex = 0.75,
  write_numbers = FALSE,
  t.cex = 1.25,
  cex.main = 1
)

Arguments

Value

a 3-way venn diagram with points inside of each segment representing the number of records on each list overlap.

Author(s)

Kristian Lum kl@hrdag.org

Examples

overlap.counts <- rpois(8, 30)
venn3(overlap.counts, main = "example diagram")

Four List Venn Diagram

Description

A function that plots a venn diagram of 4 lists. One point is plotted in each region for each record that falls into the corresponding list overlap.

Usage

venn4(
  overlap.counts,
  main = NULL,
  num.test.points = 1e+05,
  p.cex = 0.75,
  cex.main = 1
)

Arguments

Value

A venn diagram of the list overlap structure for four lists. Each region of the plot contains points representing each record in that list intersection.

Author(s)

Kristian Lum kl@hrdag.org

Examples

overlap.counts <- rpois(16, 50)
venn4(overlap.counts, main = "example diagram")