The 2023 data-encoding and likelihood bugs, and how the corrected analysis fixes them

• Status: Accepted
• Date recorded: 2026-10-06
• Decision in one line: Keep the 2023 code and its mis-encoded inputs reproducible as the "As submitted (2023)" version, and add a "Corrected" version that calls the very same original functions with seven grouped binomial observations and an intercept-plus-slope design, so the fix is a change of data, not of algorithms.

Context

The assignment asked for a binomial logistic regression of improvement on dose (seven doses, ten patients each, flat prior), fitted four ways. When I re-ran the 2023 R Markdown in 2026, each method turned out to have been given a differently broken version of the data.

• dose was converted to a seven-level factor, so every method fitted seven coefficients (an intercept and six dose offsets) instead of an intercept and a slope.
• The GLM response cbind(count, count - ifelse(improvement == 1, count, 0)) counted every patient as a success and counted the non-improvers again as failures. At dose 0 that reads as 10 successes and 7 failures.
• The Metropolis likelihood used only the first two of the seven sampled coefficients, with the factor level index 1 to 7 as the covariate. The other five coefficients never entered the likelihood and drifted as an unconstrained random walk under the flat prior.
• The Metropolis proposal covariance came from the double-counted GLM, so steps for the slope were about ten times wider than its posterior sd. Only 5.0% of proposals were accepted.
• HMC and EP were given the 0/1 outcome label as the number improved and the cell count as the number of trials, which describes 7 improvers in 70 patients instead of 37.
• The comparison figure labelled HMC as "MHC" and EP as "ABC".

The algorithms themselves (Metropolis.fn, HMC.fn, EP.logit, adapted from the lecture code) were sound.

Decision

The website keeps two versions side by side. "As submitted (2023)" feeds each method exactly the inputs the Rmd built, bugs included, and reproduces the printed 2023 numbers. "Corrected" (the default) passes the original, unmodified functions the grouped data, $y = (3, 5, 4, 4, 7, 6, 8)$ improvers out of $n = 10$ at doses 0 to 6, with the design model.matrix(~ dose) on numeric dose. The submitted Rmd in coursework/ is never edited.

Options considered

1. Fix the Rmd in place. This is the quickest way to a correct answer, but it destroys the record of what was submitted, which matters for academic integrity and for honesty about my own work.
2. Show only the corrected analysis. This is cleaner for a visitor, but it hides the mistakes and loses the most instructive part of the project.
3. Rewrite the four methods in a modern framework (Stan, PyMC). This gives better samplers, but then nothing could be compared with the 2023 output, and the "four roads" would no longer be the roads I actually built.
4. Keep both versions and change only the data (chosen). The bugs become visible as differences in inputs, and every corrected number still comes from the 2023 code.

Why

Option 4 is the only one that keeps the original results faithful while still producing a correct analysis. It also makes the lesson precise. Each bug was a data-encoding bug, and fixing the encoding is enough to make all four methods agree. Because the corrected calls use the same functions, the browser port can be tested against R for both versions to about $10^{-11}$.

What happened

• With the grouped data, the four estimates of the dose slope agree to within 0.016 (the original single chains, seed 90125). Laplace gives 0.302 with 95% interval [0.048, 0.555], Metropolis-Hastings 0.317, HMC 0.315 and EP 0.314.
• The 2026 Bayesian-workflow checks show which tools would have caught which bug, and this is the most useful thing the revival taught me. Run with four chains (seeds 90125 to 90128), the submitted Metropolis sampler fails R-hat for the five drifting coefficients (R-hat between 1.72 and 2.76, bulk ESS below 10 from 200,000 draws). The submitted HMC sampler passes every convergence check (all R-hat below 1.01), because it converges perfectly to the posterior of the wrong data. Only the posterior predictive check exposes it. That posterior predicts about 7 improvers in a 70-patient trial where 37 improved, and the total-improved check gives p = 0 (95% upper bound below 0.008).
• The weak point of the corrected analysis is the data, not the code. Seventy patients leave the slope uncertain (95% interval about 0.06 to 0.58 on the log-odds scale), and the posterior predictive checks cannot rule out models with more structure, because there is so little data to contradict them.

What I'd change

• Build one data object and one design matrix, and pass the same object to every method, so that a method cannot silently see different data.
• Simulate a dataset with known coefficients first and check that all four methods recover them before touching the real data.
• Always run at least four chains and check rank-normalised R-hat and ESS before looking at any estimate.
• Run a posterior predictive check against the raw table. Convergence diagnostics say whether a sampler found its posterior, not whether that posterior describes the data.
• Never call as.numeric() on a factor. It returns level indices, not the values.