Is this how to convert odds ratio intervals to risk ratios

Question

I have a binary logistic regression with just one binary factor predictor. I would like to report the 'risk ratio' (probability ratio) with a confidence interval. I'm using R 4.0.4 for statistical computing.

Here's some data:

> dat <- data.frame(record=c(rep(1,11), rep(0,50-11), rep(1,20), 
             rep(0,50-20)), group=as.factor(c(rep("A", 50), 
              rep("B", 50))))
> tapply(dat$record, dat$group, mean)

   A    B 
0.22 0.40

The risk ratio equals 0.40 / 0.22 = 1.82.

A logistic regression (with logit link) returns an odds ratio (OR):

> mod <- glm(record ~ group, data=mydata, family=binomial)
> antilogit(c(coef(mod)[1], sum(coef(mod)))) # probabilities
(Intercept)             
       0.22        0.40 
> exp(c(coef(mod)[2], confint(mod)[2,])) # OR and confidence interval
Waiting for profiling to be done...
  groupA    2.5 %   97.5 % 
2.363636 0.998533 5.826774 

The antilogit function converts logits to probabilities:

> antilogit <- function(x) {
  exp(x)/(1+exp(x))
} 

I can convert the OR to a risk ratio (RR) using the formula (e.g. Zhang & Yu 1998. What's the relative risk? JAMA 280) RR = OR / (1 - p0 + p0 * OR):

> 2.363636 / (1-0.22+0.22*2.363636)
[1] 1.818182

Converting the OR confidence interval to a RR one is tricky because the intercept (p0) changes with the slope (OR). This is wrong:

> ORINT=exp(confint(mod)[2,]) # OR interval
Waiting for profiling to be done...
> c(ORINT[1]/(1-0.22+0.22*ORINT[1]), ORINT[2]/(1-0.22+0.22*ORINT[2]))
    2.5 %    97.5 % 
0.9988554 2.8259379

I thought to estimate the intercepts for the OR interval using glm():

> dat$group=as.numeric(dat$group)-1 # convert factor to dummy variable
> ORINT=confint(mod)[2,] # log(OR) interval
> # estimate p0 for log(OR) interval
> p0INT=c(coef(glm(record ~ offset(ORINT[1]*group), data=dat, 
                   family=binomial)),
   coef(glm(record ~ offset(ORINT[2]*group), data=dat, 
    family=binomial))) 
> p0INT=antilogit(p0INT) # convert to probabilities
> p0INT
(Intercept) (Intercept) 
  0.3101570   0.1377992 
> ORINT=exp(ORINT) # convert to OR
> c(ORINT[1]/(1-p0INT[1]+p0INT[1]*ORINT[1]), 
    ORINT[2]/(1-p0INT[2]+p0INT[2]*ORINT[2]))
    2.5 %    97.5 % 
0.9989875 3.4992996

What I have calculated is RR = 1.82 with a 95% confidence interval (1.00, 3.50). Notice that these results are quite different to those for OR because p0 is not small (see the equation above: if p0 is small then OR ~ RR).

For comparison, here are some direct risk ratio estimates from a log binomial and a quasipoisson model (e.g. Naimi & Whitcomb 2020. Estimating risk ratios and risk differences using regression. AJE 189):

> mod=glm(record ~ group, data=dat, family=binomial(log))
> exp(c(coef(mod)[1], sum(coef(mod)))) # probabilities
(Intercept)             
       0.22        0.40 
> exp(coef(mod)[2]) # RR
   group 
1.818182 
> exp(confint(mod))[2, ] # RR interval
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.9989799 3.5578943 
> mod=glm(record ~ group, data=dat, family=quasipoisson) # quasipoisson is used to adjust for underdispersion
> exp(c(coef(mod)[1], sum(coef(mod)))) # probabilities
(Intercept)             
       0.22        0.40 
> exp(coef(mod)[2]) # RR
   group 
1.818182 
> exp(confint(mod))[2, ] # RR interval
Waiting for profiling to be done...
    2.5 %    97.5 % 
0.9940394 3.4519283

My conversion from OR interval to RR interval seems OK. Any comments? I haven't thought about how this method could be implemented for more complicated models (and it's not a priority now).

Yes, this is a lot of work for a simple dataset (but it's easy to write a function to do the conversion). The reason I want to convert ORs to RRs is because I sometimes encounter complete separation (p = 0 for one of the groups) and then use penalized-likelihood (Firth 1993. Bias reduction of maximum likelihood estimates. Biometrika 80) to obtain sensible estimates. The R packages logistf 1.24 and brglm 0.7.2 do penalized likelihood but will not accept a log link (and can't estimate RR directly). The R package brglm2 0.7.1 does penalized likelihood and integrates well with the glm() function (it can accept a log link) but doesn't compute profile penalized-likelihood intervals (and Wald tests and intervals for logistic regression are quite unreliable). The R package brglm2 0.7.1 also will not accept quasipoisson. I also had a look at the probratio() function in the R package epitools 0.5.10.1, which converts OR to RR but returns Wald intervals.

I could have simply used a chi-squared test or prop.test() to estimate differences:

> tab=table(dat)
> tab
      group
record  0  1
     0 39 30
     1 11 20
> chisq.test(tab, correct=FALSE)

        Pearson's Chi-squared test

data:  tab
X-squared = 3.7868, df = 1, p-value = 0.05166

> prop.test(c(11,20), c(50,50), correct=FALSE)

        2-sample test for equality of proportions without continuity correction

data:  c(11, 20) out of c(50, 50)
X-squared = 3.7868, df = 1, p-value = 0.05166
alternative hypothesis: two.sided
95 percent confidence interval:
 -0.357828257 -0.002171743
sample estimates:
prop 1 prop 2 
  0.22   0.40 

However, RR has been reported in previous work and is the most informative measure of the contrast between groups for my study.

Answer 1

The above proposal is statistical overkill.

I finally discovered that unconditional score statistic methods can be applied to the difference in two independent binomial proportions AND relative risk AND the odds ratio. Unconditional score statistic methods have been demonstrated to perform well for all of these estimates, e.g. Agresti & Min 2002, Biostatistics 3, Newcombe 1998 Stat. Med. 17.

Unconditional score statistic methods (and others) can be computed using the PropCIs R package. I found that my logistic regression riskratio intervals were more conservative (wider than) those from the riskscoreci function in PropCIs.

Answer 2

The OP started with odds ratios and I assume a logit link best fits the data (log links are unlikely to fit as well in general because they require constraints to keep probabilities legal). Once you recognize that risk ratios are covariate-dependent (covariate-specific) it's easy. Just pick settings for all the covariates except for the one you are changing, let the one change, computed estimated risks directly from the logistic model, and take ratios. Getting uncertainties is the only hard part. For Bayesian models, posterior intervals for covariate-specific risk ratios may be computed exactly from the posterior draws of the underlying log-odds-scale parameters. In the frequentist domain you could use the bootstrap or the $\delta$ method.

The OP started with the no-covariate case but it's best to learn more general methods IMHO. For one thing there is likely to be risk heterogeneity due to baseline variable distributions, and simple ratios were intended only for the homogeneous case.