# Simple Emax model fit with Stan

#### Kenta Yoshida

#### 2026-06-12

    library(rstanemax)
    library(dplyr)
    library(ggplot2)
    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.05    0.22  6.04 80.16 87.90 92.15 96.12 103.97 758.47    1
    #> e0     5.68    0.18  4.79 -4.17  2.56  5.61  8.91  15.17 736.90    1
    #> ec50  76.28    0.83 22.46 42.90 60.92 72.93 87.67 125.76 738.39    1
    #> gamma  1.00     NaN  0.00  1.00  1.00  1.00  1.00   1.00    NaN  NaN
    #> sigma 16.57    0.06  1.65 13.71 15.44 16.44 17.68  20.03 654.58    1
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `extract_param()` function to extract posterior draws of key parameters
    #> * 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)

Output of `plot()` function (for `stanemax` object) is a `ggplot` object, so you can apply additional settings as you would do in `ggplot2`.
Here is an example of using log scale for x axis (note that exposure == 0 is hanging at the very left, making the curve a bit weird).

    plot(fit.emax) + scale_x_log10() + expand_limits(x = 1)
    #> Warning in scale_x_log10(): log-10 transformation introduced infinite values.
    #> log-10 transformation introduced infinite values.
    #> log-10 transformation introduced infinite values.

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")
    #> Warning: The `returnType` argument of `posterior_predict()` is deprecated as of
    #> rstanemax 0.1.8.
    #> This warning is displayed once per session.
    #> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
    #> generated.

    response.pred %>% select(mcmcid, exposure, respHat, response)
    #> # A tibble: 3,000 × 4
    #>    mcmcid exposure   respHat response
    #>     <int>    <dbl> <dbl[1d]>    <dbl>
    #>  1      1        0      5.34     7.51
    #>  2      1      100     65.7     92.8
    #>  3      1     1000     93.8     86.8
    #>  4      2        0      8.20    -5.50
    #>  5      2      100     50.7     47.3
    #>  6      2     1000     90.3     83.1
    #>  7      3        0      1.41   -14.3
    #>  8      3      100     60.7     60.8
    #>  9      3     1000     92.7     85.5
    #> 10      4        0      4.98    40.1
    #> # ℹ 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 × 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.29   5.61    13.1  -24.2   5.75
    #>  2      100 1       1       1     ""          52.2   58.7     64.5   30.0  58.3
    #>  3      200 1       1       1     ""          67.6   72.7     77.4   46.5  72.8
    #>  4      300 1       1       1     ""          74.5   79.4     83.8   51.4  78.9
    #>  5      400 1       1       1     ""          78.5   83.1     87.7   57.1  84.0
    #>  6      500 1       1       1     ""          80.8   85.6     90.6   59.9  87.3
    #>  7      600 1       1       1     ""          82.2   87.4     92.6   60.6  87.3
    #>  8      700 1       1       1     ""          83.4   88.6     94.2   61.0  88.5
    #>  9      800 1       1       1     ""          84.2   89.6     95.3   61.9  88.8
    #> 10      900 1       1       1     ""          84.8   90.4     96.4   63.1  91.2
    #> # ℹ 41 more rows
    #> # ℹ 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")

Posterior draws of Emax model parameters can be extracted with `extract_param()` function.

    posterior.fit.emax <- extract_param(fit.emax)
    posterior.fit.emax
    #> # A tibble: 1,000 × 6
    #>    mcmcid  emax    e0  ec50     gamma     sigma
    #>     <int> <dbl> <dbl> <dbl> <dbl[1d]> <dbl[1d]>
    #>  1      1  93.3  5.34  54.7         1      16.8
    #>  2      2  91.6  8.20 115.          1      16.5
    #>  3      3  97.2  1.41  63.7         1      15.2
    #>  4      4  83.3  4.98  47.2         1      16.1
    #>  5      5  94.2  2.13  63.9         1      14.1
    #>  6      6  86.7  9.15  70.1         1      15.0
    #>  7      7 104.  -5.00  54.1         1      15.5
    #>  8      8  96.5  7.98  93.0         1      18.0
    #>  9      9  94.5  3.74  98.4         1      16.0
    #> 10     10  94.5  5.91  94.1         1      17.8
    #> # ℹ 990 more rows

# Fix parameter values in Emax model

You can fix parameter values in Emax model for Emax, 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 Emax and 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
    )

    fit.emax.sigmoidal
    #> ---- Emax model fit with rstanemax ----
    #>
    #>        mean se_mean    sd  2.5%   25%   50%   75%  97.5%  n_eff Rhat
    #> emax  88.81    0.51  9.35 73.14 82.60 87.95 93.93 109.34 337.17    1
    #> e0     6.98    0.22  4.86 -2.89  3.84  7.05 10.28  16.50 502.11    1
    #> ec50  75.20    1.64 26.84 43.68 58.93 69.93 83.93 138.88 267.26    1
    #> gamma  1.19    0.02  0.34  0.65  0.96  1.16  1.36   2.03 365.20    1
    #> sigma 16.69    0.07  1.67 13.80 15.45 16.58 17.77  20.26 661.42    1
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `extract_param()` function to extract posterior draws of key parameters
    #> * 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.with.cov)

    test.data <-
      mutate(exposure.response.sample.with.cov,
        SEX = ifelse(cov2 == "B0", "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.25    0.19  4.22 73.02 78.38  81.15  84.05  89.15 510.40 1.00
    #> emax[MALE]    87.36    0.23  5.15 77.35 83.76  87.46  90.82  97.83 519.75 1.00
    #> e0            15.17    0.11  2.43 10.29 13.70  15.27  16.70  19.66 504.33 1.00
    #> ec50         106.13    0.86 21.42 67.85 91.97 104.61 118.72 155.08 617.15 1.01
    #> gamma          1.00     NaN  0.00  1.00  1.00   1.00   1.00   1.00    NaN  NaN
    #> sigma         10.56    0.04  1.06  8.81  9.77  10.47  11.28  12.92 749.70 1.00
    #>
    #> * Use `extract_stanfit()` function to extract raw stanfit object
    #> * Use `extract_param()` function to extract posterior draws of key parameters
    #> * 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)

You can extract MCMC samples from raw stanfit and evaluate differences between groups.

    fit.cov.posterior <-
      extract_param(fit.cov)

    emax.posterior <-
      fit.cov.posterior %>%
      select(mcmcid, SEX, emax) %>%
      tidyr::pivot_wider(names_from = SEX, values_from = emax) %>%
      mutate(delta = FEMALE - MALE)

    ggplot2::qplot(delta, data = emax.posterior, bins = 30) +
      ggplot2::labs(x = "emax[FEMALE] - emax[MALE]")
    #> Warning: `qplot()` was deprecated in ggplot2 3.4.0.
    #> This warning is displayed once per session.
    #> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
    #> generated.

    # Credible interval of delta
    quantile(emax.posterior$delta, probs = c(0.025, 0.05, 0.5, 0.95, 0.975))
    #>       2.5%         5%        50%        95%      97.5%
    #> -14.643119 -13.342322  -6.111658   1.025680   2.219764

    # Posterior probability of emax[FEMALE] < emax[MALE]
    sum(emax.posterior$delta < 0) / nrow(emax.posterior)
    #> [1] 0.914