Skip to content

Instantly share code, notes, and snippets.

@vasishth vasishth/ranking.Rnw

Created Nov 25, 2014
Embed
What would you like to do?
hospital rankings
<<>>=
# 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)
@
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.