0
$\begingroup$

I have performed a multiple logistic regression because I wanted to see the association between Death and Cardiovascular disease. I adjusted using age, sex, risk factors. The result came in ODDS RATIO with CONFIDENCE INTERVALS.

How do I do this using RISK RATIO with CONFIDENCE INTERVAL instead?

If I just do the ratio between exposed and non-exposed then I don't adjusted for age, sex, risk factors anymore and also I don't have confidence interval.

Please help with codes and method I use R

$\endgroup$
7
  • $\begingroup$ What is your study design? Cohort? Cross sectional? Case control? These all have implications on what measures of association you should calculate. From your question it is not clear the method from which your data arose. I would caution translating these measures unless it’s truly appropriate $\endgroup$
    – jpsmith
    Feb 24, 2023 at 12:37
  • $\begingroup$ @jpsmith it is a retrospective case control study I would say. It is based on recruited patients for a certain disease and I wanted to see if the presence of cardiovascular disease in these patients is associated to the outcome (death) $\endgroup$
    – user19745561
    Feb 24, 2023 at 12:42
  • $\begingroup$ In a case control you can’t calculate a risk ratio $\endgroup$
    – jpsmith
    Feb 24, 2023 at 12:43
  • $\begingroup$ @jpsmith maybe it is cohort? I think it is cohort but I got confused by looking on the internet for the difference $\endgroup$
    – user19745561
    Feb 24, 2023 at 12:44
  • $\begingroup$ If it is reasonable to consider that the outcome of interest (death) is rare in the population (i.e. the universe, which is most likely not known in your case) behind the study, then the odds ratio is similar to the risk ratio. More details with a clarifying example here: en.wikipedia.org/wiki/Odds_ratio $\endgroup$
    – mastropi
    Feb 24, 2023 at 15:52

2 Answers 2

3
$\begingroup$

This is fairly easily done with a marginal effect in R.

However, we should warn you that the design of the study is incredibly important. If the design is a case control (where cases are purposefully oversampled) then the relative risk calculation is biased precisely because the frequency of cases was biased (by design).

Let's assume you ran a cohort study and hence relative risks are indeed allowable. The first thing to understand is that logistic regression can predict the risk conditional on age and sex. Let age be represented by $x$ and sex by $w$. Your model is

$$ p(x, w) = \dfrac{1}{1 + \exp(-(\beta_0 + \beta_1x + \beta_2w))} $$

Here, the $\beta$ are the log odds ratios. So given someone who is $\delta$ years older than some reference patient, the relative risk is

$$ \dfrac{p(x+\delta, w)}{p(x, w)} $$

Note that the relative risk is then going to depend on the denominator, so there will not be a single relative risk, there will be a distribution of them, and the relative risk depends on the combination of age and sex. However, we can calculate an average marginal effect and report a confidence interval for that. Here is how we would do that in R (you will need to install the {marginaleffects} package).

library(tidyverse)
library(marginaleffects)
#> Warning: package 'marginaleffects' was built under R version 4.2.2
set.seed(0)
N <- 250
# Imagine a rescaled age variable
age <- rnorm(N, 0, 1)
sex <- rbinom(N, 1, 0.49)

p <- plogis(-2 + 0.2*age + 0.1*sex)
y <- rbinom(N, 1, p)

# Your model
fit <- glm(y ~ age + sex, family = binomial())


avg_comparisons(
  model=fit,
  # This next part computes the relative risk
  # Else, the risk difference is returned.
  transform_pre = 'lnratioavg',
  transform_post = exp
)
#> 
#>  Term              Contrast Estimate Pr(>|z|) 2.5 % 97.5 %
#>   age mean(+1)                  1.22    0.166 0.920   1.62
#>   sex ln(mean(1) / mean(0))     1.24    0.461 0.701   2.19
#> 
#> Prediction type:  response 
#> Columns: type, term, contrast, estimate, p.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo

Created on 2023-02-24 by the reprex package (v2.0.1)

It is worth breaking this down. What this output is telling me is that when I increase age by 1 unit (here 1 standard deviation since I've used a rescaled age with 0 mean and standard deviation 1) then the relative risk is 1.22 (or 22% increase to the risk). The 95% CI is also provided. However, that is the AVERAGE. Remember, the relative risk depends on your age and your sex.

Here is a histogram of estimated relative risks for each patient in these data.

comparisons(
  fit,
  variables = 'age',
  transform_pre = 'lnratio',
  transform_post = exp
) %>% 
  as.data.frame() %>% 
  ggplot(aes(estimate)) + 
  geom_histogram()

enter image description here

We can also estimate the relative risk conditional on age and sex at the same time like this

comparisons(
  fit,
  newdata = datagrid(age=-seq(-3, 3, 0.1), sex=0:1),
  variables = 'age',
  transform_pre = 'lnratio',
  transform_post = exp
) %>% 
  as.data.frame() %>% 
  ggplot(aes(age, estimate, color=factor(sex))) + 
  geom_line()

So given the age of a patient and their sex, you can report the relative risk associated with (in this case) a 1 year increase to their age.

enter image description here

Its worth repeating that the validity of these approaches depends almost entirely on the design of the study. All of this is useless if the study is a case control since the baseline risk of the outcome is biased.

$\endgroup$
5
  • 1
    $\begingroup$ Very useful post! Two minor points. By default, the p values are computed for a null of 0, so it can make sense to set hypothesis=1. Also, the as.data.frame() calls are superfluous, because the marginaleffects objects are already data frames. $\endgroup$
    – Vincent
    Mar 2, 2023 at 11:55
  • $\begingroup$ @Vincent Thanks for the tips. regarding the hypothesis comment, are you reffering to the initial code block? Would you do avg_comparisons(model=fit, transform_pre = 'lnratioavg', transform_post = exp, hypothesis = 1)? $\endgroup$ Mar 2, 2023 at 12:54
  • $\begingroup$ Huh, maybe I wrote that too fast. The tricky thing here is that the p value is computed for a null of 0 on the scale before applying transform_post. Maybe marginaleffects should only report the confidence interval in those cases to avoid confusion. What do you think? $\endgroup$
    – Vincent
    Mar 2, 2023 at 13:38
  • $\begingroup$ @Vincent I actually think the behaviour of marginal effects is fine. Most other approaches for computing the relative risk compute the p value on the log risk scale, so I assumed this was the case too. It also made sense comparing p values to estimated relative risk confidence intervals. Don't change it IMO $\endgroup$ Mar 2, 2023 at 13:55
  • $\begingroup$ makes sense. thanks for the input. $\endgroup$
    – Vincent
    Mar 2, 2023 at 14:00
0
$\begingroup$

Estimation of relative risks for binomial outcome can also be done through a modified Poisson model, especially useful when the outcome is common (>10%, let's say). It basically combines a log Poisson regression with a robust variance (sandwich estimator from general estimating equation). Here is a example using Titanic data:

library(titanic)
library(geepack)
library(marginaleffects)

dat<- na.omit(as.data.frame(titanic_train))

fit1<- summary(geeglm(as.numeric(Survived)~as.factor(Pclass)+Sex+Age,data=dat,family=poisson("log"),id=PassengerId,corstr="independence"))
cbind(RiskRatio = exp(coef(fit1)))

 
 #                 RiskRatio.Estimate RiskRatio.Std.err RiskRatio.Wald RiskRatio.Pr(>|W|)
#(Intercept)                     1.648              1.11       4.79e+09                  1
#as.factor(Pclass)2              0.689              1.08       9.51e+10                  1
#as.factor(Pclass)3              0.379              1.12       1.22e+34                  1
#Sexmale                         0.323              1.10       2.22e+58                  1
#Age                             0.987              1.00       3.58e+08                  1

Complete code can be found in Gallis and al and other comparison with odds ratios here. Interestingly, they also compare with marginal standardization that utilizes output from logit regression, the method described here above.

$\endgroup$

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.