How to calculate SE of an odds ratio

Question

If one is calculating odds ratio with a,b,c and d counts, I believe variance of log(OR) is given by

var_log_OR = (1/a + 1/b + 1/c + 1/d)

Hence one can calculate 95% confidence intervals of OR as follows:

SE_log_OR = sqrt(var_log_OR)

CI_lower_log_OR = log(OR) - 1.96*SE_log_OR
CI_upper_log_OR = log(OR) + 1.96*SE_log_OR

CI_lower_OR = exp(CI_lower_log_OR)
CI_upper_OR = exp(CI_upper_log_OR)

But how can we calculate SE of OR?

Answer 1

@FrankHarrell is right that the standard error for an odds ratio is a problematic number in the sense that you can do better by testing on the corresponding log(odds ratio) scale, as the sampling distribution of the log(odds ratio) is more likely to be normally distributed.

Nonetheless, the standard error of the odds ratio does exist, even if it is not that useful. One possible estimate is to use the delta method to move from the standard error of the log(odds ratio) to an approximation of the standard error of the odds ratio.

$\sqrt{(1/a + 1/b + 1/c + 1/d)}\times\frac{a\times d}{b\times c}$

Answer 2

OR is not a valid quantity to compute a SE of in the sense that it cannot have a symmetric distribution. Applying +/- SE to it may lead to negative ORs.

Answer 3

Instead of a standard error why not compute the standard deviation of the posterior distribution of the OR? You can solve for it numerically very easily using an MCMC sampler.

Here is some R and JAGS code to do so.

################################################################
###                                                          ###
###      Contingency Table Analysis for Obestity Data        ###
###                                                          ###
################################################################

# Required Pacakges
library("ggplot2")
library("runjags")
library("parallel") # sets parallelization for MCMC


# set up the model
mod = 'model {

################################################
###    Greater than 80th BMI Percentile      ###
################################################

  # marginal likelihood functions
  n11_G ~ dbin(pi_one_G, n1_plus_G)
  n21_G ~ dbin(pi_two_G, n2_plus_G)

  #priors
  pi_one_G ~ dbeta(1,1)
  pi_two_G ~ dbeta(1,1)

  # transformations
  rho_G <- pi_two_G/pi_one_G
  theta_G <- pi_two_G*(1-pi_one_G)/(pi_one_G*(1-pi_two_G))
  delta_G <- pi_two_G-pi_one_G



}'

# set up the data
Dat = list(n1_plus_G = 108,
           n2_plus_G = 88,
           n11_G = 68,
           n21_G = 44)

# Monitor these variables
Vars = c("pi_one_G","pi_two_G","rho_G","theta_G","delta_G")


# set up MCMC parameters
inits1=list(.RNG.name= "base::Wichmann-Hill",
            .RNG.seed= 12341)
inits2=list(.RNG.name= "base::Marsaglia-Multicarry",
            .RNG.seed= 12342)
inits3=list(.RNG.name= "base::Super-Duper",
            .RNG.seed= 12343)
inits4=list(.RNG.name= "base::Mersenne-Twister",
            .RNG.seed= 12344)

chains = 4
burn = 5000
samp = 10000
adapt = 5000
thin = 1

# parallel chains
cl = makeCluster(4)

# MCMC estimation
HjagsOut = run.jags(model = mod, monitor = Vars, data=Dat, n.chains=chains, thin = thin,
                    burnin = burn, sample = samp, adapt=adapt, method="rjparallel",method.options=list(cl=cl),
                    inits=list(inits1,inits2,inits3,inits4))

#summarize results
summary(HjagsOut)

plot(HjagsOut, layout=c(4,2))

The odds ratio parameter ($\theta_\text{G}$) is simply a function of the samples from the binomial parameters. This of course assumes certain study design. In this case I was looking at the difference in children's BMI percentile group (80th and above or below 80th) from a control and experimental group, pre and post intervention treatment. Therefore, the row totals (number of children in the experimental group and control group respectively) were fixed.

The beta(1,1) prior is equivalent to a Uniform(0,1) prior and could easily be changed to the Jeffreys's beta(0.5,0.5) prior or anything you desire.