What went wrong in 2023, and what changed

data <- data.frame(dose = as.factor(doses),
                   improvement = improvement_status,
                   count = counts)
X <- model.matrix(~ dose, data = data)   # 14 x 7

grouped <- data.frame(dose = 0:6,
                      y = c(3, 5, 4, 4, 7, 6, 8),
                      n = rep(10, 7))
X_c <- model.matrix(~ dose, data = grouped)  # 7 x 2

The model in the question is logit p = x′β with dose as the covariate. As a factor, every method instead estimates an intercept plus six separate dose offsets, a saturated model with no notion of “a bigger dose helps more”. The data were also kept in long format, one row per dose and outcome, which set up the encoding bugs below.

glm(cbind(count, count - ifelse(improvement == 1, count, 0)) ~ dose,
    family = binomial, data = data)

glm(cbind(y, n - y) ~ dose, family = binomial, data = grouped)

The first column (successes) is the full count on every row, and non-improvers are added again as failures. At dose 0 that reads as 10 successes and 7 failures, so the fitted probability of improvement is 10/17 = 0.59 instead of the observed 0.30. This GLM also supplied the proposal covariance for Metropolis-Hastings.

log_odds <- as.numeric(data$dose) * params[2] + params[1]
# ... while the sampler proposes all p = ncol(X) = 7 parameters

loglik <- function(params, data) {
  p <- plogis(params[1] + params[2] * data$dose)   # dose = 0..6
  sum(dbinom(data$y, size = data$n, prob = p, log = TRUE))
}

Coefficients 3 to 7 never touch the likelihood, so under the flat prior they drift as an unconstrained random walk (posterior sd up to 54, effective sample size as low as 1.2 from 50,000 draws). They cannot change the acceptance ratio, though. The low acceptance rate of 5.0% comes from the proposal covariance, which was borrowed from the double-counted GLM above: its dose-1 standard error (0.74) gives a proposal sd of 1.25 for the second parameter, about 10 times the posterior sd of the slope it actually plays in the likelihood (0.13). Most jumps overshoot, which is why the six rows printed by head(results_mh) are identical. as.numeric() on a factor also gives the level index 1 to 7, not the dose 0 to 6, shifting the intercept. Interestingly, the success/trial encoding in this one function was right.

HMC.fn(y = data$improvement, n = data$count, X = X, ...)
EP.logit(response = data$improvement, n = data$count, X = X, ...)

HMC.fn(y = grouped$y, n = grouped$n, X = X_c, L = 20, ...)
EP.logit(response = grouped$y, n = grouped$n, X = X_c, ...)

Each long-format row says “1 success out of count trials” or “0 out of count”: seven improvements in seventy patients instead of 37. Both methods therefore agree with each other on a posterior centred near an intercept of -2.65 with dose effects near zero, which is a correct answer to the wrong question.

legend = c('Metropolis Hasting', 'GLM', 'MHC', 'ABC')

legend = c('Metropolis-Hastings', 'Laplace (GLM)', 'HMC', 'EP')

A cosmetic slip, but a telling one: “ABC” usually means approximate Bayesian computation, a different method. The figure also took its axis limits from the Metropolis density, which in this version was very wide for the drifting coefficients.