---
title: "BIC model selection and bootstrap inference"

output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{BIC model selection and bootstrap inference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Introduction

BIC provides a standard criterion for selecting among log-linear models
for multiple systems estimation. Inference conditional on a single selected
model, however, fails to reflect uncertainty about the model selection.
Silverman, Chan and Vincent (2024) investigate bootstrap inference in
which model selection is repeated within each resample.

This vignette describes the supported estimation and inference facilities
in `MultipleSystemsEstimation`, including enhancements beyond the paper.
Detailed reproductions of its numerical investigations are in the separate
vignette [*Reproducing the 2024 BIC bootstrap investigations*](reproduce-2024-bic-bootstrap.html).

## BIC model selection

The Bayesian information criterion (BIC) provides a standard method
for selecting among statistical models, balancing goodness of fit
against model complexity. For multiple systems estimation, the
candidate models are hierarchical log-linear models for the observed
capture histories. The model with the lowest BIC is selected, and
the corresponding fitted model is used to estimate the unobserved
population.

The package generates the candidate hierarchical models as required,
rather than relying on a stored catalogue. The maximum interaction
order can be specified using `maxorder`, allowing the search to be
restricted to models of a specified complexity. The intercept and
all main effects are always included.

A novel feature of the implementation is the way it takes advantage
of sparsity in the observed data by excluding hierarchical models
containing masked sparse interactions. To see why this is justified,
consider a model containing such interactions. Replacing it by its
non-masked reduction, obtained by removing all masked sparse
interactions, leaves the maximised log likelihood unchanged but
reduces the number of parameters, and hence decreases the BIC.

This reduces the number of candidate models and avoids considering
models containing parameters that cannot be identified from the
observed data. Models containing minimal sparse interactions remain eligible.
The distinction between minimal and masked sparse
interactions, and the construction of non-masked reductions, are
explained in the vignette [*Capture histories and hierarchical models*](captures-and-models.html).


The public interface for BIC model selection is
`estimate_population_bic()`, or equivalently
`estimate_population(..., method = "bic")`. For example:

```{r bic-example}
data(Korea)

bic_fit <- estimate_population_bic(
  Korea, maxorder = 2, return_details = TRUE
)

bic_fit$estimate
bic_fit$fitted_model
```

The result contains the estimated population size, the selected
model and other information about the estimation procedure.
The restriction `maxorder = 2` allows pairwise interactions but
excludes higher-order interactions.


### Inspecting the fit

The input is a capture-history data frame or matrix. Its first columns
indicate membership of each list using zeros and ones, and its final column
contains the number of observed individuals with each capture history.
The all-zero history represents individuals observed on no list; their
number is inferred from the fitted model rather than supplied as an
observed count.

The component `estimate` distinguishes the estimated number observed on
no list (`dark_figure`) from the estimated total population (`total`).
The latter adds the observed population to the dark figure. The component
`fitted_model` identifies the selected hierarchy of interactions.

Specifying `return_details = TRUE`, as in the example above, makes the
original-data model comparisons available:

```{r inspect-bic-fit}
head(bic_fit$details$bic_results)
```

This table is ordered by increasing BIC and includes the population
estimate and maximum interaction order for each eligible model. It allows
the selected model to be viewed alongside its competitors. Models with
similar BIC values can nevertheless give different population estimates;
reporting only the estimate from the best model conceals that aspect of
model uncertainty.

### Choosing the interaction order

The argument `maxorder` defines the model space before BIC selection.
For example, `maxorder = 1` permits only the main-effects model, whereas
`maxorder = 2` permits hierarchical models with pairwise interactions.
Increasing the order allows more complex dependence between lists, but
can greatly increase the number of candidate models. The choice should
reflect both the data and the scientific context; BIC selects within the
specified space, rather than establishing that the space itself is adequate.

If `maxorder` is omitted, the package uses order 1 for two lists, order 2
for three lists, order 3 for four or five lists, and order 2 for six or more
lists. These defaults provide a starting point rather than a substitute
for considering the dependence structure.

The argument `maxmodel` limits the number of candidate hierarchical models generated. If the number would exceed this limit, the maximum interaction order is reduced, if necessary down to 1, until the number of models is within the limit. Thus the maximum order actually searched may be lower than the requested `maxorder`. The resulting search includes all eligible hierarchical models at the reduced maximum order, rather than an arbitrary subset of models at the originally requested order.
When `return_details = TRUE`, the requested and actual maximum orders are recorded in `details$requested_maxorder` and `details$maxorder`, respectively.

## Bootstrap inference

Confidence intervals based on a single selected model do not account
for the uncertainty arising from model selection. The bootstrap
methodology developed in Silverman, Chan and Vincent (2024) addresses
this by repeating model selection within each bootstrap sample,
rather than treating the model selected from the original data as fixed.

However, a full BIC search over all candidate models for every
bootstrap sample would often be computationally prohibitive.
The key computational device within the paper and this package
is to rank the candidate models by
their BIC values for the original data, and then restrict model
selection within each bootstrap sample to the first `ntop` models
in this ranking.

With `ntop = 1`, the same model is used for every bootstrap sample,
so no allowance is made for model selection uncertainty. Increasing
`ntop` allows more models to be selected within bootstrap samples,
while avoiding repeated searches over the entire model space.
The choice `ntop = Inf` corresponds to unrestricted model selection
within each bootstrap sample.

The bootstrap distribution of the resulting population estimates
is used to construct confidence intervals. The methodology uses
bias-corrected and accelerated (BCa) intervals, which adjust for
bias and asymmetry, making use of jackknife calculations in which
observations are omitted in turn.

The set of `ntop` models is determined from the original data and
is retained throughout the bootstrap and jackknife calculations.
Resampling can introduce additional masked sparse interactions, but
cannot remove any already present in the original data. Although
some of the retained models may therefore contain newly masked
interactions, the candidate set is not revised: removing such models
would introduce an additional data-dependent selection step, and
their non-masked reductions have the same maximised likelihood.
All `ntop` models are therefore considered in every resample.

Bootstrap inference is provided by `estimate_population_bic()`, or via `estimate_population()` with `method= "bic"` and `nboot > 0`. The specialised internal computations used for the paper are documented in [*Reproducing the 2024 BIC bootstrap investigations*](reproduce-2024-bic-bootstrap.html).

The package also supports bootstrap inference after stepwise model selection; see [*Sparse data, stepwise model selection and the bootstrap*](sparse-data-and-stepwise.html).

The choice of `ntop` involves a trade-off between computational
effort and the extent to which the bootstrap reproduces full
model selection. The effect of this trade-off is investigated by Silverman, Chan and Vincent (2024).

## Practical bootstrap inference

The basic example above uses `nboot = 0`, the default, and therefore
returns a point estimate without bootstrap inference. To obtain BCa
confidence limits, specify a positive number of bootstrap replications:

```{r bootstrap-example}
bic_boot <- estimate_population_bic(
  Korea,
  maxorder = 2,
  nboot = 10,
  iseed = 1234
)

bic_boot$estimate
bic_boot$uncertainty
```

Only 10 bootstrap replications are used here to keep the example quick.
A substantially larger number, such as 1000, should be used for substantive
inference. The argument `ntopmodels` controls how many models from the
original-data BIC ranking are retained for bootstrap and jackknife
calculations.

### Reading the confidence limits

The component `uncertainty` is a matrix with rows for the dark figure and
total population. Its columns correspond to the cumulative probability
levels supplied through `alpha`. The default levels, 0.025, 0.1, 0.9 and
0.975, give endpoints for 95% and 80% BCa confidence intervals. To request
only a 95% interval, use `alpha = c(0.025, 0.975)`.

The point estimate continues to be based on the model selected from the
original data. The confidence limits, however, use population estimates
obtained after repeating selection within the retained model set for each
bootstrap sample. They therefore describe a different source of uncertainty
from an interval that conditions on the original selected model throughout.

For substantive analyses, the reported endpoints should be checked for
numerical stability. The argument `iseed` makes the random sampling
reproducible for a given implementation and set of settings, but
reproducibility with a fixed seed does not establish numerical stability.

### Choosing the retained model set

The notation `ntop` above corresponds to the public argument `ntopmodels`.
It controls the models retained for bootstrap and jackknife calculations,
not the initial original-data search or its point estimate.

* `ntopmodels = 1` retains only the best original-data model, and therefore
  does not allow for model selection uncertainty.
* A larger finite value allows selection among that many of the
  best-ranked original-data models, or all available models if fewer exist.
* `ntopmodels = Inf` retains all eligible models in the chosen model space.
* The default, `ntopmodels = NULL`, retains all available models for two or
  three lists, up to 20 for four lists, and up to 100 for five or more lists.

For example, the following displays the call for retaining all models and
requesting a 95% interval with 1000 replications. It is not evaluated during
the vignette build:

```{r full-bootstrap-call, eval=FALSE}
bic_full <- estimate_population_bic(
  Korea,
  maxorder = 2,
  ntopmodels = Inf,
  nboot = 1000,
  alpha = c(0.025, 0.975),
  iseed = 1234,
  return_details = TRUE
)

bic_full$uncertainty
```

For a small model space, retaining all models avoids the approximation
introduced by a finite restriction. For larger spaces, a finite value can
make the bootstrap practical. Comparing confidence limits over increasing
values of `ntopmodels` can help assess sensitivity to the restriction.
Confidence limits need not change monotonically as models are added, and
agreement between two choices alone does not guarantee agreement with
unrestricted selection. Increasing `nboot` addresses Monte Carlo variability;
increasing `ntopmodels` addresses the restriction of the candidate set.

When `return_details = TRUE`, the result also includes the bootstrap
population estimates in `details$bootstrap_estimates` and the BCa
acceleration in `details$bca_acceleration`. These are available for further
inspection without rerunning the fit. The intervals account for selection
within the specified model space and retained set; they do not remove the
need to assess the underlying capture-history data and modelling assumptions.

## Summary

The package applies standard BIC selection within an on-the-fly generated
hierarchical model space, and exploits sparsity by excluding models
containing masked sparse interactions while retaining minimal sparse
interactions. For inference that reflects model selection uncertainty,
it implements the restricted-model bootstrap and BCa methodology of
Silverman, Chan and Vincent (2024). The [reproduction vignette](reproduce-2024-bic-bootstrap.html)
provides code for the numerical investigations underlying the paper.

## References

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