--- title: "Quantifying uncertainty with fitPS" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Quantifying uncertainty with fitPS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5 ) suppressPackageStartupMessages(library(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. ```{r roux-data} data("Psurveys") roux = Psurveys$roux roux ``` ## 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. ```{r fit-mle} mleFit = fit( roux, model = zizModel() ) ``` ```{r fit-bootstrap, eval=FALSE} bootFit = fit( roux, model = zizModel(), method = "bootstrap", B = 1000, seed = 1234, silent = TRUE ) ``` ```{r fit-bayes} bayesFit = fit( roux, model = zizModel(), method = "bayes" ) ``` ```{r fit-bayesian-bootstrap, eval=FALSE} bayesBoot = fit( roux, model = zizModel(), method = "bayesianBootstrap", B = 1000, seed = 1234 ) ``` ```{r load-bootstrap-results, include=FALSE} uncertaintyCache = readRDS(file.path("figures", "uncertainty-bootstrap-cache.rds")) bootFit = uncertaintyCache$bootFit bayesBoot = uncertaintyCache$bayesBoot ``` ## 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: ```{r mle-regions, results='hide'} mleRegions = confint( mleFit, level = c(0.80, 0.95) ) ``` `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: ```{r mle-region-plot, fig.cap="Profile-likelihood confidence regions for the ZIZ inflation and shape parameters."} plotUncertainty(mleFit, level = c(0.80, 0.95)) ``` For a Bayesian fit, the analogous calculation is a posterior credible region. `credint()` computes the requested region and `plotUncertainty()` displays it: ```{r bayes-regions, results='hide'} 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. ```{r uncertainty-region-comparison-code, eval=FALSE} 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)) ``` ```{r uncertainty-region-comparison, echo=FALSE, out.width="100%", fig.cap="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."} knitr::include_graphics(file.path("figures", "uncertainty-comparison.png")) ``` 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}). \] ```{r mle-probabilities} fitted(mleFit, n = 6, type = "plugIn") ``` 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: ```{r bootstrap-probabilities} bootstrapProbs(bootFit, n = 6) ``` 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: ```{r posterior-probabilities} posteriorProbs(bayesFit, n = 6) ``` 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: ```{r bayesian-bootstrap-probabilities} bayesBoot$probabilities[1:6, ] ``` ## 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. ```{r compare-four-probability-estimates} 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" ) ) ``` 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. ```{r probability-uncertainty-comparison, echo=FALSE, fig.height=5, fig.cap="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."} bootProb = bootstrapProbs(bootFit, n = 6) posteriorProb = posteriorProbs(bayesFit, n = 6) bayesBootProb = bayesBoot$probabilities[1:6, ] terms = bootProb$term methodNames = c( "MLE plug-in", "Ordinary bootstrap", "Parametric Bayes", "Bayesian Bootstrap" ) methodOffset = c(-0.27, -0.09, 0.09, 0.27) termPosition = seq_along(terms) allLower = c( fitted(mleFit, n = 6, type = "plugIn"), bootProb$lower, posteriorProb$lower, bayesBootProb$lower ) allUpper = c( fitted(mleFit, n = 6, type = "plugIn"), bootProb$upper, posteriorProb$upper, bayesBootProb$upper ) plot( NA, xlim = range(c(allLower, allUpper), finite = TRUE), ylim = c(0.5, length(terms) + 0.5), yaxt = "n", ylab = "", xlab = "Fitted probability" ) axis(2, at = termPosition, labels = terms, las = 1) mleValues = unname(fitted(mleFit, n = 6, type = "plugIn")) points(mleValues, termPosition + methodOffset[1], pch = 16) intervalSets = list(bootProb, posteriorProb, bayesBootProb) for (methodIndex in seq_along(intervalSets)) { intervalData = intervalSets[[methodIndex]] y = termPosition + methodOffset[methodIndex + 1L] segments(intervalData$lower, y, intervalData$upper, y) points(intervalData$estimate, y, pch = methodIndex) } legend( "topright", legend = methodNames, pch = c(16, 1, 2, 3), bty = "n", cex = 0.8 ) ``` 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: ```{r probability-prediction-interface} predict(mleFit, newdata = 0:5, type = "plugIn") predict( bootFit, newdata = 0:5, type = "bootstrapMean", interval = "percentile" ) predict( bayesFit, newdata = 0:5, type = "posteriorMean", interval = "credible" ) ``` ## 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). \] ```{r practical-inflation} posteriorInflation(bayesFit, epsilon = 0.01) ``` 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. - Use **maximum likelihood and profile likelihood** when a frequentist analysis is required and joint parameter uncertainty is important. - Use the **ordinary bootstrap** when you want sampling uncertainty propagated directly through the fitting procedure to the derived probabilities. - Use **parametric Bayes** when posterior inference is appropriate and the prior model is scientifically defensible. - Use **Rubin's Bayesian Bootstrap** when the Bayesian uncertainty of interest is attached primarily to the empirical distribution rather than to a parametric prior on the model parameters. 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.