Show the code
library(tidyverse)
library(BayesERtools)
library(posterior)
library(tidybayes)
library(bayesplot)
library(loo)
library(here)
library(gt)
theme_set(theme_bw(base_size = 12))9 Simulation from fitted model
This page showcase the model simulation using the Emax model with no covariate.
9.1 Load data
Show the code
d_example_emax_nocov <- read_csv(here("data", "d_example_emax_nocov.csv"))9.2 Fit models
Show the code
set.seed(1234)
ermod_sigemax <- dev_ermod_emax(
data = d_example_emax_nocov,
var_resp = "Y",
var_exposure = "Conc",
gamma_fix = NULL
)
ermod_sigemax
── 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)
Another model without sigmoidal component; will be used when we do model comparison.
Show the code
set.seed(1234)
ermod_emax <- dev_ermod_emax(
data = d_example_emax_nocov,
var_resp = "Y",
var_exposure = "Conc",
gamma_fix = 1
)
ermod_emax
── 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 22.52 0.08 2.46 18.73 20.84 22.17 23.81 28.09 976.89 1.00
e0 0.98 0.04 1.36 -2.23 0.32 1.18 1.91 2.98 912.27 1.00
ec50 31.28 0.39 11.65 13.08 22.74 29.70 38.10 57.39 875.72 1.00
gamma 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
sigma 0.89 0.00 0.17 0.63 0.77 0.86 0.98 1.28 1216.16 1.01
* 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)
10 Extrapolation
This represents uncertainty in the model parameters.
p(f(theta)|xnew, yobs)
Show the code
new_conc_vec <- seq(0, 100, by = 1)
ersim_sigemax <-
sim_er_new_exp(ermod_sigemax, new_conc_vec)
ersim_sigemax_med_qi <-
ersim_sigemax |>
calc_ersim_med_qi(qi_width = c(0.5, 0.95))
ggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .epred)) +
geom_ribbon(
data = ersim_sigemax_med_qi |> filter(.width == 0.95),
aes(ymin = .epred.lower, ymax = .epred.upper),
fill = "yellow3", alpha = 0.5
) +
geom_ribbon(
data = ersim_sigemax_med_qi |> filter(.width == 0.5),
aes(ymin = .epred.lower, ymax = .epred.upper),
fill = "orange1"
) +
geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.5), col = "darkred") +
geom_point(data = d_example_emax_nocov, aes(y = Y)) +
coord_cartesian(ylim = c(-1, 22)) +
scale_y_continuous("Response", breaks = 5 * 0:4) +
labs(x = "Exposure",
title = "Extrapolation (no residual)",
subtitle = "Represents uncertainty in the model parameters (Credible interval)",
caption = "95% CI in yellow, 50% CI in orange")
This represents uncertainty in the model parameters plus the residual error.
p(ynew|xnew, yobs)
Show the code
ggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .prediction)) +
geom_ribbon(
data = ersim_sigemax_med_qi |> filter(.width == 0.95),
aes(ymin = .prediction.lower, ymax = .prediction.upper),
fill = "yellow3", alpha = 0.5
) +
geom_ribbon(
data = ersim_sigemax_med_qi |> filter(.width == 0.5),
aes(ymin = .prediction.lower, ymax = .prediction.upper),
fill = "orange1"
) +
geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.5), col = "darkred") +
geom_point(data = d_example_emax_nocov, aes(y = Y)) +
coord_cartesian(ylim = c(-1, 22)) +
scale_y_continuous("Response", breaks = 5 * 0:4) +
labs(x = "Exposure",
title = "Extrapolation (incl. residual)",
subtitle = "Represents uncertainty + residual error (Prediction interval)",
caption = "95% PI in yellow, 50% PI in orange")
No discernible difference between the two models.
Show the code
ersim_emax <-
sim_er_new_exp(ermod_emax, new_conc_vec)
ersim_emax_med_qi <-
ersim_emax |>
calc_ersim_med_qi(qi_width = c(0.5, 0.95))
ggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .epred)) +
geom_ribbon(
data = ersim_sigemax_med_qi |> filter(.width == 0.95),
aes(ymin = .epred.lower, ymax = .epred.upper),
fill = "orange1", alpha = 0.5
) +
geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.95), col = "darkred") +
geom_ribbon(
data = ersim_emax_med_qi |> filter(.width == 0.95),
aes(ymin = .epred.lower, ymax = .epred.upper),
fill = "turquoise3", alpha = 0.5
) +
geom_line(data = ersim_emax_med_qi |> filter(.width == 0.95), col = "steelblue3") +
geom_point(data = d_example_emax_nocov, aes(y = Y)) +
coord_cartesian(ylim = c(-1, 22)) +
scale_y_continuous("Response", breaks = 5 * 0:4) +
labs(x = "Exposure",
title = "Extrapolation (no residual)",
subtitle = "Represents uncertainty in the model parameters (Credible interval)",
caption = "Sigmoidal Emax: Orange, Emax: Blue, 50% CI")