Sparse data, stepwise model selection and the bootstrap

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

the maximum likelihood estimate may fail to exist, even allowing extended parameter values; and
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.

data(Artificial_3)
Artificial_3
#>   A B C  n
#> 1 1 0 0 40
#> 2 0 1 0 30
#> 3 0 0 1 20
#> 4 1 1 0  6

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.

MultipleSystemsEstimation::check_extended_MLE(Artificial_3, "[12,13,23]")
#> [1] 2

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

MultipleSystemsEstimation::check_extended_MLE(Artificial_3, "[13,23]")
#> [1] 0

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:

fit <- estimate_population_stepwise(Korea, pthresh = 0.02, maxorder = 2)

fit$estimate
#> dark_figure       total 
#>    145.7778    268.7778
fit$fitted_model
#> [1] "[12,3]"

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

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.

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

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"]
)
#>     pairwise.total unrestricted.total 
#>           14341.66           18393.31

c(
  pairwise = Kosovo_pairwise$fitted_model,
  unrestricted = Kosovo_unrestricted$fitted_model
)
#>           pairwise       unrestricted 
#> "[12,13,14,23,34]"      "[134,12,23]"

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.

fit_boot <- estimate_population_stepwise(
  Korea,
  nboot = 100,
  pthresh = 0.02,
  maxorder = 2,
  iseed = 1234
)

fit_boot$estimate
#> dark_figure       total 
#>    145.7778    268.7778
fit_boot$uncertainty
#>                   0.025       0.1      0.9   0.975
#> dark_figure   0.3274105  70.82153 235.3031 276.828
#> total       123.3274105 193.82153 358.3031 399.828

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.

data(NewOrl)

ncol(NewOrl) - 1
#> [1] 8
sum(NewOrl[, ncol(NewOrl)])
#> [1] 185

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

Applying the stepwise procedure with pthresh = 0.02:

NewOrl_fit <- estimate_population_stepwise(
  NewOrl,
  pthresh = 0.02, 
  maxorder = 2
)

NewOrl_fit$estimate
#> dark_figure       total 
#>    998.6923   1183.6923
NewOrl_fit$fitted_model
#> [1] "[45,1,2,3,6,7,8]"

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:

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:

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.

data(NewOrl_5)

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

NewOrl5_fit$estimate
#> dark_figure       total 
#>    849.1511   1034.1511
NewOrl5_fit$fitted_model
#> [1] "[1,2,3,4,5]"

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:

NewOrl5_boot <- estimate_population_stepwise(
  NewOrl_5,
  pthresh = 0.02,
  maxorder = 2, 
  nboot = 100
)

NewOrl5_boot$uncertainty
#>                0.025      0.1      0.9    0.975
#> dark_figure 416.0089 416.0089 1110.799 1651.069
#> total       601.0089 601.0089 1295.799 1836.069

Western site data: Section 4.2

The Western site data contain five lists.

data(Western)

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

Western_fit$estimate
#> dark_figure       total 
#>    2138.384    2483.384
Western_fit$fitted_model
#> [1] "[15,2,3,4]"

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:

Western_boot <- estimate_population_stepwise(
  Western,
  pthresh = 0.02,
  maxorder = 2, 
  nboot = 100
)

Western_boot$uncertainty
#>                0.025      0.1      0.9   0.975
#> dark_figure 1205.954 1205.954 2524.327 3064.25
#> total       1550.954 1550.954 2869.327 3409.25

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.

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.