Show the code
library(tidyverse)
library(BayesERtools)
library(posterior)
library(ggdist)
library(bayesplot)
library(loo)
library(here)
library(gt)
theme_set(theme_bw(base_size = 12))9 Model diagnostics & comparison
This page showcase the model diagnosis and comparison for the Emax model
9.1 Load data
Show the code
d_sim_emax# A tibble: 300 × 9
dose exposure response_1 response_2 cnt_a cnt_b cnt_c bin_d bin_e
9.2 Fit models
Show the code
set.seed(1234)
ermod_sigemax <- dev_ermod_emax(
data = d_sim_emax,
var_resp = "response_1",
var_exposure = "exposure",
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
emax 11.70 0.29 5.84 4.65 7.35 10.19 14.52 27.20 414.18
e0 6.45 0.21 4.40 -5.11 4.58 7.58 9.63 11.59 436.99
ec50 4892.89 129.05 3858.03 649.97 2787.66 4537.64 6219.22 12022.48 893.78
gamma 1.07 0.02 0.51 0.34 0.70 0.96 1.34 2.32 575.85
sigma 1.27 0.00 0.05 1.18 1.24 1.27 1.31 1.38 1513.21
Rhat
emax 1.01
e0 1.01
ec50 1.00
gamma 1.01
sigma 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_sim_emax,
var_resp = "response_1",
var_exposure = "exposure",
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
emax 10.17 0.06 1.66 7.96 9.02 9.79 10.96 14.39 704.66
e0 7.24 0.07 1.95 2.56 6.20 7.58 8.65 9.98 697.52
ec50 4213.96 57.75 1720.01 1743.37 2943.96 3917.32 5202.32 8117.56 887.11
gamma 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN
sigma 1.27 0.00 0.05 1.17 1.24 1.27 1.31 1.39 1246.77
Rhat
emax 1.01
e0 1.01
ec50 1.01
gamma NaN
sigma 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)
9.3 Parameter summary table
Show the code
d_draws_sigemax_summary <-
summarize_draws(ermod_sigemax)
d_draws_sigemax_summary |>
gt() |>
fmt_number(decimals = 2)| variable | mean | median | sd | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|
| ec50 | 4,892.89 | 4,537.64 | 3,858.03 | 2,538.85 | 949.54 | 8,655.60 | 1.01 | 680.72 | 672.71 |
| sigma | 1.27 | 1.27 | 0.05 | 0.05 | 1.19 | 1.36 | 1.00 | 1,516.16 | 1,780.38 |
| gamma | 1.07 | 0.96 | 0.51 | 0.45 | 0.41 | 2.06 | 1.01 | 513.16 | 772.22 |
| e0 | 6.45 | 7.58 | 4.40 | 3.49 | −2.57 | 11.30 | 1.01 | 499.98 | 646.45 |
| emax | 11.70 | 10.19 | 5.84 | 4.85 | 5.10 | 23.64 | 1.01 | 466.31 | 761.51 |
Here is the example of highest-density continuous interval (HDCI) for the median of ED50. See here for more details.
Show the code
# HDCI of median.ED50
as_draws_df(ermod_sigemax) |>
median_hdci(ec50)Warning: Dropping 'draws_df' class as required metadata was removed.
Warning: Dropping 'draws_df' class as required metadata was removed.
Warning: Dropping 'draws_df' class as required metadata was removed.
# A tibble: 1 × 6
ec50 .lower .upper .width .point .interval
9.4 Fitted values
Fitted values without residual errors (i.e. PRED in NONMEM term) can be extracted with sim_er() function. .epred is the expected value prediction. See ?sim_er for detail.
Show the code
sim_er(ermod_sigemax) |> head()# A tibble: 6 × 16
dose exposure response_1 response_2 cnt_a cnt_b cnt_c bin_d bin_e .row
, .iteration , .draw ,
# .epred , .linpred , .prediction
You can specify output_type = "median_qi" to get median and quantile intervals of the prediction.
Show the code
ersim_sigemax_med_qi <-
sim_er(ermod_sigemax, output_type = "median_qi")
ersim_sigemax_med_qi |>
arrange(.row) |>
head() |>
gt(rownames_to_stub = TRUE) |>
fmt_number(decimals = 2, columns = -.row)| dose | exposure | response_1 | response_2 | cnt_a | cnt_b | cnt_c | bin_d | bin_e | .row | .epred | .epred.lower | .epred.upper | .linpred | .linpred.lower | .linpred.upper | .prediction | .prediction.lower | .prediction.upper | .width | .point | .interval | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 100.00 | 4,151.15 | 12.80 | 1.00 | 5.71 | 2.33 | 7.83 | 0.00 | 1.00 | 1 | 12.63 | 12.34 | 12.90 | 12.63 | 12.34 | 12.90 | 12.61 | 10.17 | 15.18 | 0.95 | median | qi |
| 2 | 100.00 | 8,067.26 | 14.58 | 1.00 | 4.92 | 4.66 | 6.74 | 1.00 | 1.00 | 2 | 14.14 | 13.94 | 14.35 | 14.14 | 13.94 | 14.35 | 14.15 | 11.72 | 16.64 | 0.95 | median | qi |
| 3 | 100.00 | 4,877.67 | 12.76 | 1.00 | 4.88 | 4.21 | 4.68 | 1.00 | 1.00 | 3 | 13.00 | 12.75 | 13.23 | 13.00 | 12.75 | 13.23 | 13.01 | 10.50 | 15.52 | 0.95 | median | qi |
| 4 | 100.00 | 9,712.93 | 16.55 | 1.00 | 8.42 | 6.56 | 1.29 | 0.00 | 1.00 | 4 | 14.53 | 14.32 | 14.75 | 14.53 | 14.32 | 14.75 | 14.53 | 11.96 | 16.99 | 0.95 | median | qi |
| 5 | 100.00 | 11,490.74 | 14.41 | 0.00 | 4.37 | 3.96 | 3.55 | 0.00 | 1.00 | 5 | 14.86 | 14.64 | 15.08 | 14.86 | 14.64 | 15.08 | 14.83 | 12.33 | 17.31 | 0.95 | median | qi |
| 6 | 100.00 | 2,451.55 | 12.60 | 1.00 | 8.69 | 7.60 | 3.64 | 0.00 | 0.00 | 6 | 11.46 | 10.86 | 12.07 | 11.46 | 10.86 | 12.07 | 11.51 | 8.92 | 14.00 | 0.95 | median | qi |
9.5 Diagnostic plots
We use the bayesplot package (Cheat sheet) to visualize the model fit.
9.5.1 Convergence
Good fit results in:
- Parameter distributions from MCMC chains should overlap
- Trace plots should not show any trend
Rhatclose to 1 (e.g. < 1.1)
Show the code
d_draws_sigemax <- as_draws_df(ermod_sigemax)
mcmc_dens_overlay(d_draws_sigemax)
mcmc_trace(d_draws_sigemax)
mcmc_rhat(rhat(ermod_sigemax$mod$stanfit))
9.5.2 Parameter estimates distribution
Show the code
mcmc_hist(d_draws_sigemax)`stat_bin()` using `bins = 30`. Pick better value `binwidth`.
Show the code
mcmc_pairs(d_draws_sigemax,
off_diag_args = list(size = 0.5, alpha = 0.25))
9.5.3 GOF plots
Show the code
plot_er(ermod_sigemax, show_orig_data = TRUE)
Show the code
# Preparation for diagnostic plots
ersim_sigemax_med_qi <- sim_er(ermod_sigemax, output_type = "median_qi")
ersim_sigemax_med_qi |>
ggplot(aes(x = .epred, y = response_1)) +
geom_smooth(method = "loess", formula = y ~ x) +
geom_abline(linetype = 2, color = "grey") +
geom_point() +
geom_errorbar(
mapping = aes(xmin = .epred.lower, xmax = .epred.upper),
orientation = "y"
) +
labs(
title = "Observed vs Predicted",
x = "Predicted",
y = "Observed",
caption = "Symbol: median and 95% credible interval"
)
Show the code
ersim_sigemax_w_resid <-
sim_er(ermod_sigemax) |>
mutate(.residual = response_1 - .epred) # Add residuals for plotting
ersim_sigemax_w_resid_med_qi <-
ersim_sigemax_w_resid |>
group_by(across(!c(
.chain, .iteration, .draw, .epred, .linpred, .prediction, .residual
))) |>
median_qi(.epred, .residual)
ersim_sigemax_w_resid_med_qi |>
ggplot(aes(x = .epred, y = .residual)) +
xlab("Predicted (linear)") +
ylab("Residuals") +
geom_smooth(method = "loess", formula = y ~ x) +
geom_point() +
geom_errorbar(
mapping = aes(ymin = .residual.lower, ymax = .residual.upper),
width = 0
) +
geom_hline(aes(yintercept = 2), lty = 2, colour = "grey70") +
geom_hline(aes(yintercept = -2), lty = 2, colour = "grey70")
Show the code
ersim_sigemax_w_resid_med_qi |>
ggplot(aes(sample = .residual)) +
geom_qq() +
geom_qq_line(colour = "steelblue", lty = 2, alpha = 0.4) +
coord_equal() +
labs(x = "Theoretical", y = "Sample")
9.6 Model comparison
You can perform model comparison based on expected log pointwise predictive density (ELPD). ELPD is the Bayesian leave-one-out estimate (see ?loo-glossary).
Higher ELPD is better, therefore Emax model with γ fixed to be 1 appears better. However, elpd_diff is tiny and smaller than se_diff (see here), therefore we can consider the difference to be not meaningful.
Show the code
loo_sigemax <- loo(ermod_sigemax)
loo_emax <- loo(ermod_emax)
loo_compare(list(sigemax = loo_sigemax, emax = loo_emax)) elpd_diff se_diff
emax 0.0 0.0
sigemax -0.8 0.4