Relative risk models

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")
#> Fitting RC
#> Obtaining profile likelihood CI for dose
#> Obtaining profile likelihood CI for dose_squared
#> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
#> P-value for dose_squared upper bound more than 0.005 away from 0.05, reducing
#> tol.profCI and/or increasing maxit.profCI is recommended
#> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
#> P-value for dose_squared lower bound more than 0.005 away from 0.05, reducing
#> tol.profCI and/or increasing maxit.profCI is recommended
summary(fit.ameras.exp)
#> Call:
#> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars, 
#>     X = c("X1", "X2"), methods = "RC", deg = 2, doseRRmod = "EXP")
#> 
#> Total run time: 2.3 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     2.3
#> 
#> Summary of coefficients by method:
#> 
#>  Method         Term Estimate      SE CI.lowerbound CI.upperbound
#>      RC  (Intercept) -0.94461 0.08409            NA            NA
#>      RC           X1  0.44552 0.07667            NA            NA
#>      RC           X2 -0.33376 0.09601            NA            NA
#>      RC         dose  0.37904 0.10388       0.17336       0.57642
#>      RC dose_squared  0.01943 0.02750      -0.03282       0.07785

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")
#> Fitting RC
#> Obtaining profile likelihood CI for dose
#> Warning in ameras.rc(family = family, dosevars = dosevars, data = data, :
#> WARNING: Lower bound for dose is < 0 and may not exist if rescaling the
#> variable does not help
#> Obtaining profile likelihood CI for dose_squared
summary(fit.ameras.err)
#> Call:
#> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars, 
#>     X = c("X1", "X2"), methods = "RC", deg = 2, doseRRmod = "ERR")
#> 
#> Total run time: 3.4 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     3.4
#> 
#> Summary of coefficients by method:
#> 
#>  Method         Term Estimate      SE CI.lowerbound CI.upperbound
#>      RC  (Intercept) -0.87359 0.09759            NA            NA
#>      RC           X1  0.44587 0.07672            NA            NA
#>      RC           X2 -0.33552 0.09610            NA            NA
#>      RC         dose  0.04878 0.21283        0.0000        0.5115
#>      RC dose_squared  0.28763 0.08100        0.1325        0.4108

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")
#> Fitting RC
#> Obtaining profile likelihood CI for dose_linear
#> Obtaining profile likelihood CI for dose_exponential
summary(fit.ameras.linexp)
#> Call:
#> ameras(data = data, family = "binomial", Y = "Y.binomial", dosevars = dosevars, 
#>     X = c("X1", "X2"), methods = "RC", doseRRmod = "LINEXP")
#> 
#> Total run time: 4 seconds
#> 
#> Runtime in seconds by method:
#> 
#>  Method Runtime
#>      RC     4.0
#> 
#> Summary of coefficients by method:
#> 
#>  Method             Term Estimate      SE CI.lowerbound CI.upperbound
#>      RC      (Intercept)  -0.9326 0.08592            NA            NA
#>      RC               X1   0.4456 0.07668            NA            NA
#>      RC               X2  -0.3343 0.09603            NA            NA
#>      RC      dose_linear   0.3255 0.11919        0.1473        0.6339
#>      RC dose_exponential   0.3455 0.10814        0.1452        0.5770

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)