Bayesian Hierarchical Modeling via STAN

Deaths after heart attack

1 The dataset

The following pieces of explanatory text were taken from Jim Albert’s github page at https://bayesball.github.io/BOOK/bayesian-hierarchical-modeling.html. Particularly, Section 10.3 on Hierarchical Beta-Binomial Modeling, and Section 10.3.1 Example: Deaths after heart attack.

The New York State (NYS) Department of Health collects and releases data on mortality after Acute Myocardial Infarction (AMI), commonly known as a heart attack. Their 2015 report was the initial public data release by the NYS Department of Health on risk- adjusted mortality outcomes for AMI patients at hospitals across the state. We focus on 13 hospitals in Manhattan, New York City - Manhattan in 2015, with the goal of learning about the percentages of resulted deaths from heart attack for hospitals in this sample. Below we have the records for each hospital the number of heart attack cases, the corresponding number of resulted deaths, and their computed percentage of resulted deaths. NYP stands for New York Presbyterian.

• Bellevue Hospital Center
• Harlem Hospital Center
• Lenox Hill Hospital
• Metropolitan Hospital Center
• Mount Sinai Beth Israel
• Mount Sinai Hospital
• Mount Sinai Roosevelt
• Mount Sinai St. Luke’s
• NYU Hospitals Center
• NYP Hospital - Allen Hospital
• NYP Hospital - Columbia Presbyterian Center
• NYP Hospital - New York Weill Cornell Center
• NYP/Lower Manhattan Hospital

y = c(4,1,18,7,24,16,6,19,15,13,25,11,4)
n = c(129,35,228,84,291,270,46,293,241,105,353,250,41)
K = length(y)
phat = y/n


par(mfrow=c(1,1))
plot(y/n,ylim=c(0,0.13),xlab="Number of cases per hospital",ylab="Proportion of death",pch=16,cex=1.5,axes=FALSE)
axis(2);axis(1,at=1:13,label=n);box()
text(1:13,rep(0,13),y)
abline(h=sum(y)/sum(n),col=6)

Treating “cases” as trials and “deaths” as successes, the Binomial sampling model is a natural choice for this data, and the objective is to learn about the death probability \(p\) of the hospitals. If one looks at the actual death percentages, some hospitals have much higher death rates than other hospitals. For example, the highest death rate belongs to Mount Sinai Roosevelt, at 13.043% which is more than four times the rate of Harlem Hospital Center at 2.857%.

1. Pooling the data from the 13 hospitals: If one assumes a common probability p for all 13 hospitals, this model does not allow for possible differences between the death rates among these hospitals.

2. Analysing each hospital separately: If one creates 13 separate Binomial sampling models, one for each hospital, and conducts separate inferences, one loses the ability to use potential information about the death rate from hospital \(j\) when making inference about that of a different hospital \(i\). Since these are all hospitals in Manhattan, New York City, they may share attributes in common related to death rates from heart attack. The separate modeling approach does not allow for the sharing of information across hospitals.

3. Hierarchical modeling A hierarchical model provides a compromise between the combined and separate modeling approaches. In this setting, one builds a hierarchical model by assuming the hospital death rate parameters a priori come from a common distribution. Specifically, one builds a hierarchical model based on a common Beta distribution that generalizes the Beta-Binomial conjugate model. This modeling setup provides posterior estimates that partially pool information among hospitals.

2 The STAN file

Save the following chunck of code into a stan file named, for instance, “hierarchical-beta-binomial.stan”.

data {
  int<lower=1> K;               // number of observations
  int<lower=0> y[K];            // successes
  int<lower=0> n[K];            // trials
  real<lower=0> n1;             // fixed constant for eta prior
}

parameters {
  real<lower=0,upper=1> mu;     // mean parameter of Beta
  real<lower=0> eta;            // precision parameter (a+b)
  vector<lower=0,upper=1>[K] p; // individual probabilities
}

model {
  // ---- Prior for mu ----
  // mu ~ Uniform(0,1) is implicit since Stan parameters are uniform
  // on bounded intervals unless otherwise specified

  // ---- Prior for eta ----
  // p(eta) ∝ n1 / (n1 + eta)^2
  target += log(n1) - 2 * log(n1 + eta);

  // ---- Hierarchical model ----
  // p[i] ~ Beta(mu * eta, (1 - mu) * eta)
  p ~ beta(mu * eta, (1 - mu) * eta);

  // y[i] ~ Binomial(n[i], p[i])
  y ~ binomial(n, p);
}

3 MCMC via STAN

#install.packages("rstan")
library("rstan")
options(mc.cores=4)

stan_data = list(
  K = K,
  n = n,
  y = y,
  n1 = 10   
)

# Running stan code
model = stan_model("hierarchical-beta-binomial.stan")

fit = sampling(model,data=stan_data,seed=31415,iter=15000,warmup=10000,thin=1,chains=4)

params = extract(fit)

mu  = params$mu
eta = params$eta
a = mu*eta
b = (1-mu)*eta
sd = sqrt(a*b/((a+b)^2*(a+b+1)))

4 Posterior summary

print(fit)

## Inference for Stan model: anon_model.
## 4 chains, each with iter=15000; warmup=10000; thin=1; 
## post-warmup draws per chain=5000, total post-warmup draws=20000.
## 
##          mean se_mean      sd    2.5%     25%     50%     75%   97.5% n_eff
## mu       0.07    0.00    0.01    0.06    0.07    0.07    0.08    0.09  7980
## eta    555.72   36.42 1008.17   37.56  114.55  225.74  516.44 3943.77   766
## p[1]     0.06    0.00    0.01    0.03    0.05    0.06    0.07    0.08  6140
## p[2]     0.06    0.00    0.02    0.03    0.05    0.07    0.07    0.10 15146
## p[3]     0.07    0.00    0.01    0.05    0.07    0.07    0.08    0.10 14375
## p[4]     0.07    0.00    0.02    0.05    0.06    0.07    0.08    0.11 14691
## p[5]     0.08    0.00    0.01    0.06    0.07    0.08    0.08    0.10 10353
## p[6]     0.06    0.00    0.01    0.04    0.06    0.06    0.07    0.09 12927
## p[7]     0.08    0.00    0.02    0.05    0.07    0.08    0.09    0.14  7696
## p[8]     0.07    0.00    0.01    0.05    0.06    0.07    0.07    0.09 17380
## p[9]     0.07    0.00    0.01    0.04    0.06    0.07    0.07    0.09 14623
## p[10]    0.09    0.00    0.02    0.06    0.07    0.09    0.10    0.14  5308
## p[11]    0.07    0.00    0.01    0.05    0.06    0.07    0.08    0.09 13859
## p[12]    0.06    0.00    0.01    0.03    0.05    0.06    0.07    0.08  6272
## p[13]    0.08    0.00    0.02    0.04    0.06    0.07    0.09    0.13 13754
## lp__  -597.18    0.15    5.49 -607.04 -600.86 -597.60 -594.00 -584.58  1260
##       Rhat
## mu       1
## eta      1
## p[1]     1
## p[2]     1
## p[3]     1
## p[4]     1
## p[5]     1
## p[6]     1
## p[7]     1
## p[8]     1
## p[9]     1
## p[10]    1
## p[11]    1
## p[12]    1
## p[13]    1
## lp__     1
## 
## Samples were drawn using NUTS(diag_e) at Thu Oct 30 17:58:17 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

5 MCMC chains and marginal posteriors

par(mfrow=c(2,3))
ts.plot(mu)
ts.plot(eta)
ts.plot(sd)
hist(mu,prob=TRUE,main=expression(mu))
hist(eta,prob=TRUE,main=expression(eta))
hist(sd,prob=TRUE,main="Standard deviation")

6 Marginal posterior summaries

tab = apply(cbind(mu,eta,sd),2,quantile,c(0.05,0.5,0.95))
colnames(tab) = c("mu","eta","sd")
rownames(tab) = c("5%","50%","95%")
tab

##             mu        eta          sd
## 5%  0.05951586   49.21486 0.005366085
## 50% 0.07151268  225.73636 0.017042285
## 95% 0.08742623 2234.30645 0.038096447

7 Posterior for death proportions

p = params$p

pmean = apply(p,2,mean)

par(mfrow=c(1,1))
boxplot(p,outline=FALSE,ylim=c(0,0.15))
points(1:13,phat,col=2,pch=16)
abline(h=mean(phat),col=2,lwd=2)

8 Borrowing strength and shrinkage effect

plot(rep(0.1,K),phat,xlim=c(0,0.3),axes=FALSE,xlab="",ylab="",pch=16,ylim=c(0,0.14))
points(rep(0.2,K),pmean,col=2,pch=16)
for (i in 1:K)
  segments(0.1,phat[i],0.2,pmean[i],col=gray(0.8),lwd=2)
points(rep(0.1,K),phat,pch=16)
points(rep(0.2,K),pmean,col=2,pch=16)
text(0.1,0.01,"y/n")
text(0.2,0.01,"E(p)")
box()