Relative Excess Risk Due to Interaction (RERI) has been used to quantify the joint effects of 2 exposures in epidemiology. RERI is the proportion of disease among those with both exposures that is attributable to their interaction.

Given the example dataset from Nie et al, 2010:

enter image description here

This formula calculates the RERI:

$RERI = RR(AB) - RR(\bar{A}B) - RR(A\bar{B}) + 1$


"Following the notation of Hosmer and Lemenshow, let A and B denote the presence of, and $\bar{A}$ and $\bar{B}$ denote the absence of, 2 binary exposures. Let the quantity RR denote the relative risk".

My question is:

How can I implement in practice in R the estimation of the confidence interval for RERI using:

  1. bootstrap approaches with a continuity correction as recommended by Nie et al, 2010, or
  2. likelihood-based approaches, or
  3. delta method?

Any input would be much appreciated.

  • $\begingroup$ Thanks, @Heteroskedastic. Does not need to be bootstrap necessarily. Implementation of any of these approaches will be much appreciated. $\endgroup$ – Krantz Dec 7 '18 at 17:05

This is simply the interaction effect as far as I can tell. Data such as this can be analyzed as binomial since we have successes out of number of trials, so logistic regression is a reasonable option. The problem is logistic regression provides us with odds ratios not risk ratios.

However, if the risk is low in all cells, the odds ratio is a good enough approximation to the risk ratio. And this is helpful because it is easier to estimate odds ratios and conduct inference on them. Inference is the more difficult part with risk ratios.

And in relation to inference, the delta method results in the same solution as Wald confidence intervals. And profile likelihood confidence intervals are readily obtainable from a logistic regression model.

The complicated part is that we have to parameterize the model such that the estimated parameters result in the excess interaction. The standard parameterization for GLMs (and design matrix) uses reference group coding such that the estimated parameters are deviations from the reference group. Of course, we can compute these values from the GLM, but we lose the ability to compute profile likelihood CI automatically in software for derived quantities.

Richardson & Kaufman (2009) provide a reparameterization of the model that helps:

$$\frac{\pi}{1-\pi}=e^{\beta_0}\times \big(1 + \beta_1A + \beta_2B + \beta_3(A\times B)\big)$$

where $\pi/(1-\pi)$ are the odds of success $(\tau)$, and solving for $\pi$, we get $\pi=\tau/(1+\tau)$. And $\beta_3$ in particular is the parameter OP cares for. Since this is a custom parameterization, we can write the likelihood function, supply it to the binomial density function and maximize the log-likelihood for the parameters $\{\beta_0, \beta_1, \beta_2, \beta_3\}$.

Here is a pass at the problem in R:

(dat <- data.frame(
  A = c(1, 1, 0, 0), B = c(1, 0, 1, 0),
  Pass = c(278, 100, 153, 79), Fail = c(743, 581, 1232, 1731)
#   A B Pass Fail
# 1 1 1  278  743
# 2 1 0  100  581
# 3 0 1  153 1232
# 4 0 0   79 1731

Next to the likelihood function, precisely negative log-likelihood because optimizers have an easier time finding the minimum instead of the maximum of functions:

ll <- function (b0, b1, b2, b3) {
  odds <- exp(b0) * (1 + b1 * dat$A + b2 * dat$B + b3 * dat$A * dat$B)
    x = dat$Pass, size = rowSums(dat[, 3:4]),
    prob = odds / (1 + odds), log = TRUE))

Next, we use the bbmle package for maximum likelihood estimation, it will provide profile likelihood CIs.

# Call mle, provide starting values
(fit <- mle2(ll, start = list(b0 = 0, b1 = 1, b2 = 1, b3 = 1)))
# Coefficients:
#        b0        b1        b2        b3 
# -3.087007  2.771328  1.721142  2.705877 

You will probably get some warning messages. So $\beta_3$ is the relative excess odds ratio for AB.

One can arrive at the final value using standard logistic regression and the formula in the question:

(coef.def <- unname(coef(glm(cbind(Pass, Fail) ~ A * B, binomial, dat))))
# [1] -3.0870067  1.3274261  1.0010505 -0.2245448
exp(sum(coef.def[2:4])) - exp(coef.def[2]) - exp(coef.def[3]) + 1
# [1] 2.705878

This is the same value from the custom parameterization, but since it is derived, we can not obtain profile likelihood CIs from R.

For Wald CI which are equivalent to the delta method:

cbind(fit@coef - qnorm(.975) * coef(summary(fit))[, 2],
      fit@coef + qnorm(.975) * coef(summary(fit))[, 2])[4, ]
#          [,1]      [,2]
# b3  1.1980204  4.213734

And the profile likelihood CI:

confint(fit)[4, ]
#        2.5 %    97.5 %
# b3  1.257277  4.356140

These are the Wald and profile likelihood CI respectively for OP's quantity of interest if you accept the odds ratio as an approximation to the risk ratio.

Richardson, D. B., & Kaufman, J. S. (2009). Estimation of the relative excess risk due to interaction and associated confidence bounds. American Journal of Epidemiology, 169(6), 756–760. https://doi.org/10.1093/aje/kwn411

EDIT, log link for RR

Try likelihood:

ll.rr <- function (b0, b1, b2, b3) {
  p <- exp(b0) * (1 + b1 * dat$A + b2 * dat$B + b3 * dat$A * dat$B)
  -sum(dbinom(dat$Pass, rowSums(dat[, 3:4]), p, TRUE))
# Use a highly negative intercept to start, should work.
coef(fit <- mle2(ll.rr, start = list(b0 = -4, b1 = 0, b2 = 0, b3 = 0)))
#        b0        b1        b2        b3 
# -3.131679  2.364507  1.531149  1.342961 
  • $\begingroup$ Thanks a lor fot this @Heteroskedastic Jim. Just one question. Just two questions: $\endgroup$ – Krantz Dec 7 '18 at 19:44
  • 1
    $\begingroup$ Well, note that I do not make use of glm here to arrive at the results. So I'm not sure how would work. $\endgroup$ – Heteroskedastic Jim Dec 7 '18 at 19:46
  • 1
    $\begingroup$ It is possible to reparameterize but if using the binomial distribution, the log link has no guarantee that it will produce probabilities under 0, leading to optimization failure. If using a pre-packaged function, it will work to use a log link since the predictors are categorical. There are some approaches like censoring the probabilities at 1 that will help with estimation, you can call punif() function. You can also initialize the optimization with a highly -ve intercept. With regard to weighting, it is possible to incorporate weights into the likelihood by mimicing whatever svyglm does. $\endgroup$ – Heteroskedastic Jim Dec 7 '18 at 20:40
  • 1
    $\begingroup$ I edited the question to help with RR. $\endgroup$ – Heteroskedastic Jim Dec 7 '18 at 20:52
  • 1
    $\begingroup$ You'd have to replicate then modify the svyglm likelihood for the specific types of survey weights. I cannot make that effort. $\endgroup$ – Heteroskedastic Jim Dec 7 '18 at 21:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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