---
title: "Sparse data, stepwise model selection and the bootstrap"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Sparse data, stepwise model selection and the bootstrap}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Introduction

The published paper Chan, Silverman and Vincent (2021) develops methodology for multiple systems estimation when the observed capture table is sparse.

The paper has two principal themes. First, sparse tables require explicit attention to the existence and identifiability of maximum likelihood estimates, including the possibility of extended maximum likelihood estimates in which some log-linear parameters are equal to \(-\infty\). Second, the paper develops a BCa bootstrap approach that takes account of model selection, in the context of a forward stepwise model-selection procedure.

This vignette describes the treatment of sparse capture data and the stepwise model-selection and BCa inference procedure implemented in `MultipleSystemsEstimation`, using examples from the paper to illustrate the methods and reproduce selected results.

The package has developed these approaches further since the original paper was published. Where possible, this vignette uses the current public interface, and the original simulation and supplementary routines are not all carried forward into `MultipleSystemsEstimation`. For full reproduction of the original computational material, including routines specific to the paper, see the historical `SparseMSE` package (now deprecated for current analyses). 

## Sparse capture data and extended maximum likelihood

The general treatment of capture histories, hierarchical log-linear models and sparse interactions is given in [Capture histories and hierarchical models](captures-and-models.html). Here we concentrate on the consequences of sparsity for estimation and stepwise model selection.

Two lists are said to be *nonoverlapping* if no individual appears on both. The corresponding pairwise interaction is then minimal sparse, and its extended maximum likelihood estimate is minus infinity. The package handles such cases explicitly, rather than relying on numerical iteration towards increasingly negative values.

This treatment builds on the general theory of Fienberg and Rinaldo (2012), specialised to multiple systems estimation in Chan, Silverman and Vincent (2021).

### Existence and identifiability

Section 2.6 of the paper distinguishes two problems:

1. the maximum likelihood estimate may fail to exist, even allowing extended parameter values; and
2. the parameters may be nonidentifiable.

These are illustrated by the artificial three-list dataset in Table 2, which is supplied with the package as `Artificial_3`.

```{r}
data(Artificial_3)
Artificial_3
```

The current function `check_extended_MLE()` checks both conditions for a specified hierarchical model. It returns 0 if both conditions are satisfied, 1 for failure of the Fienberg--Rinaldo existence condition, 2 for failure of identifiability, and 3 if both fail.

The model containing all three two-list interactions is nonidentifiable, while the model containing only the AB interaction fails the Fienberg--Rinaldo existence condition. The model containing AC and BC, but not AB, passes both tests.

```{r}
MultipleSystemsEstimation::check_extended_MLE(Artificial_3, "[12,13,23]")

MultipleSystemsEstimation::check_extended_MLE(Artificial_3, "[12,3]")

MultipleSystemsEstimation::check_extended_MLE(Artificial_3, "[13,23]")

```

These examples reproduce the three distinct behaviours displayed in Table 3 of the paper.

## The stepwise model-selection procedure

Section 3.2 of the paper proposes a forward stepwise procedure.

Starting with the main-effects model, every two-list effect not already included is considered. A candidate is discarded if adding it would lead to nonexistence or nonidentifiability. 
Among the remaining candidates, the effect with the smallest p-value is selected if that p-value is no greater than a specified threshold. The procedure stops when no remaining candidate meets this criterion.

The paper ultimately recommends the threshold `pthresh = 0.02`.

The current high-level implementation is
`estimate_population_stepwise()`. The package extends the
original pairwise procedure to hierarchical models containing
higher-order interactions. By default, `maxorder = Inf`,
allowing interactions of any order to be considered.

For a point estimate without bootstrapping, restricting the model to pairwise interactions:

```{r}
fit <- estimate_population_stepwise(Korea, pthresh = 0.02, maxorder = 2)

fit$estimate
fit$fitted_model
```

The same method is also available from the overall package interface:

```{r}
fit <- estimate_population(
  Korea,
  method = "stepwise",
  pthresh = 0.02,
  maxorder = 2
)
```

### Allowing higher-order interactions

At each iteration, an interaction is eligible for inclusion only
if all its lower-order interactions are already included in the
model. Masked sparse interactions are excluded from consideration.
Eligible interactions of different orders are considered together,
rather than completing selection at one order before proceeding
to the next. Thus the procedure preserves hierarchy throughout.

For further details of hierarchy and masked sparse interactions,
see [Capture histories and hierarchical models](captures-and-models.html).

For the Kosovo data, the pairwise and unrestricted searches give different
results:

```{r stepwise-higher-order}
data(Kosovo)

Kosovo_pairwise <- estimate_population_stepwise(
  Kosovo,
  pthresh = 0.02,
  maxorder = 2
)

Kosovo_unrestricted <- estimate_population_stepwise(
  Kosovo,
  pthresh = 0.02,
  maxorder = Inf
)

c(
  pairwise = Kosovo_pairwise$estimate["total"],
  unrestricted = Kosovo_unrestricted$estimate["total"]
)

c(
  pairwise = Kosovo_pairwise$fitted_model,
  unrestricted = Kosovo_unrestricted$fitted_model
)
```

The pairwise search selects `[12,13,14,23,34]`, whereas the unrestricted
search selects `[134,12,23]`. The latter model includes interactions 13, 14
and 34 implicitly, so its additional maximal interaction is the three-list
term 134. The corresponding population estimates are approximately 14,342
and 18,393 respectively.

## Bootstrap inference

Section 3.3 points out that confidence intervals conditional on the model selected from the observed data do not allow for uncertainty due to model selection.

The paper therefore repeats the complete stepwise selection procedure for each bootstrap sample.

The bootstrap samples are obtained by multinomial resampling of the observed capture histories. The BCa acceleration is calculated using a delete-one jackknife.

```{r Korea_boot_1}
fit_boot <- estimate_population_stepwise(
  Korea,
  nboot = 100,
  pthresh = 0.02,
  maxorder = 2,
  iseed = 1234
)

fit_boot$estimate
fit_boot$uncertainty
```

To keep the vignette quick, the bootstrap calculations within this vignette use only 100 replications. A substantially larger number of bootstrap replications should be used for substantive inference. The resulting BCa confidence intervals are included only to demonstrate the operation of the code and may differ considerably from the published intervals quoted in the text, which are based on 1000 replications.  

## Applications and published results

### New Orleans: Section 4.1

The New Orleans dataset in Section 4.1 contains eight lists.

```{r NewOrl_boot_1}
data(NewOrl)

ncol(NewOrl) - 1
sum(NewOrl[, ncol(NewOrl)])
```

There are 28 possible pairs of lists, of which 18 are nonoverlapping.

Applying the stepwise procedure with `pthresh = 0.02`:

```{r}
NewOrl_fit <- estimate_population_stepwise(
  NewOrl,
  pthresh = 0.02, 
  maxorder = 2
)

NewOrl_fit$estimate
NewOrl_fit$fitted_model
```

The paper reports that the selected model contains the two-list effect DE and gives a total population estimate of 1184.

The published BCa 95% confidence interval, based on 1000 bootstrap replications, 
is (717, 1657).

The corresponding calculation with the current package is:

```{r NewOrl_boot2, eval = FALSE}
NewOrl_boot <- estimate_population_stepwise(
  NewOrl,
  pthresh = 0.02,
  maxorder = 2, 
  nboot = 1000
)

NewOrl_boot$uncertainty
```
The bootstrap code for the full eight-list New Orleans data is displayed but not run at all, because even a small value of `nboot` leads to slow computation.

The paper also considers the main-effects-only model, obtained when the threshold is 0.01 or smaller. It reports a point estimate of 997 and a 95% confidence interval of (644, 1618).

Setting `pthresh = 0` ensures that no interactions are selected, giving the main-effects model. This can be examined as follows:

```{r NewOrl_boot_3, eval = FALSE}
NewOrl_main <- estimate_population_stepwise(
  NewOrl,
  pthresh = 0,
  maxorder = 2, 
  nboot = 1000
)

NewOrl_main$estimate
NewOrl_main$uncertainty
```

### Five-list New Orleans data

Because several individual lists in the New Orleans data are very small, the paper also combines the four smallest lists to obtain a five-list version. This is supplied as `NewOrl_5`.

```{r NewOrl_boot_4}
data(NewOrl_5)

NewOrl5_fit <- estimate_population_stepwise(
  NewOrl_5,
  pthresh = 0.02,
  maxorder = 2
)

NewOrl5_fit$estimate
NewOrl5_fit$fitted_model
```

The paper reports that no two-list parameter is significant even at the 5% level. The point estimate is 1034, with BCa 95% confidence interval (589, 1703).

The bootstrap calculation is:

```{r NewOrl_boot_6}
NewOrl5_boot <- estimate_population_stepwise(
  NewOrl_5,
  pthresh = 0.02,
  maxorder = 2, 
  nboot = 100
)

NewOrl5_boot$uncertainty
```

### Western site data: Section 4.2

The Western site data contain five lists.

```{r Western_boot_7}
data(Western)

Western_fit <- estimate_population_stepwise(
  Western,
  pthresh = 0.02,
  maxorder = 2
)

Western_fit$estimate
Western_fit$fitted_model
```

The paper reports that the selected model contains the two-list effect AE and that the resulting point estimate is 2483.

The reported BCa 95% confidence interval is (1293, 3670).

It can be recalculated using:

```{r Western_boot_9}
Western_boot <- estimate_population_stepwise(
  Western,
  pthresh = 0.02,
  maxorder = 2, 
  nboot = 100
)

Western_boot$uncertainty
```

## Simulation studies

### Comparison with the BIC approach: Section 3.4

Section 3.4 compares the stepwise/BCa procedure with the BIC approach then implemented in `Rcapture`.

The simulation starts from the model selected by the stepwise procedure for the five-list UK data. The fitted population size and capture-history probabilities are treated as the parameters of a simulated population, and 500 datasets are generated.

For every simulated dataset:

- the BIC procedure selects a model and calculates the confidence interval available in `Rcapture`;
- the stepwise procedure selects a model afresh;
- a BCa interval is constructed for the stepwise estimate, repeating the selection within its bootstrap samples.

The original computational routines used for this simulation are retained in the now-deprecated `SparseMSE` package rather than in `MultipleSystemsEstimation`. `SparseMSE` also retains other routines used for full reproduction of the computational and supplementary material associated with the 2021 paper.

The paper reports population-scale root mean square errors of 3057 for the stepwise method and 5834 for BIC. On the logarithmic scale the corresponding values are 0.19 and 0.34.

The paper also reports that the nominal 95% BCa interval covered the true value in 90% of the 500 simulations, compared with 61.4% for the `Rcapture` method. For nominal 80% intervals, the corresponding coverage rates were approximately 70% and 42.8%.

The BIC procedure used in this comparison was that implemented in `Rcapture` at the time of the original investigation. It should not be confused with the BIC model-selection and bootstrap procedure implemented in `MultipleSystemsEstimation`, which incorporates subsequent developments; see [BIC model selection and bootstrap inference](bic-model-selection-and-bootstrap-inference.html).

The purpose of this simulation was not to establish universal superiority of one procedure, but to demonstrate the importance of allowing model-selection uncertainty to enter the inference.

### Choosing the threshold: Section 4.3

The paper uses a larger simulation study to investigate the choice of the stepwise threshold.

Seven datasets are used:

- the UK data, in full and five-list versions;
- the Netherlands data, in full and five-list versions;
- the New Orleans data, in full and five-list versions;
- the Western site data.

For each dataset four generating models are considered:

- the full two-list model;
- the main-effects-only model;
- the model selected at threshold 0.001;
- the model selected at threshold 0.05.

This gives 28 simulation scenarios.

For each scenario the paper generates 1000 datasets and obtains estimates using thresholds
0, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, 0.1, 1,
where zero corresponds to the main-effects model and one to the full two-list model.

Accuracy is assessed using the mean squared error of the logarithm of the estimated dark figure. Because the overall level of mean squared error differs substantially between scenarios, the paper combines the results by averaging the logarithm of the mean squared error across the 28 scenarios.
The minimum overall score occurs at `pthresh = 0.02`, which motivates the default threshold used by the package.

The full threshold simulation is deliberately not run as part of the vignette. It involves 28 generating scenarios, 1000 datasets per scenario, and repeated model fitting across multiple thresholds, so it is better regarded as a published simulation study than as an ordinary package example.

The original routines used to carry out this simulation study are retained in the now-deprecated `SparseMSE` package for full reproducibility of the paper.

## Relation to the current package

The original `SparseMSE` package accompanied Chan, Silverman and Vincent (2021). It is now deprecated for current analysis, but is retained for historical reproducibility. The present `MultipleSystemsEstimation` package provides a unified interface to this and other estimation methods.

The default interface is `estimate_population(..., method = "auto")`, which uses the BIC-based approach for up to five lists and the stepwise approach for six or more. The main user guide discusses the choice of methods in more detail.


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