# Simple Emax model fit with Stan

#### Kenta Yoshida

#### 2019-11-19

    library(rstanemax)
    #> Warning: package 'Rcpp' was built under R version 3.5.3
    library(dplyr)
    #> Warning: package 'dplyr' was built under R version 3.5.3
    library(ggplot2)
    #> Warning: package 'ggplot2' was built under R version 3.5.3
    set.seed(12345)

This vignette provide an overview of the workflow of Emax model analysis using this package.

# Typical workflow

## Model run with `stan_emax` function

`stan_emax()` is the main function of this package to perform Emax model analysis on the data. This function requires minimum two input arguments - `formula` and `data`. In the `formula` argument, you will specify which columns of `data` will be used as exposure and response data, in a format similar to `stats::lm()` function, e.g. `response ~ exposure`.

    data(exposure.response.sample)

    fit.emax <- stan_emax(response ~ exposure, data = exposure.response.sample,
                          # the next line is only to make the example go fast enough
                          chains = 2, iter = 1000, seed = 12345)

    fit.emax
    #> ---- Emax model fit with rstanemax ----
    #>
    #>        mean se_mean    sd  2.5%   25%   50%   75%  97.5%  n_eff Rhat
    #> emax  92.32    0.31  6.04 81.25 88.17 92.27 96.35 103.96 384.95 1.01
    #> e0     5.63    0.22  4.67 -2.97  2.42  5.78  8.85  14.39 444.72 1.01
    #> ec50  76.11    0.87 21.05 43.81 61.75 73.35 86.40 124.79 587.83 1.00
    #> gamma  1.00     NaN  0.00  1.00  1.00  1.00  1.00   1.00    NaN  NaN
    #> sigma 16.65    0.06  1.56 14.12 15.56 16.50 17.48  20.37 656.97 1.01
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `plot()` function to visualize model fit
    #> * Use `posterior_predict()` or `posterior_predict_quantile()` function to get
    #>   raw predictions or make predictions on new data
    #> * Use `extract_obs_mod_frame()` function to extract raw data
    #>   in a processed format (useful for plotting)

`plot()` function shows the estimated Emax model curve with 95% credible intervals of parameters.

    plot(fit.emax)

Raw output from `rstan` is stored in the output variable, and you can access it with `extract_stanfit()` function.

    class(extract_stanfit(fit.emax))
    #> [1] "stanfit"
    #> attr(,"package")
    #> [1] "rstan"

## Prediction of response with new exposure data

`posterior_predict()` function allows users to predict the response using new exposure data. If `newdata` is not provided, the function returns the prediction on the exposures in original data. The default output is a matrix of posterior predictions, but you can also specify “dataframe” or “tibble” that contain posterior predictions in a long format. See help of `rstanemax::posterior_predict()` for the description of two predictions, `respHat` and `response`.

    response.pred <- posterior_predict(fit.emax, newdata = c(0, 100, 1000), returnType = "tibble")

    response.pred %>% select(mcmcid, exposure, respHat, response)
    #> # A tibble: 3,000 x 4
    #>    mcmcid exposure respHat response
    #>     <int>    <dbl>   <dbl>    <dbl>
    #>  1      1        0    5.89    -7.17
    #>  2      1      100   52.0     57.7
    #>  3      1     1000   86.9     91.8
    #>  4      2        0   -4.11   -22.1
    #>  5      2      100   61.8     57.0
    #>  6      2     1000   93.0     73.2
    #>  7      3        0   -4.90   -50.4
    #>  8      3      100   66.2     55.4
    #>  9      3     1000  100.     109.
    #> 10      4        0    5.95   -10.7
    #> # ... with 2,990 more rows

You can also get quantiles of predictions with `posterior_predict_quantile()` function.

    resp.pred.quantile <- posterior_predict_quantile(fit.emax, newdata = seq(0, 5000, by = 100))
    resp.pred.quantile
    #> # A tibble: 51 x 11
    #>    exposure covemax covec50 cove0 Covariates ci_low ci_med ci_high pi_low pi_med
    #>       <dbl> <fct>   <fct>   <fct> <chr>       <dbl>  <dbl>   <dbl>  <dbl>  <dbl>
    #>  1        0 1       1       1     ""          -2.09   5.78    13.0  -22.8   6.60
    #>  2      100 1       1       1     ""          52.3   58.8     64.5   30.4  59.5
    #>  3      200 1       1       1     ""          67.8   72.9     77.7   45.3  72.9
    #>  4      300 1       1       1     ""          75.0   79.5     84.2   52.9  79.0
    #>  5      400 1       1       1     ""          78.8   83.4     88.1   56.7  83.2
    #>  6      500 1       1       1     ""          81.2   85.9     90.6   60.3  87.0
    #>  7      600 1       1       1     ""          82.6   87.6     92.6   57.5  87.2
    #>  8      700 1       1       1     ""          83.7   88.9     94.1   61.5  89.3
    #>  9      800 1       1       1     ""          84.6   89.9     95.4   63.2  89.9
    #> 10      900 1       1       1     ""          85.2   90.7     96.4   60.5  90.3
    #> # ... with 41 more rows, and 1 more variable: pi_high <dbl>

Input data can be obtained in a same format with `extract_obs_mod_frame()` function.

    obs.formatted <- extract_obs_mod_frame(fit.emax)

These are particularly useful when you want to plot the estimated Emax curve.

    ggplot(resp.pred.quantile, aes(exposure, ci_med)) +
      geom_line() +
      geom_ribbon(aes(ymin=ci_low, ymax=ci_high), alpha = .5) +
      geom_ribbon(aes(ymin=pi_low, ymax=pi_high), alpha = .2) +
      geom_point(data = obs.formatted,
                 aes(y = response)) +
      labs(y = "response")

# Fix parameter values in Emax model

You can fix parameter values in Emax model for E0 and/or gamma (Hill coefficient). See help of `stan_emax()` for the details. The default is to fix gamma at 1 and to estimate E0 from data.

Below is the example of estimating gamma from data.

    data(exposure.response.sample)

    fit.emax.sigmoidal <- stan_emax(response ~ exposure, data = exposure.response.sample, gamma.fix = NULL,
                                    # the next line is only to make the example go fast enough
                                    chains = 2, iter = 1000, seed = 12345)
    #> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
    #> Running the chains for more iterations may help. See
    #> http://mc-stan.org/misc/warnings.html#tail-ess

    fit.emax.sigmoidal
    #> ---- Emax model fit with rstanemax ----
    #>
    #>        mean se_mean    sd  2.5%   25%   50%   75%  97.5%  n_eff Rhat
    #> emax  91.77    0.99 12.29 73.44 84.14 89.72 96.48 123.18 154.54 1.00
    #> e0     6.43    0.23  4.94 -3.17  3.00  6.52  9.50  15.63 446.65 1.00
    #> ec50  83.26    3.79 46.60 44.55 61.03 72.14 87.97 223.54 151.51 1.01
    #> gamma  1.12    0.02  0.34  0.53  0.90  1.09  1.31   1.88 322.49 1.00
    #> sigma 16.74    0.06  1.66 13.93 15.56 16.60 17.73  20.31 735.46 1.00
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `plot()` function to visualize model fit
    #> * Use `posterior_predict()` or `posterior_predict_quantile()` function to get
    #>   raw predictions or make predictions on new data
    #> * Use `extract_obs_mod_frame()` function to extract raw data
    #>   in a processed format (useful for plotting)

You can compare the difference of posterior predictions between two models (in this case they are very close to each other):

    exposure_pred <- seq(min(exposure.response.sample$exposure),
                         max(exposure.response.sample$exposure),
                         length.out = 100)

    pred1 <-
      posterior_predict_quantile(fit.emax, exposure_pred) %>%
      mutate(model = "Emax")
    pred2 <-
      posterior_predict_quantile(fit.emax.sigmoidal, exposure_pred) %>%
      mutate(model = "Sigmoidal Emax")

    pred <- bind_rows(pred1, pred2)

    ggplot(pred, aes(exposure, ci_med, color = model, fill = model)) +
      geom_line() +
      geom_ribbon(aes(ymin=ci_low, ymax=ci_high), alpha = .3) +
      geom_ribbon(aes(ymin=pi_low, ymax=pi_high), alpha = .1, color = NA) +
      geom_point(data=exposure.response.sample, aes(exposure, response),
                 color = "black", fill = NA, size=2) +
      labs(y = "response")

# Set covariates

You can specify categorical covariates for Emax, EC50, and E0. See help of `stan_emax()` for the details.

    data(exposure.response.sample.test)

    test.data <-
      mutate(exposure.response.sample.test,
             SEX = ifelse(cov2 == 0, "MALE", "FEMALE"))

    fit.cov <- stan_emax(formula = resp ~ conc, data = test.data,
                         param.cov = list(emax = "SEX"),
                         # the next line is only to make the example go fast enough
                         chains = 2, iter = 1000, seed = 12345)

    fit.cov
    #> ---- Emax model fit with rstanemax ----
    #>
    #>                mean se_mean    sd  2.5%   25%    50%    75%  97.5%  n_eff Rhat
    #> emax[FEMALE]  81.34    0.15  3.89 73.91 78.82  81.13  83.78  89.05 651.49    1
    #> emax[MALE]    87.66    0.20  5.00 78.19 84.22  87.72  91.15  97.42 621.72    1
    #> e0            15.07    0.09  2.35 10.30 13.51  15.13  16.70  19.32 638.51    1
    #> ec50         106.28    0.74 20.07 70.69 93.03 104.42 116.83 150.64 740.38    1
    #> gamma          1.00     NaN  0.00  1.00  1.00   1.00   1.00   1.00    NaN  NaN
    #> sigma         10.53    0.04  1.06  8.77  9.80  10.42  11.16  12.82 573.60    1
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `plot()` function to visualize model fit
    #> * Use `posterior_predict()` or `posterior_predict_quantile()` function to get
    #>   raw predictions or make predictions on new data
    #> * Use `extract_obs_mod_frame()` function to extract raw data
    #>   in a processed format (useful for plotting)
    plot(fit.cov)