Relative risk models
In ameras: Analyze Multiple Exposure Realizations in Association Studies

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(ameras)
library(ggplot2)
data(data, package="ameras")
dosevars <- paste0("V", 1:10)

Introduction

For non-Gaussian families, three relative risk models for the main exposure are supported, the usual exponential model $$RR_i=\exp(\beta_1 D_i+\beta_2 D_i^2+ \mathbf{M}i^T \mathbf{\beta}{m1}D_i + \mathbf{M}i^T \mathbf{\beta}{m2} D_i^2),$$ the linear(-quadratic) excess relative risk (ERR) model $$RR_i= 1+\beta_1 D_i+\beta_2 D_i^2 + \mathbf{M}i^T \mathbf{\beta{m1}}D_i + \mathbf{M}i^T \mathbf{\beta}{m2}D_i^2,$$ and the linear-exponential model $$ RR_i= 1+(\beta_1 + \mathbf{M}i^T \mathbf{\beta}{m1}) D_i \exp{(\beta_2+ \mathbf{M}i^T \mathbf{\beta}{m2})D_i}. $$ This vignette illustrates fitting the three models using regression calibration for logistic regression, but the same syntax applies to all other settings.

Exponential relative risk

The usual exponential relative risk model is given by $RR_i=\exp(\beta_1 D_i+\beta_2 D_i^2+ \mathbf{M}i^T \mathbf{\beta}{m1}D_i + \mathbf{M}i^T \mathbf{\beta}{m2} D_i^2)$, where the quadratic and effect modification terms are optional (not fit by setting deg=1 and not passing anything to M, respectively). This model is fit by setting doseRRmod="EXP" as follows:

fit.ameras.exp <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="binomial", deg=2, doseRRmod = "EXP", methods="RC")
summary(fit.ameras.exp)

Linear excess relative risk

The linear excess relative risk model is given by $RR_i=1+\beta_1 D_i+\beta_2 D_i^2+ \mathbf{M}i^T \mathbf{\beta}{m1}D_i + \mathbf{M}i^T \mathbf{\beta}{m2} D_i^2$, where again the quadratic and effect modification terms are optional. This model is fit by setting doseRRmod="ERR" as follows:

fit.ameras.err <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="binomial", deg=2, doseRRmod = "ERR", methods="RC")
summary(fit.ameras.err)

Linear-exponential relative risk

The linear-exponential relative risk model is given by $RR_i= 1+(\beta_1 + \mathbf{M}i^T \mathbf{\beta}{m1}) D_i \exp{(\beta_2+ \mathbf{M}i^T \mathbf{\beta}{m2})D_i}$, where the effect modification terms are optional. This model is fit by setting doseRRmod="LINEXP" as follows:

fit.ameras.linexp <- ameras(Y="Y.binomial", dosevars=dosevars, X=c("X1","X2"), data=data, 
                            family="binomial", doseRRmod = "LINEXP", methods="RC")
summary(fit.ameras.linexp)

Comparison between models

To compare between models, it is easiest to do so visually:

ggplot(data.frame(x=c(0, 5)), aes(x))+
  theme_light(base_size=15)+
  xlab("Exposure")+
  ylab("Relative risk")+
  labs(col="Model", lty="Model") +
  theme(legend.position = "inside", 
        legend.position.inside = c(.2,.85),
        legend.box.background = element_rect(color = "black", fill = "white", linewidth = 1))+
  stat_function(aes(col="Linear-quadratic ERR", lty="Linear-quadratic ERR" ),fun=function(x){
    1+fit.ameras.err$RC$coefficients["dose"]*x + fit.ameras.err$RC$coefficients["dose_squared"]*x^2
  }, linewidth=1.2) + 
  stat_function(aes(col="Exponential", lty="Exponential"),fun=function(x){
    exp(fit.ameras.exp$RC$coefficients["dose"]*x + fit.ameras.exp$RC$coefficients["dose_squared"]*x^2)
  }, linewidth=1.2) +
  stat_function(aes(col="Linear-exponential", lty="Linear-exponential"),fun=function(x){
    1+fit.ameras.linexp$RC$coefficients["dose_linear"]*x * exp(fit.ameras.linexp$RC$coefficients["dose_exponential"]*x)
  }, linewidth=1.2)