Surgical: Institutional Ranking

An example from OpenBUGS [38] concerning mortality rates in 12 hospitals performing cardiac surgery in infants.

Model

Number of deaths are modelled as

r_i &\sim \text{Binomial}(n_i, p_i) \quad\quad i=1,\ldots,12 \\
\operatorname{logit}(p_i) &= b_i \\
b_i &\sim \text{Normal}(\mu, \sigma) \\
\mu &\sim \text{Normal}(0, 1000) \\
\sigma^2 &\sim \text{InverseGamma}(0.001, 0001),

where r_i is the number of deaths, out of n_i operations, at hospital i.

Analysis Program

using Mamba

## Data
surgical = (Symbol => Any)[
  :r => [0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24],
  :n => [47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360]
]
surgical[:N] = length(surgical[:r])


## Model Specification

model = Model(

  r = Stochastic(1,
    @modelexpr(n, p, N,
      Distribution[Binomial(n[i], p[i]) for i in 1:N]
    ),
    false
  ),

  p = Logical(1,
    @modelexpr(b,
      invlogit(b)
    )
  ),

  b = Stochastic(1,
    @modelexpr(mu, s2,
      Normal(mu, sqrt(s2))
    ),
    false
  ),

  mu = Stochastic(
    :(Normal(0, 1000))
  ),

  pop_mean = Logical(
    @modelexpr(mu,
      invlogit(mu)
    )
  ),

  s2 = Stochastic(
    :(InverseGamma(0.001, 0.001))
  )

)

## Initial Values
inits = [
  [:r => surgical[:r], :b => fill(0.1, surgical[:N]), :s2 => 1, :mu => 0],
  [:r => surgical[:r], :b => fill(0.5, surgical[:N]), :s2 => 10, :mu => 1]
]


## Sampling Scheme
scheme = [NUTS([:b]),
          Slice([:mu, :s2], [1.0, 1.0])]
setsamplers!(model, scheme)


## MCMC Simulations
sim = mcmc(model, surgical, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

Results

Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750

Empirical Posterior Estimates:
             Mean         SD      Naive SE       MCSE        ESS
      mu -2.550741500 0.1567451 0.0018099366 0.0032183095 2372.0965
pop_mean  0.073072811 0.0104407 0.0001205589 0.0002104121 2462.1738
    p[1]  0.052805725 0.0198027 0.0002286619 0.0004412988 2013.6467
    p[2]  0.104001890 0.0219549 0.0002535130 0.0003555323 3813.3232
    p[3]  0.070669708 0.0176251 0.0002035168 0.0002098008 7057.4433
    p[4]  0.059118030 0.0079059 0.0000912898 0.0001290882 3750.8750
    p[5]  0.051504980 0.0132372 0.0001528496 0.0002512701 2775.2895
    p[6]  0.069199820 0.0147703 0.0001705526 0.0002018471 5354.6699
    p[7]  0.067117588 0.0160309 0.0001851085 0.0001921212 6962.4719
    p[8]  0.123097365 0.0219326 0.0002532561 0.0003878670 3197.5350
    p[9]  0.070410452 0.0148534 0.0001715122 0.0001833296 6564.2635
   p[10]  0.078694851 0.0198460 0.0002291623 0.0002492771 6338.4458
   p[11]  0.102801934 0.0176951 0.0002043250 0.0002605003 4614.1109
   p[12]  0.068715136 0.0119246 0.0001376936 0.0001584527 5663.5606
      s2  0.194478287 0.1732356 0.0020003522 0.0056335715  945.5980

Quantiles:
            2.5%       25.0%        50.0%       75.0%      97.5%
      mu -2.886808 -2.641133981 -2.542990063 -2.4516929 -2.26518316
pop_mean  0.052810  0.066537569  0.072898833  0.0793148  0.09404783
    p[1]  0.016802  0.038901688  0.052307730  0.0654430  0.09367208
    p[2]  0.067526  0.087996977  0.101695407  0.1173628  0.15350268
    p[3]  0.039639  0.058578620  0.069377099  0.0811090  0.10938201
    p[4]  0.044485  0.053534322  0.058953960  0.0642819  0.07502250
    p[5]  0.027655  0.042219334  0.050835609  0.0601217  0.07939497
    p[6]  0.043011  0.058938078  0.068255749  0.0784358  0.10189272
    p[7]  0.038490  0.056074218  0.066107926  0.0769301  0.10139561
    p[8]  0.084454  0.107222812  0.121811458  0.1372187  0.16871609
    p[9]  0.044307  0.059573839  0.069798505  0.0799364  0.10182827
   p[10]  0.044685  0.064689236  0.076790598  0.0905553  0.12280446
   p[11]  0.072230  0.090317164  0.101265471  0.1140400  0.14092283
   p[12]  0.047187  0.060534280  0.068033520  0.0762521  0.09388019
      s2  0.030929  0.090703538  0.147636980  0.2405679  0.62218529