BIC model selection and bootstrap inference • MultipleSystemsEstimation

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.

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.

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

data(Korea)

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

bic_fit$estimate
#> dark_figure       total 
#>    34.16667   157.16667
bic_fit$fitted_model
#> [1] "[12,23]"

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:

head(bic_fit$details$bic_results)
#> $res
#>         abundance       BIC modelsorder
#> [12,23]  157.1667  57.14081           2
#> [12,3]   268.7778  58.96860           2
#> [13,23]  123.4630  91.89711           2
#> [23,1]   126.1944 138.06368           2
#> [13,2]   126.9394 155.16018           2
#> [1,2,3]  141.9926 184.14431           1
#> 
#> $xdata
#>      b c d Count
#> [1,] 0 0 1    41
#> [2,] 0 1 0     5
#> [3,] 1 0 0     5
#> [4,] 1 1 0    54
#> [5,] 1 0 1     6
#> [6,] 0 1 1     0
#> [7,] 1 1 1    12
#> 
#> $maxorder
#> [1] 2
#> 
#> $models
#> $models$`[12,23]`
#> [1] 1 2 3 4 5 7
#> 
#> $models$`[12,3]`
#> [1] 1 2 3 4 5
#> 
#> $models$`[13,23]`
#> [1] 1 2 3 5 6 7
#> 
#> $models$`[23,1]`
#> [1] 1 2 3 5 7
#> 
#> $models$`[13,2]`
#> [1] 1 2 3 5 6
#> 
#> $models$`[1,2,3]`
#> [1] 1 2 3 5
#> 
#> 
#> $requested_maxorder
#> [1] 2
#> 
#> $enumeration_maxorder
#> [1] 2

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.

The package also supports bootstrap inference after stepwise model selection; see Sparse data, stepwise model selection and the bootstrap.

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:

bic_boot <- estimate_population_bic(
  Korea,
  maxorder = 2,
  nboot = 10,
  iseed = 1234
)

bic_boot$estimate
#> dark_figure       total 
#>    34.16667   157.16667
bic_boot$uncertainty
#>                 0.025       0.1      0.9 0.975
#> dark_figure  11.72727  11.82012 197.3872   213
#> total       134.72727 134.82012 320.3872   336

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:

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