4
$\begingroup$

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.

Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ A lot more information is necessary such as: the regression model (e.g. how many variables and observations you have), the prior distributions, a link to the paper you mention ... Otherwise the answer is the Reference Manual. $\endgroup$
    – user10525
    Commented Jul 13, 2012 at 16:08
  • $\begingroup$ Start with johnmyleswhite.com/notebook/2010/08/20/… and let us know about any problem you face. 1. Also see rpubs.com/jeromyanglim/rjags_normal_distribution $\endgroup$
    – Stat-R
    Commented Jul 13, 2012 at 16:24
  • $\begingroup$ Hi Omar, I just approved your anonymous edit - note that you have full edit privileges over your own questions so, if you log in, you do not need to wait for approval when editing your question. $\endgroup$
    – Macro
    Commented Jul 14, 2012 at 13:07
  • $\begingroup$ What do you want to do with the standard error ? Can't you be satisfied with a confidence/credibility interval ? $\endgroup$
    – Stéphane Laurent
    Commented Jul 14, 2012 at 13:26
  • $\begingroup$ yes it's for the credibility interval,I thought mcmc provide standard error only, either would do for me . It just that I have no experience in programming MCMC. the paper mentioned that mLwin does it but I don't have the software, that's why I thought of rJags. This is for interpreting result with odds ratio at the level 2 rather using ICC which is not reliable for logistic regression as it is with linear multilevel. Thanks. $\endgroup$
    – Omar
    Commented Jul 15, 2012 at 11:15

1 Answer 1

Reset to default
4
$\begingroup$

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))
Improve this answer
$\endgroup$
1
  • $\begingroup$ This was super helpful conceptually and the code too. I feel a bit troubled by what I found when I did it: I have interaction terms in a random effects model that have lots of unobserved combinations (upwards of 10k/ 90%) and few observations for many of the levels that were observed. It seems the true effect for these interactions should be highly uncertain. But the "CI" for these terms are very small! Is there a risk of pseudoreplication here? $\endgroup$
    – Michael Roswell
    Commented Dec 12, 2019 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.