hospital rankings

ranking.Rnw

<<>>=
# n: no of operations
# x: no of deaths
# N: no of hospitals
dat<-list(n=c(47,211,810,148,196,360,119,207,97,
              256,148,215),
          x=c(0,8,46,9,13,24,8,14,8,29,18,31),
          N=12)

cat("model
    {
    for(i in 1:N){
    x[i] ~ dbinom(p[i],n[i])
    logit(p[i]) <- gamma[i] 
    gamma[i] ~ dnorm(mu,tau)
    }
    ## priors:
    mu ~ dnorm(0,1.0E-3)
    #tau ~ dgamma(0.001,0.001)
    sigma ~ dunif(0,10)
    ## half-normal prior:
    ##sigma ~ dnorm(0,0.001)T(0,)
    ##sigma <- 1/sqrt(tau)
    tau <- 1/pow(sigma,2)
    for(i in 1:N){
    theta.pred[i] <- exp(gamma[i])/(1+exp(gamma[i]))
    }
}",file="JAGSmodels/hospitaldeaths.jag")

# Number of steps to "tune" the samplers.
n.adapt <- 1000 ## was 10000
# Number of steps to "burn-in" the samplers.
n.burnin <- 2000 ## was 20000            
# Number of chains to run.
n.chains <- 2                   
# Total number of steps in chains to save.
n.savedsteps <- 4000 ## 80000         
# Number of steps to "thin" (1=keep every step).
n.thin<-1
# Steps per chain.
n.perchain = ceiling( ( n.savedsteps * n.thin) / n.chains) 

## vary inits:
#inits<-list(list(mu=0,alpha=-3,beta=-4,tau=0.001),
#            list(mu=10,alpha=5,beta=4,tau=10))

mod <- jags.model( 
  file="JAGSmodels/hospitaldeaths.jag",
  data = dat,
#  inits= inits,
  n.chains = n.chains,
  n.adapt =n.adapt , quiet=T)

update(mod,n.iter=n.burnin)

track.variables=c("mu","sigma","theta.pred")

res <- coda.samples(mod,
                      var = track.variables,
                      n.iter = n.savedsteps,
                      thin = n.thin)

plot(res)

gelman.diag(res)

summary(res,quantiles=c(0.025,0.5,0.975))

mcmcChain<-as.matrix(res)

hosp<-LETTERS[1:12]
hosp.res<-data.frame(hosp,summary(res)$statistics[3:14,1:2])

## sample from 12 distributions, derive 
## distrn. of ranks
n.sim<-100000

samp<-matrix(c(rep(NA,12*n.sim)),ncol=12)

for(i in 1:n.sim){
  for(j in 1:12){
    samp[i,j] <- rnorm(1,mean=hosp.res[j,2],
                       sd=hosp.res[j,3])
  }
}

samp<-as.data.frame(samp)

colnames(samp)<-LETTERS[1:12]

ranks<-matrix(c(rep(NA,12*n.sim)),ncol=12)

for(i in 1:n.sim){
ranks[i,]<-rank(samp[i,])
}
@

<<fig=T>>=
op <- par(mfrow=c(4,3))
for(i in 1:12){
hist(ranks[,i],freq=F,main=(paste("Rank ",i,sep="")))
}

probs<-rep(NA,12)
for(i in 1:12){
counts<-table(ranks[,i]==i)
##P(rank.I==9)
probs[i]<-counts[2]/sum(counts)
}
@

<<fig=T>>=
## probabilities of each rank:
barplot(probs)
@

<<>>=
## 95\% region for I:
counts<-table(2<ranks[,9] & ranks[,9]>11)
1-counts[2]/sum(counts)
@