# Is this how to convert odds ratio intervals to risk ratios

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),sum(coef(mod)))) # probabilities
(Intercept)
0.22        0.40
> exp(c(coef(mod), 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.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-0.22+0.22*ORINT), ORINT/(1-0.22+0.22*ORINT))
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*group), data=dat, family=binomial)),
+   coef(glm(record ~ offset(ORINT*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-p0INT+p0INT*ORINT), ORINT/(1-p0INT+p0INT*ORINT))
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), sum(coef(mod)))) # probabilities
(Intercept)
0.22        0.40
> exp(coef(mod)) # 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), sum(coef(mod)))) # probabilities
(Intercept)
0.22        0.40
> exp(coef(mod)) # 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=F)

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=F)

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.

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.