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
where is the number of deaths, out of operations, at hospital .
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