Why Reproducible Computation Matters
A Bayesian result is produced by a chain of computational choices: which data enter the analysis, how those data are represented, what model and prior are used, and how posterior quantities are calculated. A reported posterior probability or interval is difficult to review if the submission does not preserve enough information to reconstruct that chain. Reproducible computation makes the chain inspectable: an appropriately equipped reviewer can run the documented analysis and obtain the same result, or understand a justified numerical difference.
Reproducibility is not the same as correctness. Re-running code can show that a result follows from the supplied data, model, and implementation, but it cannot by itself show that the data are accurate, the model is appropriate, or the clinical interpretation is justified. Software validation is a related but distinct activity; the earlier material on validation of Bayesian analysis software addresses whether software performs as intended. Here the focus is the documentation of a particular analysis and its computational execution.
What a Reproduction Record Must Specify
A useful record is more than a copy of the final script. It should make explicit what the script reads, what it does, and how its output is judged. The following items are typically needed for key analyses in a regulatory submission:
- Data specification: the precise analysis dataset or an authorized route to it; variable names, formats, units, permitted values, and coding conventions; the analysis population; and all inclusion, exclusion, derivation, and transformation rules.
- Model specification: the likelihood, prior distributions, parameterization, links or transformations, and the quantities reported. Record enough detail to distinguish, for example, a prior on a probability from a prior on its log-odds.
- Code and dependencies: the analysis code, any functions or macros it calls, the order of execution, and the exact versions of software libraries that affect the calculation.
- Computational environment: the software and operating-system versions, relevant numerical libraries, and any settings that can affect results. Where practical, preserve a locked dependency file or a controlled environment image as well as a human-readable version list.
- Computational settings: for simulation-based calculations, document the algorithm, initialization, number of chains, warm-up and retained iterations, thinning if used, random-number generator and seeds, and convergence and precision settings. For deterministic calculations, document numerical tolerances, integration settings, or optimizer controls as applicable.
- Expected outputs and checks: identify the key results to compare, their units and rounding rules, acceptable numerical tolerances, and checks that confirm the data and analysis ran as intended.
- Provenance: preserve a version identifier for the code and inputs, execution date, and a record of changes. A cryptographic checksum can help verify that a file has not changed; it does not establish that the file contains the correct data.
The specification should be concrete enough to prevent a reviewer from having to infer an important choice. “Use the final dataset” is not sufficient if several files could qualify. “Use the locked analysis extract identified by its filename, version, and checksum” is more precise. Similarly, “run the model with defaults” is not reproducible when defaults vary by software version.
Worked Example: Reconstructing a Conjugate Posterior
Consider a hypothetical trial analysis with 40 participants and a binary outcome. There are 28 participants with the outcome and 12 without it. The model treats the outcome indicator for each participant as Bernoulli with a common probability \(p\), and assigns \(p\) a Beta(1, 1) prior. Assume the participants are conditionally independent given \(p\), the outcome is coded consistently as 1 for an event and 0 otherwise, and all 40 outcomes are included.
For \(x\) events among \(n\) participants, the likelihood is proportional to \(p^x(1-p)^{n-x}\). Combining it with a Beta(\(a,b\)) prior gives a Beta(\(a+x,b+n-x\)) posterior. Here:
The posterior mean is \(29/(29+13)=29/42=0.6905\). Thus, under this model and prior, the posterior mean event probability is 0.6905, or 69.05%. That numerical result is reproducible only if the analysis record fixes the 40-person analysis population, the event coding, the prior, the model, and the calculation.
A concise computational record would state the input counts (28 events, 12 non-events), the likelihood and prior above, the software and version, the code or calculation used, and the expected posterior parameters and mean. A reviewer can check the arithmetic independently and compare it with the submitted output. The record should also preserve the participant-level analysis data or an appropriately controlled means of accessing them; the aggregate counts alone cannot reveal whether the analysis population was constructed as specified.
Data and Execution Checks
Before comparing posterior results, verify that the same analysis was actually run. Data checks should establish that the submitted file is the intended file and that its contents agree with the specification. Useful checks include record counts, unique participant identifiers, allowable values, missingness summaries where relevant, and reconciliation of derived totals. These checks complement—not replace—appropriate data-quality processes.
Execution checks should record whether the run completed as intended and whether the outputs contain the expected parameters, summaries, and units. For stochastic methods, retain diagnostic summaries and relevant computational settings alongside the key estimates. If output differs from the expected result, investigate the input data, transformations, code version, environment, and numerical behavior rather than silently adjusting the reported value.
Worked Example: Documenting an Analysis Population
Suppose a hypothetical binary-outcome analysis begins with 10 screened participants. The documented rule includes participants who received at least one dose and met the prespecified eligibility criteria. One participant received no dose, and one did not meet the eligibility criteria; both are excluded. The remaining eight participants all have an observed outcome, of whom six are event-positive and two are event-negative.
The population reconciliation is \(10-1-1=8\) included participants. The outcome reconciliation is \(6+2=8\), so the included records are accounted for by the two outcome categories. If the analysis uses a Beta(1, 1) prior and the Bernoulli model described above, the posterior is:
Its posterior mean is \(7/(7+3)=0.7000\). A reproducible data specification would preserve the eligibility rule, the reason for each exclusion, the event coding, and the resulting count checks. A reviewer should be able to confirm both the eight-person analysis set and the six-event, two-non-event split, rather than being asked to trust only the final posterior output.
If the input file instead contained nine included records, or if the event totals summed to seven, the mismatch would be a failed check. The analysis should not proceed to an unqualified report until the discrepancy is resolved and its resolution recorded.
Computational Settings and Numerical Agreement
Not every result can be reproduced by a single closed-form calculation. Numerical integration, optimization, and posterior simulation can introduce small differences due to algorithm settings, finite precision, or random variation. The submission should define what counts as agreement before comparing runs. Exact agreement is appropriate for quantities that should be identical under the same deterministic calculation and environment. When exact agreement is not expected, specify a tolerance appropriate to the quantity and its precision, and explain why that tolerance is adequate for the clinical decision.
For independent Monte Carlo draws, the estimated probability of an event is a sample proportion. If \(I_j\) is 1 when draw \(j\) satisfies the event and 0 otherwise, then the estimate is \(\hat{q}=\sum_{j=1}^{N} I_j/N\). Under independent, identically distributed draws, its Monte Carlo standard error is approximately \(\sqrt{\hat{q}(1-\hat{q})/N}\). This quantifies simulation error, not uncertainty about the clinical parameter itself. For correlated MCMC draws, this independent-draw expression is not appropriate with the raw iteration count; autocorrelation must be accounted for, for example through an effective sample size and an appropriate Monte Carlo error estimate.
Worked Example: Reporting Monte Carlo Precision
In a hypothetical simulation-based calculation, 10,000 retained independent draws are used to estimate a posterior event probability. The event occurs in 7,420 draws, so the estimated probability is:
Under the stated independent-draw assumption, the Monte Carlo standard error is:
The simulation estimate is therefore 0.7420, with Monte Carlo standard error about 0.0044. A reproducible report would identify the event being counted, the number of retained draws, the algorithm and its settings, the software and random-number configuration, and the resulting event count. The calculation does not imply that the underlying posterior event probability is exactly 0.7420; it reports a finite-simulation estimate.
Justifying Reproduction Checks
A comparison rule should match the computational method. The justification for a deterministic calculation is different from the justification for a stochastic estimate. In the deterministic case, if the input, algorithm, and arithmetic are fixed, identical outputs are expected apart from known differences such as floating-point rounding. In a simulation, finite draws add Monte Carlo variation, so a difference between two estimates must be interpreted relative to their simulation error.
Name the specific result, its definition, units, and displayed precision. Do not compare rounded values if the unrounded output is available.
For a deterministic result, identify whether exact agreement is expected. For a simulation estimate, record the estimator and an appropriate measure of Monte Carlo error.
Choose a tolerance before reviewing the reproduction result, justify it in relation to numerical precision or simulation error, and report whether the comparison passes.
Check the data version, transformations, code and dependency versions, seeds and generator, algorithm settings, and convergence or numerical diagnostics. Document the cause and any corrective action.
For example, the two Monte Carlo estimates should not be expected to match digit for digit if independent random streams are used. By contrast, an unexpected change in the event definition or analysis population is not a Monte Carlo fluctuation. Separating computational variability from changes in the analysis is essential to a meaningful reproduction assessment.
Alternative Computational Approaches and Their Limits
When available, a closed-form result can provide a transparent reference calculation, as in the Beta-binomial example. Deterministic numerical integration or optimization can also avoid Monte Carlo sampling error for some quantities, but requires documentation of the numerical method, tolerances, and convergence behavior. Posterior simulation is flexible for complex models, but its results depend on algorithm settings and require assessment of simulation error and chain behavior. These approaches are not interchangeable in every model; the selected method should be justified by the analysis and documented well enough to reproduce its output.
A fixed random seed can help reproduce a simulation in the same software environment, but it is not a complete reproducibility plan. Results may change with a different random-number generator, software version, parallel execution pattern, or numerical library. Conversely, a seed does not ensure good mixing or adequate precision. Record the seed and generator where applicable, but also report the method, settings, diagnostics, and precision relevant to the result.
Common Mistakes and Failure Modes
- Providing only the final code. Code that reads an unspecified file or relies on undocumented defaults does not identify the analysis inputs. Preserve the data specification, dependencies, settings, and run instructions as well.
- Confusing a checksum with data verification. A matching checksum supports the claim that two files are byte-for-byte identical. It does not show that the file is the right extract or that its records are accurate.
- Treating a seed as a guarantee. A seed helps define a random stream in a specified environment; it does not guarantee identical results across environments or adequate Monte Carlo precision.
- Using agreement as proof of validity. An incorrect model can be implemented and reproduced consistently. Reproducibility supports transparency and checking, not the scientific adequacy of every assumption.
- Ignoring small discrepancies. A difference should not be dismissed as “rounding” without checking its size, source, and relevance to the reported conclusion.
The practical limitation is that exact reproduction may be constrained by access controls, software licensing, changing infrastructure, or the use of stochastic algorithms. In such cases, document the constraint and provide the most complete authorized record possible. State which elements can be independently rerun and which results are compared using a justified numerical tolerance.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A script produces the same posterior mean on two computers, but each computer reads a file named “final.csv.” What additional information would make the data input reproducible?
- In a simulation with 10,000 independent draws, an event occurs 7,420 times. What is the estimated event probability and its approximate Monte Carlo standard error?
- Why does recording a random seed not, by itself, establish that a simulation result is reliable?
- A checksum matches for the submitted dataset. What does that establish, and what does it not establish?
Solutions. (1) Identify the exact dataset version, its controlled location or authorized access route, a checksum, and the data schema and transformations used. The filename alone may refer to different files. (2) The estimate is \(7{,}420/10{,}000=0.7420\). Its standard error is \(\sqrt{0.7420(0.2580)/10{,}000}\approx0.004375\), or about 0.0044. (3) The seed does not establish adequate mixing, sufficient simulation precision, or identical behavior across software environments; the generator, algorithm, settings, and diagnostics also matter. (4) A matching checksum supports that the compared files are byte-for-byte identical. It does not establish that the file is the correct analysis extract or that its data are accurate.