---
title: "A User Guide to MultipleSystemsEstimation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A User Guide to MultipleSystemsEstimation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(MultipleSystemsEstimation)
```

## Introduction

Multiple systems estimation uses several incomplete lists of observed individuals to estimate the size of a population, including the number of individuals who appear on none of the lists, often referred
to as the *dark figure*. The `MultipleSystemsEstimation` package implements Poisson log-linear methods for this problem, with particular attention to two important practical issues:

* sparse capture data, including combinations of lists with no observed overlap; and
* uncertainty introduced by choosing a log-linear model from the observed data.

For most analyses the recommended starting point is the high-level function

```r
estimate_population(zdat)
```

whose default `method = "auto"` chooses an estimation method according to the number of lists. For data with up to five lists it uses BIC model selection; for six or more lists it uses the stepwise procedure. The individual methods can also be selected explicitly.

The BIC/bootstrap approach follows Silverman, Chan and Vincent (2024). The stepwise procedure and the treatment of sparse-data existence and identifiability issues build on Chan, Silverman and Vincent (2021).
The package also implements the Bayesian thresholding approach of Section 6 of Silverman (2020).

This vignette is intended as an overall introduction to the package. Five
companion vignettes give further details and, in some cases, paper-specific
workflows:

* [Capture histories and hierarchical models](captures-and-models.html):
  capture notation, hierarchical models and sparse interactions;
* [BIC model selection and bootstrap inference](bic-model-selection-and-bootstrap-inference.html):
  BIC methodology and routine use;
* [Reproducing the 2024 BIC bootstrap investigations](reproduce-2024-bic-bootstrap.html):
  selected computational investigations from the paper;
* [Sparse data, stepwise model selection and the bootstrap](sparse-data-and-stepwise.html):
  the sparse-data and stepwise methodology and selected published results; and
* [Bayesian-thresholding multiple systems estimation](bayesian-threshold.html):
  prior specification, thresholding and higher-order interactions.


## Data format

Suppose there are \(t\) lists. Each observed individual has a *capture history* specifying the lists on which that individual occurs.

Data supplied to the package have \(t+1\) columns. The first \(t\) columns contain zeros and ones indicating membership of the lists, and the final column contains the number of individuals observed with that capture history. Observable capture histories omitted from the data are treated as having count zero.

The all-zero capture history is fundamentally different. It represents individuals who occur on none of the lists and therefore cannot be observed directly. Estimating the size of this unobserved group, the 'dark figure', is the central objective of multiple systems estimation.

The package includes several example datasets. The three-list `Korea` data provide a compact illustration:

```{r korea-data}
data(Korea)
Korea
```

## Start with `estimate_population()`

For routine use, begin with the default call:

```{r auto-fit}
fit <- estimate_population(Korea)
fit$estimate
fit$method
```

Because `Korea` has three lists, `method = "auto"` selects the BIC method. The returned `estimate` gives both the estimated dark figure and the estimated total population.

BIC enumeration is generally practical only for up to five lists. Six-list BIC
analyses are computationally demanding even when restricted to pairwise
interactions, and enumeration beyond six lists is impractical
because of the rapid growth in the number of models. The `auto` method therefore
selects stepwise fitting for six or more lists. See the appendix
to [Capture histories and hierarchical models](captures-and-models.html)
for model counts for different numbers of lists and orders of models.

The rest of this vignette explains the available methods and the common inferential ideas behind them.

## Specifying log-linear models

For an observable capture history $\omega$, the log-linear model is

$$
N_\omega \sim \operatorname{Poisson}(\mu_\omega),
\qquad
\log \mu_\omega =
\sum_{\substack{\theta \in \mathcal{M}\\ \theta \subseteq \omega}}
\alpha_\theta.
$$

Here $\mathcal{M}$ is the set of terms included in the model.
The empty set $\varnothing$ corresponds to the intercept
$\alpha_\varnothing$. The fitted expected count for the unobserved
all-zero capture history is $\exp(\hat\alpha_\varnothing)$,
and the population estimate is the total observed count plus
this estimated unobserved count.

Throughout the package, we consider only *hierarchical* log-linear
models. In such a model, whenever an interaction is included,
all its lower-order interactions and main effects are also included.
The intercept and all main effects are always present. Thus, for
example, including the three-way interaction between lists 1, 2
and 3 requires the inclusion of all three pairwise interactions
between these lists.

Hierarchical models can be specified by listing their highest-order
terms in square brackets. For example, the model `[123, 24, 5]`
includes the three-way interaction between lists 1, 2 and 3,
the pairwise interaction between lists 2 and 4, and all five
main effects, together with the intercept and the lower-order
interactions implied by hierarchy.

The [Capture histories and hierarchical models](captures-and-models.html)
vignette explains this notation in more detail, together with
the generation of hierarchical models and the treatment of
sparse interactions.

## BIC model selection

The BIC approach is the method used automatically for data with up to five lists. It is based on the methodology of Silverman, Chan and Vincent (2024). For details of the method, see [BIC model selection and bootstrap inference](bic-model-selection-and-bootstrap-inference.html) and the help page for [estimate_population_bic()](../help/estimate_population_bic).

A list of candidate hierarchical Poisson log-linear models is generated, in principle covering every possible model up to a specified maximum order, except those that cannot be fitted or are guaranteed to be outranked by other models. The candidates are fitted to the capture-history counts and ranked by Bayesian information criterion (BIC). The model with the smallest BIC gives the point estimate for the original data.

The BIC method can be requested explicitly:

```{r bic-point}
fit_bic <- estimate_population_bic(Korea)

fit_bic$estimate
fit_bic$fitted_model
```

With the appropriate calling arguments, the `bic_results` component contains the full list of models 
considered and their BIC values:

```{r bic-results}
fit_bic <- estimate_population_bic(Korea, return_details = TRUE)

fit_bic$details$bic_results$res
```
### Which models are considered?

The argument `maxorder` controls the maximum order of interactions considered by the method. If it is left at its default `NULL`, the package chooses:

| Number of lists | Automatic `maxorder` |
|----------------:|---------------------:|
| 2               | 1 |
| 3               | 2 |
| 4--5            | 3 |
| 6               | 2 |

Candidate models are generated as needed. 
With six lists, even with `maxorder = 2` there are 32,768 possible
pairwise hierarchical models before excluding any models, which explains the
automatic choice of stepwise fitting.  


## Bootstrap inference after BIC selection

Choosing the best model and then calculating a confidence interval as though that model had been fixed in advance ignores the uncertainty arising from the model-selection step. The bootstrap approach of Silverman, Chan and Vincent (2024) addresses this by repeating BIC model selection within bootstrap samples.

Bootstrap inference is requested through `nboot`; to keep the vignette quick, the examples below use only 100 bootstrap replications. 
A much larger number, for example 1000 or more, should be used for substantive inference.

```{r bic-bootstrap}
fit_bic_boot <- estimate_population_bic(Korea,nboot = 100)

fit_bic_boot$uncertainty
```

### Restricting model selection within the bootstrap

Repeating an exhaustive model search for every bootstrap sample can be expensive. The 2024 methodology therefore ranks models by their BIC on the original data and allows bootstrap model selection to be restricted to the best-ranked models.

With `ntopmodels = NULL`, the current defaults are:

| Number of lists | Models retained for bootstrap inference |
|----------------:|:-----------------------------------------|
| 2--3            | All available models                     |
| 4               | Top 20 by original-data BIC              |
| 5--6            | Top 100 by original-data BIC             |

These defaults affect bootstrap and jackknife calculations only. The original-data point estimate is always based on the minimum-BIC model from the full candidate set allowed by `maxorder`.

The number retained can be specified directly, for example:

```{r bic-ntop}
fit_bic_5 <- estimate_population_bic( Korea,
  nboot = 100,
  ntopmodels = 5
)

fit_bic_5$uncertainty
```

The `uncertainty` component gives the endpoints of the bias-corrected and accelerated (BCa) confidence intervals obtained when model selection is allowed among all models retained for the bootstrap calculation.

## Sparse capture data

Sparse capture data require particular care. For example, suppose
two lists A and B have no observed individuals in common. If the
interaction AB is included in a fitted model, the likelihood
increases towards its supremum as the corresponding parameter
tends to $-\infty$. It is therefore natural to work in an
*extended maximum likelihood* framework, in which the estimated
parameter for AB explicitly takes the value $-\infty$.
The fitted mean for the capture history AB, and for all capture
histories containing both A and B, is then zero.

AB is called a *minimal sparse interaction* because, by assumption,
every individual list contains some observations, so neither A
nor B is sparse. (Completely empty lists should simply be removed
from consideration.)

Now suppose AB is a sparse interaction and its parameter has
been set to $-\infty$. If ABC is also included in the model,
no finite value assigned to the parameter for ABC can affect
the fitted means. ABC is called a *masked sparse interaction*
because it has a proper subset, AB, which is itself sparse.

The package handles sparse interactions explicitly. Fuller
details are given in
[Capture histories and hierarchical models](captures-and-models.html).

Two further problems can occur:

1. The maximum likelihood estimate may fail to exist, even
   allowing extended boundary values.
2. The model parameters may fail to be identifiable.

The package's fitting routines check
for both, so that unsuitable models are excluded from model
selection and population estimation.
The treatment of these issues builds on the general results
of Fienberg and Rinaldo (2012), specialised to multiple
systems estimation by Chan, Silverman and Vincent (2021).


## The stepwise method

The second principal model-selection method is an extension of the forward
stepwise procedure of Chan, Silverman and Vincent (2021). It can be requested
explicitly for any supported number of lists and is selected automatically
when there are six or more lists.

The procedure starts with the main-effects model. At each stage it considers
interactions that can be added while preserving the hierarchical structure of
the model. Candidates that fail the existence or identifiability checks
described above are discarded. Among the remaining candidates, the interaction
with the smallest p-value is added if that p-value is below the threshold
`pthresh`; otherwise the procedure stops.
The procedure uses extended maximum likelihood, allowing
minimal sparse interactions to be included in fitted models
with their parameters estimated at $-\infty$.

The default is `pthresh = 0.02`, following the simulation study reported by
Chan, Silverman and Vincent (2021). Setting `maxorder = 2` means that only
two-list interactions are considered and the procedure agrees with the method
studied in that paper.

For example:

```{r stepwise-point}
data(Kosovo)

fit_pairwise <- estimate_population_stepwise(
  Kosovo,
  maxorder = 2
)

fit_pairwise$estimate
fit_pairwise$fitted_model
```

The threshold can be changed using `pthresh`; for example, 
setting `pthresh = 0.01` makes the procedure more conservative about adding interactions.

By default, `maxorder = Inf`, so interactions of any order
can be considered. Setting `maxorder = 2` restricts the
procedure to pairwise interactions. Higher-order interactions
become eligible only when all their lower-order terms are
already present. Masked sparse interactions are excluded
from consideration, since their parameters cannot affect
the fitted means.

For the Kosovo data, allowing interactions of any order changes both the
selected model and the population estimate:

```{r stepwise-higher-order}
fit_step_unrestricted <- estimate_population_stepwise(
  Kosovo,
  maxorder = Inf
)

c(
  pairwise = fit_pairwise$estimate["total"],
  unrestricted = fit_step_unrestricted$estimate["total"]
)

c(
  pairwise = fit_pairwise$fitted_model,
  unrestricted = fit_step_unrestricted$fitted_model
)
```

The pairwise search selects `[12,13,14,23,34]`, whereas the unrestricted
search selects `[134,12,23]`. The latter hierarchy includes interactions 13,
14 and 34 implicitly and adds the three-list interaction 134. The corresponding
population estimates are approximately 14,342 and 18,393.

The vignette
[Sparse data, stepwise model selection and the bootstrap](sparse-data-and-stepwise.html)
describes the published pairwise procedure in more detail.

## Bootstrap inference after stepwise selection


Bootstrap inference for the stepwise method follows the same general
approach as for BIC model selection, described above. In each
bootstrap replicate, the stepwise selection procedure is repeated,
so that the resulting inference takes account of model-selection
uncertainty. BCa confidence intervals are calculated in the same way.


```{r stepwise-bootstrap}
fit_step_boot <- estimate_population_stepwise(
  Kosovo,
  maxorder = Inf,
  nboot = 100)

fit_step_boot$uncertainty
```



Thus the BCa 95% confidence interval (based on only 100 bootstrap replications) for the total population is
`r format(round(fit_step_boot$uncertainty["total", "0.025"]), scientific = FALSE)` to
`r format(round(fit_step_boot$uncertainty["total", "0.975"]), scientific = FALSE)`.

## Fitting a specified model


Sometimes a particular log-linear model is to be fitted without
using either BIC or stepwise model selection. In that case,
use `method = "fixed"` or the routine `estimate_population_fixed()`.


For example:

```{r fixed-point}
fit_fixed <- estimate_population(
  Korea,
  method = "fixed",
  model = "[12,23]"
)

fit_fixed$estimate
fit_fixed$fitted_model
```

Bootstrap inference for a fixed model is requested in the same way. In contrast to the BIC and stepwise procedures, the same specified model is fitted to every bootstrap and jackknife sample. There is no model-selection step inside the resampling procedure.

```{r fixed-bootstrap}
fit_fixed_boot <- estimate_population(
  Korea,
  method = "fixed",
  model = "[12,23]",
  nboot = 100,
  return_details = TRUE
)

fit_fixed_boot$uncertainty
fit_fixed_boot$details$bootstrap_failed
```

With only 100 replications, the interval illustrates the calculation
and should not be used for substantive inference. The displayed count records bootstrap replications that failed
to give a finite population estimate and were omitted. Such failures
can arise because resampling changes the pattern of zero counts,
so that the specified model no longer yields a finite estimate.

For all three likelihood-based methods (BIC, stepwise and fixed),
bootstrap inference is optional, with `nboot = 0` by default.
When bootstrap inference is requested, model selection is repeated
in each replicate for BIC and stepwise estimation, whereas the
specified model is retained for fixed-model estimation.
All three methods use BCa confidence intervals, with the
acceleration term calculated by a delete-one jackknife.
The random-number seed can be controlled using `iseed`, and
the probability levels for confidence limits can be changed
using `alpha`.



## The Bayesian thresholding method

The package also provides the Bayesian thresholding method of Silverman
(2020, Section 6). This can be selected by setting
`method = "bayesthresh"` or by calling
`estimate_population_bayesthresh()` directly. It is not selected
by `method = "auto"`.

The procedure starts with the log-linear model containing all
two-list interactions, fitted using Markov chain Monte Carlo (MCMC).
For each interaction parameter, the absolute value of its posterior
mean is divided by its posterior standard deviation. Interactions
for which this ratio is below `threshold` are discarded, and
the model containing the retained interactions is refitted.

Improper priors are used for the intercept and main effects.
By default, proper normal priors are used for the interaction
parameters, although improper priors can also be specified.
The MCMC calculations use the `MCMCpack` package.

```{r bayesthresh-point}
fit_bayes <- estimate_population_bayesthresh(
  Kosovo, maxorder = 3
)

fit_bayes$estimate
fit_bayes$fitted_model
fit_bayes$uncertainty
```

The returned `uncertainty` component gives posterior quantiles
for both the dark figure and the total population. These provide
the endpoints of Bayesian credible intervals, rather than the
bootstrap confidence intervals used by the other methods.

In this example, the 95% credible interval for the total population
is `r formatC(fit_bayes$uncertainty["total", "0.025"], format = "f", digits = 0, big.mark = ",")`
to `r formatC(fit_bayes$uncertainty["total", "0.975"], format = "f", digits = 0, big.mark = ",")`.

By default, `maxorder = 2`, so only two-list interactions are
considered. Setting `maxorder = 3` allows three-list interactions,
while `maxorder = 4` also permits four-list interactions.
A higher-order interaction becomes eligible only if all its
immediate lower-order interactions were retained at the preceding
stage. At each stage, thresholding is applied only to the
newly introduced interactions; lower-order interactions are
not reconsidered.

When improper priors are specified for the interaction parameters,
the package checks models for identifiability and the
Fienberg--Rinaldo existence criterion. If the initial model
containing all two-list interactions fails either check, the
procedure stops. If a proposed higher-order extension fails
the existence check, the procedure returns the result from
the preceding stage. These existence checks are not required
when proper priors are used for the interaction parameters.
With improper priors, masked sparse interactions are also
excluded from consideration.

For fuller details of the method, including prior specification,
higher-order interactions and the treatment of sparse data, see
[Bayesian-thresholding multiple systems estimation](bayesian-threshold.html).
The general treatment of sparse interactions is explained in
[Capture histories and hierarchical models](captures-and-models.html).


## Which method should I use?

For most users the answer is simple: begin with

```r
estimate_population(zdat)
```

and therefore allow `method = "auto"` to select the route.

There are, however, reasons to choose explicitly:

* use `method = "bic"` when BIC model selection is specifically wanted and the number of lists makes exhaustive enumeration feasible;
* use `method = "stepwise"` when the Chan--Silverman--Vincent stepwise procedure is preferred, when sensitivity to `pthresh` is of interest, or for larger numbers of lists where BIC enumeration is unavailable or unattractive;
* use `method = "fixed"` when the model has been specified independently of these data-driven selection procedures;
* use `method = "bayesthresh"` if you prefer the Bayesian thresholding approach.


## Looking inside the result

All four methods share the following top-level components:

```{r names-result}
names(fit_bic)
```

* `input`: the original call and data;
* `method`: the method actually used, including after `method = "auto"`;
* `estimate`: estimates of the dark figure and total population;
* `fitted_model`: the supplied, selected or retained hierarchy;
* `uncertainty`: confidence or credible endpoints for both estimates, or an
  explanatory message when bootstrap inference was not requested;
* `details`: method-specific supporting output, or `"not requested"`.

The fixed method also returns `masked_effects`.

Set `return_details = TRUE` to retain supporting output. For BIC this includes
the selected BIC and complete original-data enumeration; for the likelihood
methods it includes the bootstrap estimates and BCa acceleration. The complete
conditional GLM fit is returned only by the fixed method. Bayesian details
include the thresholding information and complete posterior sample of total
population size.

The individual help pages give the complete definitions of returned components:

```r
?estimate_population
?estimate_population_bic
?estimate_population_stepwise
?estimate_population_fixed
?estimate_population_bayesthresh
```

## References

Chan, L., Silverman, B. W. and Vincent, K. (2021). Multiple Systems Estimation for Sparse Capture Data: Inferential Challenges When There Are Nonoverlapping Lists. *Journal of the American Statistical Association*, **116**(535), 1297--1306. [doi:10.1080/01621459.2019.1708748.](https://doi.org/10.1080/01621459.2019.1708748).

DiCiccio, T. J. and Efron, B. (1996). Bootstrap Confidence Intervals. *Statistical Science*, **11**(3), 189--228.

Fienberg, S. E. and Rinaldo, A. (2012). Maximum Likelihood Estimation in Log-Linear Models. *The Annals of Statistics*, **40**, 996--1023.

Silverman, B. W. (2020). Multiple-systems analysis for the quantification of modern slavery: classical and Bayesian approaches. *Journal of the Royal Statistical Society: Series A (Statistics in Society)*, **183**(3), 691--736. [doi:10.1111/rssa.12505.](https://doi.org/10.1111/rssa.12505).

Silverman, B. W., Chan, L. and Vincent, K. (2024). Bootstrapping Multiple Systems Estimates to Account for Model Selection. *Statistics and Computing*, **34**, 44. [doi:10.1007/s11222-023-10346-9.](https://doi.org/10.1007/s11222-023-10346-9).