# Relative risk models

    library(ameras)
    library(ggplot2)
    data(data, package="ameras")

# 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. See Effect modification for details on reference-plus-contrast and subgroup-specific modifier coding.

# 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 supplying `modifier`, respectively). This model is fit by setting `model="EXP"` as follows:

    fit.ameras.exp <- ameras(Y.binomial~dose(V1:V10, deg=2, model="EXP")+X1+X2,
                             data=data, family="binomial", methods="RC")
    #> Fitting RC
    summary(fit.ameras.exp)
    #> Call:
    #> ameras(formula = Y.binomial ~ dose(V1:V10, deg = 2, model = "EXP") +
    #>     X1 + X2, data = data, family = "binomial", methods = "RC")
    #>
    #> Rows: 3000
    #>
    #> Total CPU runtime: 0.2 seconds
    #>
    #> CPU runtime in seconds by method:
    #>
    #>  Method   Fit  CI Total
    #>      RC 0.173 0.0 0.173
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term Estimate      SE
    #>      RC  (Intercept) -0.94461 0.08409
    #>      RC           X1  0.44552 0.07667
    #>      RC           X2 -0.33376 0.09601
    #>      RC         dose  0.37904 0.10388
    #>      RC dose_squared  0.01943 0.02750
    #>
    #> Note: confidence intervals not yet computed. Use confint() to add them.

# 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. The degree is specified with `deg`, as in the exponential model; the example below uses `deg=2`. This model is fit by setting `model="ERR"` as follows:

    fit.ameras.err <- ameras(Y.binomial~dose(V1:V10, deg=2, model="ERR")+X1+X2,
                             data=data, family="binomial", methods="RC")
    #> Fitting RC
    summary(fit.ameras.err)
    #> Call:
    #> ameras(formula = Y.binomial ~ dose(V1:V10, deg = 2, model = "ERR") +
    #>     X1 + X2, data = data, family = "binomial", methods = "RC")
    #>
    #> Rows: 3000
    #>
    #> Total CPU runtime: 0.3 seconds
    #>
    #> CPU runtime in seconds by method:
    #>
    #>  Method   Fit  CI Total
    #>      RC 0.272 0.0 0.272
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method         Term Estimate      SE
    #>      RC  (Intercept) -0.87359 0.09759
    #>      RC           X1  0.44587 0.07672
    #>      RC           X2 -0.33552 0.09610
    #>      RC         dose  0.04878 0.21283
    #>      RC dose_squared  0.28763 0.08100
    #>
    #> Note: confidence intervals not yet computed. Use confint() to add them.

# 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 `model="LINEXP"` as follows:

    fit.ameras.linexp <- ameras(Y.binomial~dose(V1:V10, model="LINEXP")+X1+X2,
                             data=data, family="binomial", methods="RC")
    #> Fitting RC
    summary(fit.ameras.linexp)
    #> Call:
    #> ameras(formula = Y.binomial ~ dose(V1:V10, model = "LINEXP") +
    #>     X1 + X2, data = data, family = "binomial", methods = "RC")
    #>
    #> Rows: 3000
    #>
    #> Total CPU runtime: 0.3 seconds
    #>
    #> CPU runtime in seconds by method:
    #>
    #>  Method   Fit  CI Total
    #>      RC 0.254 0.0 0.254
    #>
    #> Summary of coefficients by method:
    #>
    #>  Method             Term Estimate      SE
    #>      RC      (Intercept)  -0.9326 0.08592
    #>      RC               X1   0.4456 0.07668
    #>      RC               X2  -0.3343 0.09603
    #>      RC      dose_linear   0.3255 0.11919
    #>      RC dose_exponential   0.3455 0.10814
    #>
    #> Note: confidence intervals not yet computed. Use confint() to add them.

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