Quantifying uncertainty with fitPS

Two kinds of uncertainty

A fitted probability model gives point estimates, but those estimates come from a finite survey. In fitPS, it is useful to distinguish two related questions:

  1. Parameter uncertainty: which values of the fitted model parameters remain plausible after observing the survey?
  2. Probability uncertainty: how much uncertainty does that imply for the fitted P or S probabilities used in an application?

The second question is often the more important one in practice. Model parameters are intermediate quantities; P and S probabilities are usually the quantities that enter the scientific or forensic calculation.

This vignette uses the Roux et al. (2001) footwear survey and the two-parameter zero-inflated zeta (ZIZ) model. The ZIZ model is especially useful for illustrating parameter uncertainty because its inflation parameter, pi, and its zeta shape parameter can be displayed together as a two-dimensional region.

data("Psurveys")
roux = Psurveys$roux
roux
#> Number of Groups
#> 
#>   n    rn
#> ---  ----
#>   0   754
#>   1     9
#>   2     8
#>   3     4
#>   4     1
#> Roux C, Kirk R, Benson S, Van Haren T, Petterd C (2001).
#> "Glass particles in footwear of members of the public in
#> south-eastern Australia-a survey." _Forensic Science
#> International_, *116*(2), 149-156.
#> doi:10.1016/S0379-0738(00)00355-8
#> <https://doi.org/10.1016/S0379-0738%2800%2900355-8>.

Four approaches available in fitPS

The same probability model can be fitted using four different approaches to uncertainty.

Method Request What the uncertainty represents
Maximum likelihood fit(..., method = "mle") Frequentist uncertainty about fitted parameters, including profile-likelihood confidence regions
Ordinary bootstrap fit(..., method = "bootstrap") Repeated-sampling uncertainty estimated by resampling the observed survey and refitting the model
Parametric Bayes fit(..., method = "bayes") Posterior uncertainty about parameters and probabilities given the observed data and prior model
Rubin Bayesian Bootstrap fit(..., method = "bayesianBootstrap") Bayesian uncertainty in the empirical distribution, propagated through weighted refits

These methods do not attach the same interpretation to an 80% or 95% region. The purpose of comparing them is therefore not to decide which one has the “correct” contour. It is to see how different inferential treatments of the same survey affect the fitted parameters and, ultimately, the probabilities of interest.

Fit the ZIZ model

The public fit() interface keeps the model fixed while changing the inferential method.

mleFit = fit(
  roux,
  model = zizModel()
)
bootFit = fit(
  roux,
  model = zizModel(),
  method = "bootstrap",
  B = 1000,
  seed = 1234,
  silent = TRUE
)
bayesFit = fit(
  roux,
  model = zizModel(),
  method = "bayes"
)
bayesBoot = fit(
  roux,
  model = zizModel(),
  method = "bayesianBootstrap",
  B = 1000,
  seed = 1234
)

Confidence and credible regions for the parameters

For the two-parameter ZIZ model, separate intervals for pi and shape would miss an important feature: the two parameters can trade off against one another. A joint region shows which combinations of the parameters are supported together.

For a maximum-likelihood fit, confint() computes profile-likelihood confidence regions. For example, the following obtains 80% and 95% regions:

mleRegions = confint(
  mleFit,
  level = c(0.80, 0.95)
)
#> Computing contours. This may take a few seconds.

mleRegions contains the coordinates of the region boundaries so that they are available for further analysis if required. There is normally little value in printing hundreds of coordinate pairs. The region is much easier to understand graphically:

plotUncertainty(mleFit, level = c(0.80, 0.95))
Profile-likelihood confidence regions for the ZIZ inflation and shape parameters.
Profile-likelihood confidence regions for the ZIZ inflation and shape parameters.

For a Bayesian fit, the analogous calculation is a posterior credible region. credint() computes the requested region and plotUncertainty() displays it:

bayesRegions = credint(
  bayesFit,
  level = c(0.80, 0.95),
  silent = TRUE
)

A 95% confidence region and a 95% credible region are not the same probability statement. The confidence region belongs to a repeated-sampling procedure; the credible region describes posterior probability conditional on the observed data and prior model.

Comparing parameter uncertainty across the four methods

The four methods are most informative when shown on the same parameter scale.

plotUncertainty(mleFit, level = c(0.80, 0.95))
plotUncertainty(bootFit, level = c(0.80, 0.95))
plotUncertainty(bayesFit, level = c(0.80, 0.95))
plotUncertainty(bayesBoot, level = c(0.80, 0.95))
ZIZ parameter uncertainty under four inferential approaches. Top left: profile-likelihood confidence regions. Top right: ordinary-bootstrap regions with bootstrap realizations. Bottom left: parametric-Bayesian credible regions. Bottom right: Rubin Bayesian-Bootstrap regions with weighted-fit realizations.

ZIZ parameter uncertainty under four inferential approaches. Top left: profile-likelihood confidence regions. Top right: ordinary-bootstrap regions with bootstrap realizations. Bottom left: parametric-Bayesian credible regions. Bottom right: Rubin Bayesian-Bootstrap regions with weighted-fit realizations.

The plots have a common visual language, but their meanings differ. The MLE panel shows a profile-likelihood confidence region. The ordinary-bootstrap panel shows the sampling distribution of repeated fitted estimates. The parametric-Bayesian panel shows posterior probability for the parameters. The Bayesian-Bootstrap panel shows the uncertainty induced by changing the empirical weights assigned to the observed survey.

The bootstrap realization clouds are useful because they show the actual fitted values behind the smoothed regions. The shape and orientation of all four panels also make clear why reporting separate marginal intervals for pi and shape would lose information about their dependence.

From parameter uncertainty to probability uncertainty

The fitted P and S probabilities are functions of the model parameters. Uncertainty in pi and shape therefore propagates to uncertainty in the probabilities.

For maximum likelihood, the usual reported probabilities are the plug-in estimates

\[ P_k(\widehat{\theta}). \]

fitted(mleFit, n = 6, type = "plugIn")
#>           P0           P1           P2           P3           P4 
#> 0.9716490911 0.0169404164 0.0052597614 0.0022938450 0.0012050764 
#>           P5 
#> 0.0007122067

The ordinary bootstrap refits the model to resampled surveys and evaluates the probabilities for each refit. This produces both a bootstrap mean and percentile confidence intervals:

bootstrapProbs(bootFit, n = 6)
#>   term     estimate           sd        lower       upper level
#> 1   P0 0.9716375255 0.0057726138 0.9600514569 0.983246597  0.95
#> 2   P1 0.0170633473 0.0036934940 0.0099996106 0.024891696  0.95
#> 3   P2 0.0052063263 0.0010798241 0.0031081111 0.007331501  0.95
#> 4   P3 0.0022566563 0.0005132985 0.0012332454 0.003270360  0.95
#> 5   P4 0.0011838643 0.0002996450 0.0006146481 0.001782362  0.95
#> 6   P5 0.0007003115 0.0001948866 0.0003330060 0.001091135  0.95
#>   bootstrapMethod
#> 1   nonparametric
#> 2   nonparametric
#> 3   nonparametric
#> 4   nonparametric
#> 5   nonparametric
#> 6   nonparametric

For parametric Bayesian inference, each probability is averaged over the joint posterior distribution of the parameters. posteriorProbs() reports posterior means and equal-tailed credible intervals:

posteriorProbs(bayesFit, n = 6)
#>   term    estimate           sd        lower       upper level
#> 1   P0 0.970441723 0.0060790780 0.9583304216 0.981095900  0.95
#> 2   P1 0.018288884 0.0044455847 0.0102542494 0.028210446  0.95
#> 3   P2 0.005274782 0.0011371699 0.0033408441 0.007646046  0.95
#> 4   P3 0.002226757 0.0006084882 0.0012008310 0.003436022  0.95
#> 5   P4 0.001152871 0.0003918647 0.0004887586 0.001891575  0.95
#> 6   P5 0.000677721 0.0002719478 0.0002341282 0.001258605  0.95
#>   posteriorMethod
#> 1       numerical
#> 2       numerical
#> 3       numerical
#> 4       numerical
#> 5       numerical
#> 6       numerical

The posterior mean probability

\[ E\{P_k(\theta) \mid x\} \]

is not, in general, the same as evaluating the probability at the posterior mean parameter values. That distinction matters because P and S probabilities are nonlinear functions of the parameters.

Rubin’s Bayesian Bootstrap instead repeatedly changes the empirical weights on the observed survey and refits the model. Its stored probability summaries describe the resulting distribution of model-implied probabilities:

bayesBoot$probabilities[1:6, ]
#>   term     estimate           sd        lower       upper level
#> 1   P0 0.9716888212 0.0060212853 0.9588067841 0.982177349  0.95
#> 2   P1 0.0169814599 0.0037686732 0.0103735324 0.025265815  0.95
#> 3   P2 0.0052026003 0.0011254875 0.0032298046 0.007649017  0.95
#> 4   P3 0.0022611319 0.0005384930 0.0013305748 0.003437722  0.95
#> 5   P4 0.0011885914 0.0003150000 0.0006562026 0.001885156  0.95
#> 6   P5 0.0007042277 0.0002050671 0.0003654398 0.001161952  0.95

Comparing the probability estimates

Putting the four methods side by side shows how much the treatment of uncertainty changes the central probability estimates for the same ZIZ model.

probabilityComparison = data.frame(
  term = bootstrapProbs(bootFit, n = 6)$term,
  mlePlugIn = unname(fitted(mleFit, n = 6, type = "plugIn")),
  bootstrapMean = bootstrapProbs(bootFit, n = 6)$estimate,
  posteriorMean = posteriorProbs(bayesFit, n = 6)$estimate,
  bayesianBootstrapMean = bayesBoot$probabilities$estimate[1:6]
)

knitr::kable(
  probabilityComparison,
  digits = 5,
  col.names = c(
    "Term",
    "MLE plug-in",
    "Bootstrap mean",
    "Posterior mean",
    "Bayesian Bootstrap mean"
  )
)
Term MLE plug-in Bootstrap mean Posterior mean Bayesian Bootstrap mean
P0 0.97165 0.97164 0.97044 0.97169
P1 0.01694 0.01706 0.01829 0.01698
P2 0.00526 0.00521 0.00527 0.00520
P3 0.00229 0.00226 0.00223 0.00226
P4 0.00121 0.00118 0.00115 0.00119
P5 0.00071 0.00070 0.00068 0.00070

This comparison is useful because it separates two effects. The first is the fitted model itself: all four analyses use ZIZ. The second is how uncertainty is represented and averaged. Similar values indicate that the central fitted probabilities are relatively insensitive to that choice for this survey; visible differences identify terms for which the uncertainty treatment matters more.

The interval estimates tell the complementary story. The ordinary bootstrap gives percentile confidence intervals, the parametric Bayesian fit gives posterior credible intervals, and the Bayesian Bootstrap gives intervals induced by uncertainty in the empirical weights. The MLE plug-in estimate is included as a reference point.

Uncertainty in the first six fitted P probabilities. The MLE plug-in value is a point estimate; the other three methods show their stored 95% uncertainty intervals.
Uncertainty in the first six fitted P probabilities. The MLE plug-in value is a point estimate; the other three methods show their stored 95% uncertainty intervals.

The figure should be read by comparing methods within each probability term, rather than by comparing interval widths as though they had identical interpretations. In particular, a bootstrap confidence interval and a Bayesian credible interval answer different inferential questions even when their numerical endpoints are close.

The same probability summaries are available directly through the fitted-object prediction interface:

predict(mleFit, newdata = 0:5, type = "plugIn")
#>           P0           P1           P2           P3           P4 
#> 0.9716490911 0.0169404164 0.0052597614 0.0022938450 0.0012050764 
#>           P5 
#> 0.0007122067
predict(
  bootFit,
  newdata = 0:5,
  type = "bootstrapMean",
  interval = "percentile"
)
#>       predicted        lower       upper
#> P0 0.9716375255 0.9600514569 0.983246597
#> P1 0.0170633473 0.0099996106 0.024891696
#> P2 0.0052063263 0.0031081111 0.007331501
#> P3 0.0022566563 0.0012332454 0.003270360
#> P4 0.0011838643 0.0006146481 0.001782362
#> P5 0.0007003115 0.0003330060 0.001091135
predict(
  bayesFit,
  newdata = 0:5,
  type = "posteriorMean",
  interval = "credible"
)
#>      predicted        lower       upper
#> P0 0.970441723 0.9583304216 0.981095900
#> P1 0.018288884 0.0102542494 0.028210446
#> P2 0.005274782 0.0033408441 0.007646046
#> P3 0.002226757 0.0012008310 0.003436022
#> P4 0.001152871 0.0004887586 0.001891575
#> P5 0.000677721 0.0002341282 0.001258605

What does the inflation parameter mean?

For the ZIZ model, pi is more than a technical fitting parameter. It measures extra probability at the inflated value beyond that supplied by the zeta component. For a P survey this is extra mass at P0; for an S survey it is extra mass at S1.

A Bayesian analysis can therefore ask whether the inflation is practically negligible. If epsilon is an application-specific threshold, posteriorInflation() calculates

\[ \Pr(\pi < \epsilon \mid x). \]

posteriorInflation(bayesFit, epsilon = 0.01)
#>   epsilon    probBelow probAtOrAbove posteriorMethod
#> 1    0.01 2.998125e-05       0.99997       numerical

This is often more informative than asking whether pi is exactly zero. Under a continuous prior, the posterior probability of one exact parameter value is zero; a practically meaningful threshold asks whether the inflation is small enough not to matter in the application.

Which method should I use?

The choice should follow the inferential question.

Whichever method is used, report its name together with the interval or region. A percentage such as “95%” is not enough on its own because confidence regions, percentile bootstrap intervals, and posterior credible intervals do not have the same interpretation.

For two-parameter models, plotUncertainty() is particularly useful because it shows parameter dependence that separate intervals cannot. For applied work, however, the final step should usually be to examine the uncertainty in the P or S probabilities themselves. Those are the quantities through which parameter uncertainty affects the substantive calculation.

References

Curran, J. M., Buzzini, P., and Trejos, T. (2024). Estimating probability terms for the background presence of glass when considering activity in forensic casework. Forensic Science International, 364, 112221. https://doi.org/10.1016/j.forsciint.2024.112221

Efron, B. (1979). Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1), 1-26.

Rubin, D. B. (1981). The Bayesian bootstrap. The Annals of Statistics, 9(1), 130-134.

Roux, C., Kirk, R., Benson, S., Van Haren, T., and Petterd, C. I. (2001). Glass particles in footwear of members of the public in south-eastern Australia: a survey. Forensic Science International, 116(2), 149-156.