In this section, we will show a basic workflow for performing an Emax model analysis for continuous endpoint.
Setup and load
Show the code
library (tidyverse)library (BayesERtools)library (posterior)library (tidybayes)library (here)library (gt)theme_set (theme_bw (base_size = 12 ))
Load data
Show the code
<- read_csv (here ("data" , "d_example_emax_nocov.csv" ))|> head () |> gt () |> fmt_number (decimals = 2 )
47.51
15.14
12.44
7.26
14.90
8.78
13.62
7.05
29.16
10.82
45.16
13.27
Show the code
ggplot (d_example_emax_nocov, aes (Conc, Y)) + geom_point () + geom_smooth (formula = y ~ x, method = "loess" , se = F, col = "darkgrey" )
Sigmoidal Emax model
Show the code
set.seed (1234 )<- dev_ermod_emax (data = d_example_emax_nocov,var_resp = "Y" ,var_exposure = "Conc" ,gamma_fix = NULL
── Emax model ──────────────────────────────────────────────────────────────────
ℹ Use `plot_er()` to visualize ER curve
── Developed model ──
---- Emax model fit with rstanemax ----
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
emax 15.92 0.24 5.49 8.67 11.64 14.71 19.07 29.38 545.77 1.01
e0 3.13 0.09 2.00 -2.15 2.25 3.58 4.54 5.62 550.23 1.00
ec50 24.94 0.30 8.63 14.98 19.47 22.51 28.06 48.31 843.12 1.00
gamma 1.77 0.03 0.77 0.75 1.22 1.62 2.19 3.54 664.61 1.00
sigma 0.86 0.00 0.16 0.62 0.75 0.84 0.95 1.24 1348.22 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)
Observation vs model fit
Show the code
<- sim_er (ermod_sigemax, output_type = c ("median_qi" ))plot_er_gof (d_sim_ermod_sigemax)
Parameter estimates
Show the code
<- as_draws_df (ermod_sigemax)<- summarize_draws (d_draws_sigemax)<- |> filter (variable == "ec50" ) |> pull (mean)|> gt () |> fmt_number (decimals = 2 )
ec50
24.94
22.51
8.63
5.65
16.17
42.09
1.00
1,111.48
1,157.42
sigma
0.86
0.84
0.16
0.15
0.64
1.15
1.00
1,434.70
1,777.17
gamma
1.77
1.62
0.77
0.68
0.84
3.12
1.00
593.05
951.82
e0
3.13
3.58
2.00
1.62
−0.84
5.38
1.00
675.84
809.63
emax
15.92
14.71
5.49
5.13
9.35
26.51
1.00
583.13
809.81
Prediction at a certain concentrations
Show the code
<- sim_er_new_exp (ermod_sigemax,exposure_to_sim_vec = c (10 , 20 , 30 , 50 ),output_type = c ("median_qi" ))|> select (- starts_with (".linpred" ), - c (.row, .width, .point, .interval)) |> gt () |> fmt_number (decimals = 2 ) |> tab_header (title = md ("Prediction at specific concentrations" )
Conc
.epred
.epred.lower
.epred.upper
.prediction
.prediction.lower
.prediction.upper
10.00
6.57
5.95
7.16
6.56
4.77
8.40
20.00
10.04
9.39
10.77
10.04
8.24
11.97
30.00
12.38
11.78
13.02
12.38
10.47
14.34
50.00
14.76
13.55
15.93
14.75
12.50
16.78
Show the code
|> ggplot (aes (x = Conc, y = Y)) + # Emax model curve geom_vline (xintercept = ec50_mean, linetype = "dashed" , color = "deepskyblue" ) + geom_ribbon (aes (y = .epred, ymin = .epred.lower, ymax = .epred.upper),alpha = 0.3 , fill = "deepskyblue" ) + geom_line (aes (y = .epred), color = "deepskyblue" ) + # Observed data geom_point (data = d_example_emax_nocov, color = "grey" ) + # Prediction at the specified doses geom_point (data = d_sim_new_conc, aes (y = .epred), color = "tomato" , size = 3 ) + geom_errorbar (data = d_sim_new_conc,aes (y = .epred, ymin = .epred.lower, ymax = .epred.upper),width = 1 , color = "tomato" ) + coord_cartesian (ylim = c (- 1 , 17 )) + scale_fill_brewer (palette = "Greys" ) + labs (title = "Sigmoidal E~max~ model predictions at new exposure levels" ,caption = "vertical dashed line: estimated EC~50~ value<br>area: 95% credible interval" ) + theme (plot.title = ggtext:: element_markdown (),plot.caption = ggtext:: element_markdown ())
