# 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)