# 2012.11.06    from WinBUGS Examples Vol I in help menu

Surgical: Institutional ranking

#This example considers mortality rates in 12 hospitals performing 
# cardiac surgery in babies. The data are shown below.


#The number of deaths ri for hospital i are modelled as a binary 
# response variable with `true' failure probability pi:


# We first assume that the true failure probabilities are 
# independent  (i.e.fixed effects) for each hospital. This is 
# equivalent to assuming a standard non-informative prior distribution 
# for the pi's, namely:

#		pi  ~  Beta(1.0, 1.0)

model
		{
		   for( i in 1 : N ) {
		      p[i] ~ dbeta(1.0, 1.0)
			r[i] ~ dbin(p[i], n[i])
		   }
		}
		
		
#data
list(n = c(47, 148, 119, 810,	211, 196, 148, 215, 207, 97, 256, 360),
	     r  = c(0, 18, 8, 46, 8, 13, 9,	31, 14, 8, 29, 24),
	     N = 12)

	list(p = c(0.1, 0.1, 0.1, 0.1,	0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1))


######

# A more realistic model for the surgical data is to assume that the failure 
# rates across hospitals are similar in some way. This is equivalent to specifying 
# a random effects model for the true failure probabilities pi as follows:

#		logit(pi)  =  bi
		
#		bi  ~  Normal(m, t)
		
# A standard non-informative prior is then specified for the population 
# mean (logit) probability of failure, m, with two alternative "noninformative" 
# priors considered for the random effects variance: prior 1 is a uniform prior 
# on the standard deviation, and prior 2 is a gamma(0.001, 0.001) prior on the precision, t.


model
		{
			for( i in 1 : N ) {
				b[i] ~ dnorm(mu,tau)
				r[i] ~ dbin(p[i],n[i])
				logit(p[i]) <- b[i]
			  	}
			  pop.mean <- exp(mu) / (1 + exp(mu))
			  mu ~ dnorm(0.0,1.0E-6)
			  # Choice of prior on random effects variance
			  # Prior 1: uniform on SD
			  sigma~ dunif(0,100)
			  tau<-1/(sigma*sigma)
			
			  #Prior 2: 
			  #tau ~ dgamma(1.0E-3, 1.0E-3);
			  #sigma <- 1/sqrt(tau); # s.d. of random effects
		} 


list(b = c(	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1), 
	     sigma = 1, mu = 0)
	     
	     
list(b = c(	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1,	0.1), 
	     tau = 1, mu = 0)




# A particular strength of the Markov chain Monte Carlo 
# (Gibbs sampling) approach implemented in BUGS is the ability 
# to make inferences on arbitrary functions of unknown model
# parameters. For example, we may compute the rank probabilty 
# of failure for each hospital at each iteration. This yields 
# a sample from the posterior distribution of the ranks. 

# The figures below show the posterior ranks for the estimated 
# surgical mortality rate in each hospital for the random effect models.
# These are obtained by setting the rank monitor for variable p 
# (select the "Rank" option from the "Statistics" menu) after the 
# burn-in phase, and then selecting the "histogram" option from 
# this menu after a further 10000 updates. These distributions 
# illustrate the considerable uncertainty associated with 
# 'league tables': there are only 2 hospitals (H and K) whose 
# intervals exclude the median rank and none whose intervals 
# fall completely within the lower or upper quartiles. 

# Plots of distribution of ranks of true failure probability 
# for random effects model:



#############################################################  R JAGS


HospitalModel=readLines(,18)
model
		{
			for( i in 1 : N ) {
				b[i] ~ dnorm(mu,tau)
				r[i] ~ dbin(p[i],n[i])
				logit(p[i]) <- b[i]
			  	}
			  pop.mean <- exp(mu) / (1 + exp(mu))
			  mu ~ dnorm(0.0,1.0E-6)
			  # Choice of prior on random effects variance
			  # Prior 1: uniform on SD
			  sigma~ dunif(0,100)
			  tau<-1/(sigma*sigma)
			
			  #Prior 2: 
			  #tau ~ dgamma(1.0E-3, 1.0E-3);
			  #sigma <- 1/sqrt(tau); # s.d. of random effects
		}
		
#		


writeLines(HospitalModel,'HospitalModel.txt')

n = c(47, 148, 119, 810,211,196,148,215,207,97,256,360)
r = c( 0,  18,   8,  46,  8, 13,  9, 31, 14, 8, 29, 24)
N = 12
	     
INITS=list(b = c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1), 
	     sigma = 1, mu = 0)
 
library('rjags')

jags <- jags.model("HospitalModel.txt",
             data = list(n=n, r=r, N=N),
             inits=INITS,
             n.chains = 1,
             n.adapt = 4000)
 
update(jags, 10000)
 
results=jags.samples(jags, c('p','mu','sigma'), 10000)
str(results)

summary(results$mu[1,,1])
summary(results$sigma[1,,1])

apply(results$p,c(1,3),median)
apply(results$p,c(1,3),sd)

ranks=apply(results$p,c(2,3),rank)

ranks[,1:5,1]

# graph the rankings

par(mfrow = c(3,4))   # multiple panels

for(hosp in 1:12) {
	fr = table(ranks[hosp,,1])
	r = as.numeric(dimnames(fr)[[1]])
	plot(r, table(ranks[hosp,,1]), xlim=c(1,12) ,
	     type="h", xlab="Rank", ylab="prob", main=paste("Hospital",toString(hosp)) )
}