Bayesian workflow · added in 2026

Can this posterior be trusted?

The 2023 assignment ran one chain per sampler and stopped at a comparison plot. This page adds the checks a careful analysis needs: four chains per sampler with R-hat and effective sample sizes, posterior predictive checks against the observed table, a prior sensitivity analysis, and a comparison of all four methods with the exact posterior. Seeds are fixed and shown, and the convergence diagnostics are verified against R’s posterior package on the same draws.

1 · Convergence

Four chains per sampler, from dispersed starting points

Each sampler is the unmodified port of the 2023 function, run four times with set.seed(90125) to set.seed(90128). Chain 1 is exactly the chain the original analysis ran. Each chain starts from its own rnorm(p) draw, as the original code does, which is wide relative to the posterior (the slope’s posterior sd is about 0.13).

R-hat compares the spread between chains with the spread within them, and values near 1 mean the chains agree. The version here is the rank-normalised, split, folded R-hat of Vehtari et al. (2021), which also catches chains that agree on the centre but not the tails. Bulk and tail ESS count how many independent draws the correlated chains are worth for the centre and the 5% and 95% quantiles. Every posterior mean below carries its Monte Carlo standard error, and digits beyond it are not printed.

Thresholds: R-hat at most 1.01 and both ESS at least 400 (Vehtari et al. 2021). For the corrected model every coefficient passes with room to spare: the largest R-hat is 1.0002, the smallest bulk ESS 26,504 and the smallest tail ESS 29,913.

Chains

Seed, starting values and acceptance rate of each chain
Chain	set.seed()	start (β₀, β₁)	accepted: MH	HMC
12023	90125	(0.39, -0.16)	35.8%	95.9%
2	90126	(-1.05, -0.10)	35.9%	96.2%
3	90127	(0.43, 0.44)	35.2%	96.4%
4	90128	(0.36, 0.15)	36.1%	96.3%

Metropolis-Hastings 4 chains × 50,000 draws after warm-up · 100,000 iterations, first 50,000 discarded, proposal N(θ, c²·vcov(GLM)) with c = 2.4/√2

Metropolis-Hastings: posterior summaries and convergence diagnostics over 4 chains of 50,000 post-warm-up draws
Coefficient	mean ± MCSE	sd	95% interval (± MCSE)	R-hat	bulk ESS	tail ESS	Status
β₀ (Intercept)	-0.8159 ± 0.0028	0.464	[-1.758 ±0.008, 0.072 ±0.007]	1.0002	26,643	33,115	passes
β₁ dose	0.31536 ± 0.00081	0.132	[0.065 ±0.001, 0.585 ±0.002]	1.0001	26,504	33,472	passes

Hamiltonian Monte Carlo 4 chains × 9,000 draws after warm-up · 10,000 iterations, first 1,000 discarded, L = 20 leapfrog steps of ε = 1/20, identity mass matrix

Hamiltonian Monte Carlo: posterior summaries and convergence diagnostics over 4 chains of 9,000 post-warm-up draws
Coefficient	mean ± MCSE	sd	95% interval (± MCSE)	R-hat	bulk ESS	tail ESS	Status
β₀ (Intercept)	-0.8157 ± 0.0015	0.461	[-1.753 ±0.009, 0.068 ±0.008]	1.0002	94,387	29,913	passes
β₁ dose	0.31491 ± 0.00045	0.133	[0.062 ±0.002, 0.582 ±0.002]	1.0000	87,223	31,864	passes

HMC’s bulk ESS exceeds its number of draws (36,000): successive draws are negatively correlated, so their average is more precise than independent draws would give. The tail ESS, which does not benefit, is the more conservative number.

Chain 1 seed 90125
Chain 2 seed 90126
Chain 3 seed 90127
Chain 4 seed 90128

Metropolis-Hastingsβ₁ dose

Warm-up: iterations 1 to 200

After warm-up: iterations 50,001 to 100,000 (1 in every 167 draws plotted)

Rank plot after warm-up (20 bins; dashed line = what perfect mixing gives)

Chain 1

Chain 2

Chain 3

Chain 4

Hamiltonian Monte Carloβ₁ dose

Warm-up: iterations 1 to 200

After warm-up: iterations 1,001 to 10,000 (1 in every 30 draws plotted)

Rank plot after warm-up (20 bins; dashed line = what perfect mixing gives)

Chain 1

Chain 2

Chain 3

Chain 4

Optional · bring your own key

Explain these diagnostics

A language model can turn the tables above into plain English. It sees only the rounded numeric summary (open “What will be sent”), its answer is labelled AI-generated, checked automatically and logged with your decision in your browser’s AI audit log. Nothing on this page depends on it. See the AI use statement.

What will be sent (this JSON and fixed instructions, nothing else, about 2,837 characters)

{
  "summary_id": "corrected",
  "analysis": "Bayesian logistic regression, logit P(improvement) = b0 + b1 * dose, flat prior, fitted by four MCMC chains of each sampler.",
  "data": "Dose-finding study: 7 doses (0 to 6), 10 patients per dose, 37 of 70 improved.",
  "thresholds": {
    "rhat_max": 1.01,
    "ess_bulk_min": 400,
    "ess_tail_min": 400
  },
  "samplers": [
    {
      "sampler": "Metropolis-Hastings",
      "chains": 4,
      "seeds": [
        90125,
        90126,
        90127,
        90128
      ],
      "draws_per_chain": 50000,
      "warmup_discarded_per_chain": 50000,
      "acceptance_rate_by_chain": [
        0.358,
        0.359,
        0.352,
        0.361
      ],
      "parameters": [
        {
          "name": "β₀ (Intercept)",
          "mean": -0.8159,
          "sd": 0.4637,
          "mcse_mean": 0.0028,
          "rhat": 1.0002,
          "ess_bulk": 26643,
          "ess_tail": 33115,
          "failed_checks": []
        },
        {
          "name": "β₁ dose",
          "mean": 0.3154,
          "sd": 0.132,
          "mcse_mean": 0.00081,
          "rhat": 1.0001,
          "ess_bulk": 26504,
          "ess_tail": 33472,
          "failed_checks": []
        }
      ]
    },
    {
      "sampler": "Hamiltonian Monte Carlo",
      "chains": 4,
      "seeds": [
        90125,
        90126,
        90127,
        90128
      ],
      "draws_per_chain": 9000,
      "warmup_discarded_per_chain": 1000,
      "acceptance_rate_by_chain": [
        0.959,
        0.962,
        0.964,
        0.963
      ],
      "parameters": [
        {
          "name": "β₀ (Intercept)",
          "mean": -0.8157,
          "sd": 0.4614,
          "mcse_mean": 0.0015,
          "rhat": 1.0002,
          "ess_bulk": 94387,
          "ess_tail": 29913,
          "failed_checks": []
        },
        {
          "name": "β₁ dose",
          "mean": 0.3149,
          "sd": 0.1326,
          "mcse_mean": 0.00045,
          "rhat": 1,
          "ess_bulk": 87223,
          "ess_tail": 31864,
          "failed_checks": []
        }
      ]
    }
  ],
  "posterior_predictive_checks": {
    "based_on": "36,000 posterior draws of the HMC chains, one replicated trial (10 patients at each of 7 doses) per draw",
    "statistics": [
      {
        "statistic": "Total improved",
        "p_value": 0.503,
        "mcse": 0.0023
      },
      {
        "statistic": "Dose-to-dose drops",
        "p_value": 0.387,
        "mcse": 0.0018
      },
      {
        "statistic": "Improved at dose 0",
        "p_value": 0.509,
        "mcse": 0.0022
      },
      {
        "statistic": "Improved at dose 6",
        "p_value": 0.382,
        "mcse": 0.0021
      },
      {
        "statistic": "Chi-square discrepancy",
        "p_value": 0.741,
        "mcse": 0.0025
      }
    ]
  }
}

2 · Method comparison

Four roads against one exact answer

The reference is the four-chain HMC run (36,000 draws). Because the posterior has only two dimensions it can also be computed exactly on a grid, which shows how much of each distance is Monte Carlo noise in the reference itself.

Posterior summaries of each method and distances to the HMC reference and to the exact posterior
Method	β₀ mean ± MCSE	β₁ mean ± MCSE	β₁ sd	β₁ 95% interval	β₁ shift (sd)	W₁(β₁) to HMC	W₁(β₁) to exact	Gaussian KL to HMC
Laplace (GLM)Gaussian approximation (closed form)	-0.7801	0.3016	0.1292	[0.048, 0.555]	-0.101	0.0133	0.0128	0.00630
Metropolis-Hastings4 chains × 50,000 draws = 200,000	-0.8159 ± 0.0028	0.31536 ± 0.00081	0.1320	[0.065, 0.585]	0.003	0.0014	0.0011	0.00019
HMC (reference)4 chains × 9,000 draws = 36,000	-0.8157 ± 0.0015	0.31491 ± 0.00045	0.1326	[0.062, 0.582]	ref.	ref.	0.0009	ref.
Expectation propagationGaussian approximation (closed form)	-0.8139	0.3142	0.1319	[0.056, 0.573]	-0.005	0.0023	0.0027	0.00011
Exact (quadrature)2-D grid, trapezoidal rule	-0.8142	0.3143	0.1322	[0.063, 0.582]	-0.004	0.0009	–	0.00009
Noise floor: HMC chains 1–2 vs 3–4						0.0011	0.0009	0.00017

How to read it. The shift is the difference from the HMC mean in units of the posterior sd, which makes it an effect size. Laplace sits 0.10 sd low on the slope, and every other method is within 0.006 sd. W₁ is the 1-Wasserstein distance between marginals, in the units of β₁. The noise floor row is the distance between two halves of the reference. By the triangle inequality, a method’s distance to the HMC reference can differ from its distance to the exact posterior by at most the reference’s own distance to it (0.0009).

MH lands at the noise floor, as an exact sampler should. The two Gaussian approximations differ in a way the means alone hide. EP’s KL divergence from the exact posterior is 0.00274, within 0.4% of the 0.00273 of the best Gaussian possible (the one with the exact mean and covariance). Laplace’s is 0.0085, about 3.1 times larger. EP matches the posterior’s moments, while Laplace matches its curvature at the mode, and the posterior is slightly skewed.

How far is each Gaussian from the exact posterior?

KL divergence from the exact posterior to each Gaussian
Gaussian	KL(exact ‖ Gaussian)
Laplace (GLM)	0.00850
Expectation propagation	0.00274
Best possible Gaussian (the exact mean and covariance)	0.00273

Computed by quadrature on a 401 × 1,601 grid, so these carry no Monte Carlo error.

3 · Posterior predictive checks

Does the fitted model reproduce the trial?

For each of the 36,000 HMC draws, a replicated trial with the same design (10 patients at each dose) is simulated with R’s generator (seed 90125) and compared with the observed table.

The thin bar covers 90% of replicated trials and the thick bar 50%. The tick is the median and the open circle the observed count.

Posterior predictive p-values with Monte Carlo error
Test statistic	obs.	p ± MCSE	95% interval	ESS
Total improved	37	0.5028 ± 0.0023	[0.498, 0.507]	43,346
Dose-to-dose drops	2	0.3870 ± 0.0018	[0.383, 0.391]	38,331
Improved at dose 0	3	0.5086 ± 0.0022	[0.504, 0.513]	41,733
Improved at dose 6	8	0.3824 ± 0.0021	[0.378, 0.386]	40,561
Chi-square discrepancy	–	0.7406 ± 0.0025	[0.736, 0.745]	31,346

p is the posterior predictive p-value, Pr(T(y_rep) > T(y)) + ½ Pr(equal) for the counts and Pr(T(y_rep, θ) ≥ T(y, θ)) for the X² discrepancy. Values near 0 or 1 mean the model cannot reproduce that feature of the data. The MCSE uses the effective sample size of each indicator across the four chains. The X² indicator is 0 or 1, so its interval is a Wilson interval with that effective n. The mid-p indicator of a count is 0, ½ or 1, so its MCSE is its own sd over the square root of its ESS and its interval is p ± 1.96 MCSE. MCSEs here run from 0.0018 to 0.0025.

All 14 cells: observed count, replicated 90% interval, p-value

Observed count and replicated 90% interval for each of the 14 cells of the dose by outcome table
Outcome	dose 0	dose 1	dose 2	dose 3	dose 4	dose 5	dose 6
improved	31–6p 0.51	51–7p 0.25	42–7p 0.62	43–8p 0.77	73–9p 0.31	64–9p 0.68	84–10p 0.38
not improved	74–9p 0.49	53–9p 0.75	63–8p 0.38	62–7p 0.23	31–7p 0.69	41–6p 0.32	20–6p 0.62

The p-values are mid-p values like those of the count statistics. Their Monte Carlo errors, computed the same way, run from 0.0019 to 0.0022 (hover a p-value to see its own).

Total improved p = 0.503

Patients improved across all seven doses (observed 37 of 70). The intercept and slope are fitted to this, so a well-specified model should put it near the middle.

Dose-to-dose drops p = 0.387

How often the count falls from one dose to the next (observed: 5 to 4 and 7 to 6, so 2). Checks whether the noise around a smooth logistic curve is as large as binomial noise predicts.

Improved at dose 0 p = 0.509

The lowest-dose cell (observed 3 of 10): does the curve fit its lower end?

Improved at dose 6 p = 0.382

The highest-dose cell (observed 8 of 10): does the curve fit its upper end?

Chi-square discrepancy p = 0.741

Every check sits comfortably inside its replicated distribution: no cell falls outside its 90% interval, the cell p-values run from 0.227 to 0.773, and no test-statistic p-value is below 0.382. The total improved is close to 0.5 by construction (the intercept and slope are fitted to it), which is why it is useless here and decisive in the 2023 version below.

Pearson's X² = Σ (y − n p)² / (n p (1 − p)) at each draw's probabilities, computed for the observed and the replicated table. An omnibus check of overall fit. Points above the dashed diagonal are draws where a replicated trial fits worse than the real one. With p = 0.741, the observed table fits slightly better than a typical replicated one.

4 · Prior sensitivity

How much does the flat prior matter?

The assignment fixed a flat prior. Here the same posterior is recomputed under independent

\mathcal N(0, s^2)

priors on both coefficients, exactly (by quadrature), for scales from very tight to effectively flat.

A weakly informative $\mathcal N(0, 2.5^2)$ moves the slope by 0.06 posterior sd, so the flat-prior conclusions stand. Priors of scale 1 or tighter move it noticeably: 0.29 sd at scale 1, 0.75 sd at 0.5 and 1.26 sd at 0.25. Only the tightest, a prior that says a one-unit dose increase almost certainly changes the log-odds by less than about 0.5, changes the conclusion: it shrinks the slope to 0.148, and its 95% interval then includes zero.

P(β₁ > 0) stays between 0.966 and 0.993 across every prior: the evidence that dose helps is robust, but its size is estimated from only 70 patients and a tight enough prior can halve it.

Posterior of the dose slope under each prior
Prior on β₀ and β₁	β₁ mean	β₁ sd	β₁ 95% interval	P(β₁ > 0)	odds ratio per dose	shift in β₁ (sd)	shift in β₀ (sd)
N(0, 0.25²)	0.148	0.082	[-0.010, 0.310]	0.966	1.16	-1.26	1.40
N(0, 0.5²)	0.216	0.104	[0.016, 0.422]	0.983	1.24	-0.75	0.85
N(0, 1²)	0.276	0.121	[0.044, 0.518]	0.990	1.31	-0.29	0.34
N(0, 2.5²)	0.307	0.130	[0.059, 0.569]	0.993	1.36	-0.06	0.07
N(0, 5²)	0.312	0.132	[0.062, 0.579]	0.993	1.36	-0.01	0.02
N(0, 10²)	0.314	0.132	[0.062, 0.581]	0.993	1.36	0.00	0.00
flat (as assigned)	0.314	0.132	[0.063, 0.582]	0.993	1.37	ref.	ref.

Dose is on its original 0 to 6 scale, so a prior sd on β₁ is per unit of dose. The intercept gets the same prior, which pulls it towards a 50% improvement rate at dose 0.

5 · The 2023 version

Would these checks have caught the 2023 bugs?

The same four-chain diagnostics and predictive check, run on the code exactly as submitted (see the errata). One kind of check catches one bug each.

R-hat catches the Metropolis bug

The 2023 likelihood ignored five of the seven sampled coefficients, so they wandered freely. Four chains disagree about them at once. Coefficients β₂, β₃, β₄, β₅, β₆ all fail, with R-hat up to 2.76 and a bulk ESS of 5 from 200,000 draws. Checking R-hat would have flagged the problem before any plot was drawn.

Only the predictive check catches the HMC and EP bug

HMC was given the 0/1 label as the number improved. Its chains converge perfectly, to the wrong posterior, and every R-hat is below 1.01. But that posterior predicts about 7 improvers in a trial of 70 where 37 improved, so the total-improved check gives p = 0.000 (95% upper bound 0.007).

Metropolis-Hastings, as submitted

Metropolis-Hastings: posterior summaries and convergence diagnostics over 4 chains of 50,000 post-warm-up draws
Coefficient	mean ± MCSE	R-hat	bulk ESS	tail ESS	Status
β₀ (Intercept)	-1.122 ± 0.012	1.0017	2,285	3,493	passes
β₁ dose1	0.3145 ± 0.0025	1.0012	2,722	3,959	passes
β₂ dose2	-28 ± 34	2.6910	5	11	fails R-hat, bulk ESS, tail ESS
β₃ dose3	-11 ± 20	2.3074	5	11	fails R-hat, bulk ESS, tail ESS
β₄ dose4	-17 ± 18	1.7217	6	22	fails R-hat, bulk ESS, tail ESS
β₅ dose5	22 ± 24	2.2903	5	11	fails R-hat, bulk ESS, tail ESS
β₆ dose6	-15 ± 51	2.7571	5	12	fails R-hat, bulk ESS, tail ESS

Hamiltonian Monte Carlo, as submitted

Hamiltonian Monte Carlo: posterior summaries and convergence diagnostics over 4 chains of 9,000 post-warm-up draws
Coefficient	mean ± MCSE	R-hat	bulk ESS	tail ESS	Status
β₀ (Intercept)	-2.683 ± 0.068	1.0036	519	538	passes
β₁ dose1	-0.043 ± 0.076	1.0052	692	694	passes
β₂ dose2	-0.072 ± 0.075	1.0017	709	720	passes
β₃ dose3	-0.053 ± 0.077	1.0031	675	650	passes
β₄ dose4	-0.019 ± 0.076	1.0029	682	618	passes
β₅ dose5	-0.043 ± 0.077	1.0030	650	631	passes
β₆ dose6	-0.030 ± 0.078	1.0032	663	668	passes

Chain 1 seed 90125
Chain 2 seed 90126
Chain 3 seed 90127
Chain 4 seed 90128

Metropolis-Hastingsβ₆ dose6

Warm-up: iterations 1 to 200

After warm-up: iterations 50,001 to 100,000 (1 in every 167 draws plotted)

Rank plot after warm-up (20 bins; dashed line = what perfect mixing gives)

Chain 1

Chain 2

Chain 3

Chain 4

Total improved under the 2023 HMC posterior p = 0.000

Posterior predictive p-values with Monte Carlo error
Test statistic	obs.	p ± MCSE	95% interval	ESS
Total improved	37	0.000	[0.000, 0.007]	519
Dose-to-dose drops	2	0.5192 ± 0.0032	[0.513, 0.525]	13,050
Improved at dose 0	3	0.0941 ± 0.0041	[0.086, 0.102]	3,993
Improved at dose 6	8	0.000403 ± 0.000089	[< 0.001, < 0.001]	36,055
Chi-square discrepancy	–	0.000	[0.000, 0.007]	519

6 · Tuning and robustness

HMC step size and path length, and EP's site choices

Evidence behind DR-002 and DR-003. Each HMC setting runs four chains (same seeds) of 5,000 iterations with 1,000 of warm-up, on the corrected model. Efficiency counts L + 1 gradient evaluations per iteration, the number the original function makes (one at the start of each trajectory, one per leapfrog step).

HMC acceptance, convergence and efficiency for each step size and number of steps
L × ε	trajectory	acceptance (range over chains)	max R-hat	β₁ bulk / tail ESS	ESS per 1k gradients	range over chains
5 × 0.2	1.00	0.1% (0.0%–0.1%)	3.879	stuck	0	–
10 × 0.1	1.00	89.8% (89.3%–90.4%)	1.001	34,802 / 9,421	53.5	45.0–58.6
20 × 0.05original	1.00	96.4% (96.0%–96.7%)	1.000	38,517 / 14,016	41.7	38.8–42.9
40 × 0.025	1.00	99.1% (99.0%–99.2%)	1.000	37,433 / 15,369	23.4	21.5–24.2
5 × 0.05	0.25	98.0% (97.8%–98.3%)	1.005	1,604 / 3,326	16.7	11.8–19.8
10 × 0.05	0.50	96.5% (96.2%–97.0%)	1.001	4,945 / 7,934	28.1	22.7–31.4
40 × 0.05	2.00	97.1% (96.9%–97.4%)	1.000	45,163 / 12,012	18.3	15.4–20.2
10 × 0.2	2.00	0.0% (0.0%–0.0%)	n/a	stuck	0	–
3 × 0.5	1.50	0.0% (0.0%–0.0%)	7.120	stuck	0	–
3 × 0.5M = GLM precision	1.50	97.0% (96.8%–97.2%)	1.000	12,731 / 12,759	198.9	188.5–200.2
2 × 0.8M = GLM precision	1.60	91.8% (91.4%–92.2%)	1.000	15,412 / 14,355	299.1	275.2–304.8

The original L = 20, ε = 1/20 accepts 96.4% of proposals and is reliable, but it is not the most efficient choice. With the same identity mass matrix, 10 steps of 0.1 deliver about 1.2 times as much ESS per gradient evaluation (counting the smaller of bulk and tail ESS), between 1.15 and 1.41 times across 10 independent seed sets. Step sizes of 0.2 and above almost never accept a proposal, and count as zero.

The reason is geometry. With an identity mass matrix the leapfrog integrator is only stable when ε is below about twice the smallest posterior standard deviation along a principal axis, which is 0.14 here, because the intercept and slope are strongly correlated (posterior correlation -0.84). Using the GLM’s precision matrix as the mass matrix removes that limit: 2 steps of 0.8 then give about 7 times the original efficiency (6.4 to 7.7 times across the seed sets).

Does the comparison survive a change of seeds? The four settings DR-002 compares, rerun on 10 independent sets of four seeds (90125–90128, then 90129–90132, up to 90161–90164). Ratios are taken within each seed set, so each comparison is paired by seeds.

Efficiency of four HMC settings over 10 independent seed sets, and its ratio to the original setting within each seed set
L × ε	ESS per 1k gradients, median (range)	× original, median (range)
20 × 0.05original	41.7 (39.4–44.0)	1 (ref.)
10 × 0.1	52.5 (50.4–56.0)	1.23 (1.15–1.41)
3 × 0.5M = GLM precision	198.9 (188.6–212.9)	4.77 (4.41–5.11)
2 × 0.8M = GLM precision	296.1 (274.8–313.0)	6.97 (6.38–7.74)

EP: does the answer depend on how the sites start, or on damping?

EP with different initial sites and damping
EP variant	sweeps	β₁ mean	β₁ sd	largest change vs original
Original: identity sites, no damping	5	0.314201	0.131885	ref.
Near-flat initial sites (0.01 I)	5	0.314201	0.131885	2.3e-10
Strong initial sites (100 I)	5	0.314201	0.131885	5.9e-10
Damped updates (0.5)	23	0.314200	0.131885	2.8e-7

Reproducing these numbers

Everything on this page is computed at build time by the TypeScript in web/src/lib/stats/ with the seeds shown. The chains use seeds 90125, 90126, 90127, 90128 and the predictive replications use seed 90125. scripts/diagnostics-reference.R runs the original R functions with the same seeds and the posterior package. The test suite checks every R-hat, ESS and MCSE against it, and checks the prior-sensitivity quadrature against nested adaptive quadrature in R. Design choices are recorded in the decision records.

Download the numbers (JSON)Methods, model card and limitations