Why Reproducibility Requires More Than an Estimate
A selection-model likelihood can depend on both an outcome model and a model for whether outcomes are observed. When an outcome is missing, its contribution to the likelihood requires integrating over possible missing values. The reported estimate may therefore depend not only on the stated statistical model but also on parameterization, numerical integration, optimization choices, and the method used to quantify uncertainty. A coefficient and standard error alone do not tell another analyst enough to reconstruct or evaluate the analysis.
The earlier material on numerical integration, optimization, and likelihood-based inference establishes that these computational decisions can matter especially near ridges, boundaries, or competing maxima. The purpose here is to report those decisions and their evidence in a concise, auditable form—not to substitute a list of software settings for a statistical argument.
State the Likelihood Being Maximized
Let \(Y_i=(Y_i^{\mathrm{obs}},Y_i^{\mathrm{mis}})\) denote the complete outcome for subject \(i\), partitioned into observed and missing components, and let \(R_i\) indicate the observation pattern. Let \(x_i\) be the recorded covariates. A selection model specifies an outcome density \(f_\theta(y_i\mid x_i)\) and a pattern probability \(g_\gamma(r_i\mid y_i,x_i)\). The parameter vector is \(\psi=(\theta,\gamma)\). For continuous missing outcomes, the subject's observed-data likelihood contribution is the integral below; for discrete missing outcomes it is the corresponding sum.
For a fully observed outcome, there is no missing component to integrate over, so the contribution is \(f_\theta(y_i\mid x_i)g_\gamma(r_i\mid y_i,x_i)\). In a report, specify whether the displayed likelihood includes factors that are constant in \(\psi\), such as a multinomial coefficient. Such factors do not change the maximizer, but stating the convention prevents confusion when comparing reported log-likelihood values. Also identify the analysis sample, the outcome and covariates, the coding of the observation pattern, and any assumptions about independence across subjects.
Worked Example: Writing a Binary-Outcome Likelihood
Suppose a hypothetical sample contains 120 independent subjects. A binary outcome \(Y\) is observed for 90 subjects: 36 have \(Y=1\) and 54 have \(Y=0\). The outcome is missing for 30 subjects. Use the model \(P(Y=1)=p=\operatorname{expit}(\eta)\) and \(P(R=1\mid Y=y)=\pi_y=\operatorname{expit}(\alpha+\beta y)\), where \(R=1\) means \(Y\) is observed. Assume these probabilities describe the sampling process and that subjects are independent.
An observed \(Y=1\) contributes \(p\pi_1\), an observed \(Y=0\) contributes \((1-p)\pi_0\), and a missing outcome contributes the sum of its two possible complete-data contributions: \(p(1-\pi_1)+(1-p)(1-\pi_0)\). Thus, the likelihood kernel is
The full likelihood for these three observed cell counts includes the multinomial coefficient \(120!/(36!\,54!\,30!)\), which does not depend on \((\eta,\alpha,\beta)\). The log-likelihood kernel is
This expression makes the missing-outcome integration explicit: the 30 missing cases contribute the probability of missingness after summing over \(Y=0\) and \(Y=1\), not an invented outcome value. For reproducibility, a report would additionally state the covariate structure if present, the link functions, the software implementation, and the numerical and inferential procedures used to fit this likelihood.
Report the Computation as Part of the Analysis
A likelihood formula does not uniquely describe a computational result when its terms require numerical integration or numerical optimization. A useful report names the integration method, nodes or order, transformations and tolerances, and any adaptive-integration controls that affect evaluation. If estimates or likelihood values are sensitive to these settings, show the comparison used to judge adequacy rather than saying only that integration was “accurate.”
For optimization, record the parameterization actually optimized, including transformations or rescaling, the algorithm, tolerances, and how starting values were selected. Explain the convergence checks: for example, the optimizer's termination status, the scaled gradient norm, and whether the fitted point is a local maximum according to a curvature check. A small gradient alone is not proof of a global maximum. If multiple starts were used, report how many, how they were generated, how many reached the best objective value, and what happened to unsuccessful fits. Distinguish a numerical warning from a failed replicate, and never omit a fit without explaining the rule.
Worked Example: Reporting Integration and Optimization Checks
Consider a hypothetical continuous-outcome selection model whose observed likelihood requires one-dimensional quadrature for each subject. The analyst evaluates the same likelihood with Gauss–Hermite quadrature at orders 15, 25, and 35 and uses the 35-node result for the final fit. The following hypothetical summary records the estimate of a treatment contrast \(\Delta\) and the maximized log likelihood.
| Quadrature order | \(\hat{\Delta}\) (outcome units) | Maximized log likelihood |
|---|---|---|
| 15 | 0.0412 | -182.417 |
| 25 | 0.0410 | -182.413 |
| 35 | 0.0410 | -182.413 |
The comparison supports using order 35 for this fit: increasing the order from 25 to 35 leaves the displayed contrast unchanged to the shown precision and changes the log likelihood by less than 0.001. It does not establish accuracy for every parameter or every dataset; a report should also say which quantities were compared and whether the quadrature convergence criterion was prespecified.
Suppose optimization used transformed parameters, 40 dispersed starting values, and a scaled-gradient stopping tolerance of \(10^{-6}\). The best solution was reached from 12 starts; the maximum scaled-gradient component at that solution was \(3.2\times10^{-7}\), and the observed information was positive definite under the stated sign convention for a local maximum. The other starts either reached lower log likelihoods or stopped with warnings. A reproducible account reports these facts, the software and version, the optimization algorithm and tolerances, and how starts were constructed. It does not claim that agreement among starts proves a global optimum; it provides evidence that can be checked and extended.
Justifying the Observed-Data Contribution
The integration in the likelihood follows from the law of total probability. For any observed pattern \(r_i\) and observed value \(y_i^{\mathrm{obs}}\), the joint density or probability mass of what was recorded is obtained by summing or integrating the complete-data joint distribution over the unrecorded values. Densities and masses are understood with respect to the appropriate dominating measures. Under the stated model, that complete-data joint distribution factors into an outcome component and a selection component.
The second equality uses the selection-model factorization: conditional on the complete outcome and covariates, the pattern probability is \(g_\gamma\), while the complete outcome distribution is \(f_\theta\). It is this model-based factorization that makes the missing values contribute through integration. The report should identify it, along with the assumed outcome and selection distributions. If the data are clustered or repeated within subjects, the independence statement and likelihood unit must instead reflect the specified dependence structure.
Make Uncertainty Procedures Reconstructible
State the inferential target in context, such as a treatment contrast in outcome units, and explain how its uncertainty was obtained. For a likelihood-ratio interval, specify the likelihood quantity, cutoff, and any calibration used. For a parametric bootstrap, state the fitted model used to generate data, the number of datasets simulated, the refitting procedure, the interval rule, and the random-number details needed to reproduce the simulation. For every method, report failed fits and the rule for handling them.
Worked Example: Reporting a Parametric Bootstrap Without Mixing Counts
Suppose the fitted selection model gives a treatment contrast estimate of 0.041 outcome units. A hypothetical analysis simulates 2,000 bootstrap datasets from the fitted model and refits the model once to each dataset. The planned interval is the percentile interval using the 2.5th and 97.5th percentiles of the final replicate estimates.
Of the 2,000 initial fits, 37 trigger a prespecified numerical warning. The analyst refits those same 37 datasets using the same model, dispersed starting values, and tighter tolerances. All 37 then satisfy the stated convergence checks. There are therefore 2,000 simulated datasets and 2,000 final bootstrap estimates, but 2,037 fitting attempts: 2,000 initial attempts plus 37 additional attempts. Refitting does not create new datasets or extra independent bootstrap estimates.
Suppose the 2.5th and 97.5th percentiles of the 2,000 final estimates are 0.008 and 0.076. The reported 95% percentile bootstrap interval is 0.008 to 0.076 outcome units. Under repeated simulation and refitting using this procedure, the percentile rule is intended to provide an interval with approximately 95% coverage under its assumptions; it is not a statement that the fixed parameter has a 95% probability of lying in this realized interval. A reproducible report also gives the software and version, random-number generator and seed, simulation model, fit tolerances, warning and refit rules, and the number of final estimates used. If any fit remained unsuccessful, the report would state how it was handled and the resulting number of usable estimates rather than silently dropping it.
Diagnostics, Alternatives, and Limits
Warning signs include estimates or log likelihoods that change materially with quadrature order, large gradients at the reported solution, materially different solutions from different starts, boundary estimates, a nearly singular curvature matrix, and unexplained differences between software runs. Such evidence should prompt investigation and transparent reporting, not a cosmetic claim of convergence. Report parameter estimates with units and enough precision to distinguish relevant numerical changes; include the likelihood value and uncertainty method where they help audit the analysis.
There is no single reporting format that replaces the inferential method. A likelihood-ratio interval uses a likelihood cutoff and may require nonstandard calibration; a Wald interval uses local curvature and can be unreliable near ridges or boundaries; a parametric bootstrap approximates sampling behavior by repeated simulation and refitting. These methods answer uncertainty questions through different approximations. Name the method actually used and justify it for the fitted situation. Reporting several methods can be informative when they are clearly distinguished, but they should not be presented as interchangeable confirmations.
Reproducibility also has limits. A fully specified computational record makes the analysis auditable, but it cannot establish that the selection model is scientifically correct, that parameters are identified by the available data, or that conclusions are robust to unexamined assumptions. A precise report should separate numerical stability from substantive credibility and should avoid implying that a successful optimizer resolves weak identification.
Name the likelihood unit, observed outcomes, missingness pattern, covariates, outcome model, selection model, and parameterization.
Give integration method and settings, software, optimizer, parameter scaling, tolerances, and starting-value strategy.
Report convergence checks, multiple-start results, sensitivity checks, the inferential target, uncertainty method, and all simulation and refitting counts.
Common Reporting Mistakes
- Giving only a coefficient and standard error. That does not identify the likelihood, parameterization, numerical approximation, or method behind the uncertainty estimate. State these explicitly.
- Calling an optimizer's success flag proof of a solution. A termination message is one diagnostic, not evidence by itself of a local maximum or global optimum. Pair it with gradient and curvature checks, and report multiple-start evidence where used.
- Equating fitting attempts with bootstrap replicates. A retry on an already simulated dataset is another fitting attempt, not another independent replicate. Report dataset count, final estimate count, and total attempts separately.
- Claiming numerical accuracy from a single setting. State the sensitivity comparison or accuracy checks supporting the selected integration settings; do not imply that an untested approximation is exact.
Check Your Understanding
Answer each question, then compare with the solutions.
- Why does a subject with a missing outcome contribute an integral or sum to a selection-model likelihood?
- A report says it simulated 1,500 bootstrap datasets and retried 20 fits on those same datasets. How many datasets and fitting attempts should it report if every initial fit was attempted?
- Name two distinct pieces of evidence, beyond an optimizer's success message, that can support a reported solution.
- Why is a comparison of quadrature orders useful, and what does it not establish?
- What key distinction should a report make between a stable computation and a credible scientific conclusion?
Solutions. (1) The missing value is unobserved, so its possible values must be summed or integrated out of the complete-data joint distribution, including both outcome and selection probabilities. (2) It should report 1,500 datasets and, if all 1,500 initial fits were attempted, 1,520 fitting attempts; if all retries succeed, there are 1,500 final replicate estimates. (3) Examples include a small scaled gradient, a curvature check consistent with a local maximum, and agreement of the best solution across dispersed starts. None alone proves a global optimum. (4) It assesses whether estimates or likelihood values materially change as the approximation is refined; it does not prove accuracy for every quantity, dataset, or model. (5) Numerical stability concerns computation under the specified model, while scientific credibility also depends on the plausibility and adequacy of that model and its assumptions.