How to implement credible 95% interval for median odds ratio using JAGS? As described in Merlo et al (J Epidem Comm Health 2006), the 95% credible interval for MOR is calculated using MCMC. MOR is defined as $\exp(\sqrt{2\sigma^2}\times 0.675)$, where $\sigma$ is the level-2 variance of the random intercept $u$ from a null model of a hierarchical logistic regression.  
Does anyone have an idea of how to write a program for an Markov chain Monte Carlo to calculate the standard error of the  median odds ratio (MOR) using rjags?
My dependent variable is outcome(alive/dead) and the clustering (level2)variable is Hospital. There are 140 hospitals and would like to see variations in outcome between hospitals. Other risk factors will be included later as independent level1 variables.
 A: I don't know if this is a solution for you, but since the lme4 glmer function can provide random intercept posterior median estimates and their conditional variance - and under the assumption of normality (for random effect), posterior median = posterior mode - wouldn't it be valid to do a parametric bootstrap repeatedly drawing from the estimated posterior distribution to obtain a confidence interval for the MOR?
Merlo describes MOR conceptually as the median of all possible pairwise odds ratios (with the larger of the pair in the numerator in all cases) - so 
for random intercept posterior estimates on the log scale - 
this would be median of:
exp(abs(X(a) - X(b))) for all a,b pairs of clusters.
-- to test this interpretation, in a fitted model I compared this 
estimate to what was produced by the approximation based on area variance estimate alone in Merlo (eq 6).  I did get very similar results, 1.623 from the median of 1000 bootstrap replicates vs. 1.629 based on the Merlo approximation.
t<- ranef(fit.1,condVar = TRUE)

est<-as.numeric(unlist(t$clientid.x))

var<- as.numeric(unlist(attr(t$clientid.x,"postVar")))

### bootstrap 
### create empty output collection:
mor_boot<- c()
#### iterate over replicates
for(i in 1:1000){
### draw vector of area random effects from normal
drw<- rnorm(n = length(est),mean = est,sd = sqrt(var))
### create data frame with all possible pairs 
s<- combn(drw,2)
#### estimate MOR and save in output
mor_boot<- c(mor_boot, median(exp(abs(s[1,]- s[2,]))))
}
### bootstrap median and 95% CI
quantile(mor_boot, c(.025,.5,.975))

