GitHub actions to render site

.github/workflows/publish.yml

@@ -0,0 +1,25 @@
+on:
+  workflow_dispatch:
+  push:
+    branches: main
+
+name: Quarto Publish
+
+jobs:
+  build-deploy:
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+    steps:
+      - name: Check out repository
+        uses: actions/checkout@v4
+
+      - name: Set up Quarto
+        uses: quarto-dev/quarto-actions/setup@v2
+
+      - name: Render and Publish
+        uses: quarto-dev/quarto-actions/publish@v2
+        with:
+          target: gh-pages
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}

README.Rmd

@@ -132,7 +132,7 @@ tab_mod |>
 
 ## Important links
 
-- Rendered book: https://github.com/genentech/BayesERbook/
+- Rendered book: https://genentech.github.io/BayesERbook/
 
 ## Note for developer
 

README.md

@@ -95,7 +95,7 @@ install.packages('rstanemax')
 
 ## Important links
 
-- Rendered book: <https://github.com/genentech/BayesERbook/>
+- Rendered book: <https://genentech.github.io/BayesERbook/>
 
 ## Note for developer
 

_freeze/notebook/binary/mod_structure_comparison/execute-results/html.json

@@ -0,0 +1,17 @@
+{
+  "hash": "1ff2936cee198ac5f20ce8c1faa43eeb",
+  "result": {
+    "engine": "knitr",
+    "markdown": "# Model comparison between linear and E~max~\n\nThis page showcase how to compare model structures between linear and E~max~\nlogistic regression models\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(loo)\nlibrary(here)\nlibrary(gt)\n\ntheme_set(theme_bw(base_size = 12))\n```\n:::\n\n\n\n\n\n\n\n## Load data\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndata(d_sim_binom_cov)\n\nd_sim_binom_cov_2 <-\n  d_sim_binom_cov |>\n  mutate(\n    AUCss_1000 = AUCss / 1000, BAGE_10 = BAGE / 10,\n    BWT_10 = BWT / 10, BHBA1C_5 = BHBA1C / 5,\n    Dose = glue::glue(\"{Dose_mg} mg\")\n  )\n\n# Grade 2+ hypoglycemia\ndf_er_ae_hgly2 <- d_sim_binom_cov_2 |> filter(AETYPE == \"hgly2\")\n\nvar_resp <- \"AEFLAG\"\nvar_exposure <- \"AUCss_1000\"\n```\n:::\n\n\n\n\n## Fit model\n\n::: {.panel-tabset}\n\n### Linear logistic regression\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_hgly2 <- dev_ermod_bin(\n  data = df_er_ae_hgly2,\n  var_resp = var_resp,\n  var_exposure = var_exposure\n)\nermod_bin_hgly2\n```\n\n<pre class=\"r-output\"><code>\n<span style='color: #00BBBB;'>──</span> <span style='font-weight: bold;'>Binary ER model</span> <span style='color: #00BBBB;'>─────────────────────────────────────────────────────────────</span>\n<span style='color: #00BBBB;'>ℹ</span> Use `plot_er()` to visualize ER curve\n\n── <span style='font-weight: bold;'>Developed model</span> ──\n\nstan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000\n observations: 500\n predictors:   2\n------\n            Median MAD_SD\n(Intercept) -2.04   0.23 \nAUCss_1000   0.41   0.08 \n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n</code></pre>\n\n```{.r .cell-code}\nplot_er_gof(ermod_bin_hgly2, var_group = \"Dose\")\n```\n\n::: {.cell-output-display}\n![](mod_structure_comparison_files/figure-html/fig-bin-lin-1.png){#fig-bin-lin width=672}\n:::\n:::\n\n\n\n\n### E~max~ logistic regression\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_emax_hgly2 <- dev_ermod_bin_emax(\n  data = df_er_ae_hgly2,\n  var_resp = var_resp,\n  var_exposure = var_exposure\n)\nermod_bin_emax_hgly2\n```\n\n<pre class=\"r-output\"><code>\n<span style='color: #00BBBB;'>──</span> <span style='font-weight: bold;'>Binary Emax model</span> <span style='color: #00BBBB;'>───────────────────────────────────────────────────────────</span>\n<span style='color: #00BBBB;'>ℹ</span> Use `plot_er()` to visualize ER curve\n\n── <span style='font-weight: bold;'>Developed model</span> ──\n\n---- Binary Emax model fit with rstanemax ----\n       mean se_mean   sd  2.5%   25%   50%   75% 97.5%   n_eff Rhat\nemax   4.05    0.02 0.91  2.36  3.41  4.02  4.67  5.88 1563.98    1\ne0    -2.39    0.01 0.44 -3.40 -2.62 -2.33 -2.09 -1.68 1033.08    1\nec50   4.77    0.06 2.17  1.44  3.19  4.47  6.06  9.70 1325.28    1\ngamma  1.00     NaN 0.00  1.00  1.00  1.00  1.00  1.00     NaN  NaN\n* Use `extract_stanfit()` function to extract raw stanfit object\n* Use `extract_param()` function to extract posterior draws of key parameters\n* Use `plot()` function to visualize model fit\n* Use `extract_obs_mod_frame()` function to extract raw data \n  in a processed format (useful for plotting)\n</code></pre>\n\n```{.r .cell-code}\nplot_er_gof(ermod_bin_emax_hgly2, var_group = \"Dose\")\n```\n\n::: {.cell-output-display}\n![](mod_structure_comparison_files/figure-html/fig-bin-emax-1.png){#fig-bin-emax width=672}\n:::\n:::\n\n\n\n\n:::\n\n## Model comparison\n\nYou can perform model comparison based on expected log pointwise predictive density (ELPD).\nELPD is the Bayesian leave-one-out estimate (see ?`loo-glossary`).\n\nHigher ELPD is better, therefore linear logistic regression model appears better than E~max~ model.\nHowever, `elpd_diff` is small and similar to `se_diff` ([see here](https://stats.stackexchange.com/questions/608881/how-to-interpret-elpd-diff-of-bayesian-loo-estimate-in-bayesian-logistic-regress)),\ntherefore we can consider the difference to be not meaningful.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nloo_bin_emax_hgly2 <- loo(ermod_bin_emax_hgly2)\nloo_bin_hgly2 <- loo(ermod_bin_hgly2)\n\nloo_compare(list(bin_emax_hgly2 = loo_bin_emax_hgly2, bin_hgly2 = loo_bin_hgly2))\n```\n\n<pre class=\"r-output\"><code>               elpd_diff se_diff\nbin_hgly2       0.0       0.0   \nbin_emax_hgly2 -1.7       1.7   \n</code></pre>\n:::\n",
+    "supporting": [
+      "mod_structure_comparison_files"
+    ],
+    "filters": [
+      "rmarkdown/pagebreak.lua"
+    ],
+    "includes": {},
+    "engineDependencies": {},
+    "preserve": {},
+    "postProcess": true
+  }
+}

_freeze/notebook/emax/data_generation/execute-results/html.json

@@ -0,0 +1,17 @@
+{
+  "hash": "dd04ac6c0ba111a23cafe61fe93ebbcc",
+  "result": {
+    "engine": "knitr",
+    "markdown": "# Data generation\n\nThis page shows examples of data generation for Emax model with and without covariates.\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(here)\n\ntheme_set(theme_bw(base_size = 12))\n\nset.seed(1234)\n```\n:::\n\n\n\n\n## Data generation - No covariate\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nn <- 20 # number of subjects\nE0 <- 5 # effect at 0 concentration\nEmax <- 10 # maximal effect\nEC50 <- 20 # concentration at half maximal effect\nh <- 2 # Hill coefficient\n\nset.seed(130)\nc.is <- 50 * runif(n) # exposure\n\nset.seed(130)\neps <- rnorm(n) # residual error\n\ny.is <- E0 + ((Emax * c.is^h) / (c.is^h + EC50^h)) + eps\n\nd_example_emax_nocov <- tibble(Conc = c.is, Y = y.is)\n```\n:::\n\n\n\n\n### Check data\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(d_example_emax_nocov, aes(Conc, Y)) +\n  geom_point() +\n  geom_smooth(formula = y ~ x, method = \"loess\", se = F, col = \"darkgrey\") +\n  scale_x_continuous(\"Exposure\", breaks = c(3, 10, 30, 100))\n\nggplot(d_example_emax_nocov, aes(Conc, Y)) +\n  geom_point() +\n  geom_smooth(formula = y ~ x, method = \"loess\", se = F, col = \"darkgrey\") +\n  scale_x_log10(\"Exposure\", breaks = c(3, 10, 30, 100))\n```\n\n::: {.cell-output-display}\n![](data_generation_files/figure-html/fig-check-data-no-cov-1.png){#fig-check-data-no-cov-1 width=672}\n:::\n\n::: {.cell-output-display}\n![](data_generation_files/figure-html/fig-check-data-no-cov-2.png){#fig-check-data-no-cov-2 width=672}\n:::\n:::\n\n\n\n\n## Data generation - with covariate\n\nOnly one covariate (Prognostic factor positive/negative)\n\n1. Prognostic factor = positive (GpA) is more sensitive to the drug\n   - lower Emax in GpA; Emax.GpA = 10; Emax.GpB = 15\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nNgp <- 2\nN <- 20 * Ngp\nGPid <- as.factor(rep(c(\"A\", \"B\"), each = 20))\n\n# Set parameters\nE0 <- 5\nEC50 <- 15\nh <- 2\nbeta1 <- .7\n\n# Calc response\nset.seed(12345)\n\nd_example_emax_cov <-\n  tibble(GP = GPid) |>\n  mutate(\n    c.is = 50 * runif(N), eps = rnorm(N)\n  ) |>\n  mutate(\n    Emax.i = ifelse(GP == \"A\", 10, 15)\n  ) |>\n  mutate(y.is = E0 + ((Emax.i * c.is^h) / (c.is^h + EC50^h)) + eps) |>\n  mutate(Conc = c.is, Y = y.is)\n```\n:::\n\n\n\n\n### Check data\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(d_example_emax_cov, aes(Conc, Y)) +\n  geom_point(aes(colour = GP)) +\n  geom_smooth(aes(group = GP, colour = GP), se = F) +\n  scale_x_continuous(\"Exposure\") +\n  theme(legend.position = \"top\")\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\n`geom_smooth()` using method = 'loess' and formula = 'y ~ x'\n```\n\n\n:::\n\n::: {.cell-output-display}\n![](data_generation_files/figure-html/fig-check-data-cov-1.png){#fig-check-data-cov width=672}\n:::\n:::\n\n\n\n\n## Save data\n\nOnly run in an interactive session so that the data is not saved every time\nthe document is rendered (by setting `eval: FALSE`).\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nwrite_csv(d_example_emax_nocov, here(\"data\", \"d_example_emax_nocov.csv\"))\nwrite_csv(d_example_emax_cov, here(\"data\", \"d_example_emax_cov.csv\"))\n```\n:::\n",
+    "supporting": [
+      "data_generation_files"
+    ],
+    "filters": [
+      "rmarkdown/pagebreak.lua"
+    ],
+    "includes": {},
+    "engineDependencies": {},
+    "preserve": {},
+    "postProcess": true
+  }
+}

_freeze/notebook/emax/simulation/execute-results/html.json

@@ -0,0 +1,17 @@
+{
+  "hash": "f1ca9b0b653d905dab2dbe0756ef4ec3",
+  "result": {
+    "engine": "knitr",
+    "markdown": "# Simulation from fitted model\n\nThis page showcase the model simulation using the E~max~ model with no covariate.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(posterior)\nlibrary(tidybayes)\nlibrary(bayesplot)\nlibrary(loo)\nlibrary(here)\nlibrary(gt)\n\ntheme_set(theme_bw(base_size = 12))\n```\n:::\n\n\n\n\n\n\n## Load data\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_example_emax_nocov <- read_csv(here(\"data\", \"d_example_emax_nocov.csv\"))\n```\n:::\n\n\n\n\n## Fit models\n\n::: {.panel-tabset}\n\n### Sigmoidal E~max~ model\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\nermod_sigemax <- dev_ermod_emax(\n  data = d_example_emax_nocov,\n  var_resp = \"Y\",\n  var_exposure = \"Conc\",\n  gamma_fix = NULL\n)\n\nermod_sigemax\n```\n\n<pre class=\"r-output\"><code>\n<span style='color: #00BBBB;'>──</span> <span style='font-weight: bold;'>Emax model</span> <span style='color: #00BBBB;'>──────────────────────────────────────────────────────────────────</span>\n<span style='color: #00BBBB;'>ℹ</span> Use `plot_er()` to visualize ER curve\n\n── <span style='font-weight: bold;'>Developed model</span> ──\n\n---- Emax model fit with rstanemax ----\n       mean se_mean   sd  2.5%   25%   50%   75% 97.5%   n_eff Rhat\nemax  15.92    0.24 5.49  8.67 11.64 14.71 19.07 29.38  545.77 1.01\ne0     3.13    0.09 2.00 -2.15  2.25  3.58  4.54  5.62  550.23 1.00\nec50  24.94    0.30 8.63 14.98 19.47 22.51 28.06 48.31  843.12 1.00\ngamma  1.77    0.03 0.77  0.75  1.22  1.62  2.19  3.54  664.61 1.00\nsigma  0.86    0.00 0.16  0.62  0.75  0.84  0.95  1.24 1348.22 1.00\n* Use `extract_stanfit()` function to extract raw stanfit object\n* Use `extract_param()` function to extract posterior draws of key parameters\n* Use `plot()` function to visualize model fit\n* Use `posterior_predict()` or `posterior_predict_quantile()` function to get\n  raw predictions or make predictions on new data\n* Use `extract_obs_mod_frame()` function to extract raw data \n  in a processed format (useful for plotting)\n</code></pre>\n:::\n\n\n\n\n### Regular E~max~ (h fixed at 1)\n\nAnother model without sigmoidal component; will be used when we do\nmodel comparison.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\nermod_emax <- dev_ermod_emax(\n  data = d_example_emax_nocov,\n  var_resp = \"Y\",\n  var_exposure = \"Conc\",\n  gamma_fix = 1\n)\n\nermod_emax\n```\n\n<pre class=\"r-output\"><code>\n<span style='color: #00BBBB;'>──</span> <span style='font-weight: bold;'>Emax model</span> <span style='color: #00BBBB;'>──────────────────────────────────────────────────────────────────</span>\n<span style='color: #00BBBB;'>ℹ</span> Use `plot_er()` to visualize ER curve\n\n── <span style='font-weight: bold;'>Developed model</span> ──\n\n---- Emax model fit with rstanemax ----\n       mean se_mean    sd  2.5%   25%   50%   75% 97.5%   n_eff Rhat\nemax  22.52    0.08  2.46 18.73 20.84 22.17 23.81 28.09  976.89 1.00\ne0     0.98    0.04  1.36 -2.23  0.32  1.18  1.91  2.98  912.27 1.00\nec50  31.28    0.39 11.65 13.08 22.74 29.70 38.10 57.39  875.72 1.00\ngamma  1.00     NaN  0.00  1.00  1.00  1.00  1.00  1.00     NaN  NaN\nsigma  0.89    0.00  0.17  0.63  0.77  0.86  0.98  1.28 1216.16 1.01\n* Use `extract_stanfit()` function to extract raw stanfit object\n* Use `extract_param()` function to extract posterior draws of key parameters\n* Use `plot()` function to visualize model fit\n* Use `posterior_predict()` or `posterior_predict_quantile()` function to get\n  raw predictions or make predictions on new data\n* Use `extract_obs_mod_frame()` function to extract raw data \n  in a processed format (useful for plotting)\n</code></pre>\n:::\n\n\n\n\n:::\n\n# Extrapolation \n\n::: {.panel-tabset}\n\n## without residual error\n\nThis represents uncertainty in the model parameters.\n\n- `p(f(theta)|xnew, yobs)`\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nnew_conc_vec <- seq(0, 100, by = 1)\n\nersim_sigemax <-\n  sim_er_new_exp(ermod_sigemax, new_conc_vec)\n\nersim_sigemax_med_qi <-\n  ersim_sigemax |>\n  calc_ersim_med_qi(qi_width = c(0.5, 0.95))\n\nggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .epred)) +\n  geom_ribbon(\n    data = ersim_sigemax_med_qi |> filter(.width == 0.95),\n    aes(ymin = .epred.lower, ymax = .epred.upper),\n    fill = \"yellow3\", alpha = 0.5\n  ) +\n  geom_ribbon(\n    data = ersim_sigemax_med_qi |> filter(.width == 0.5),\n    aes(ymin = .epred.lower, ymax = .epred.upper),\n    fill = \"orange1\"\n  ) +\n  geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.5), col = \"darkred\") +\n  geom_point(data = d_example_emax_nocov, aes(y = Y)) +\n  coord_cartesian(ylim = c(-1, 22)) +\n  scale_y_continuous(\"Response\", breaks = 5 * 0:4) +\n  labs(x = \"Exposure\",\n    title = \"Extrapolation (no residual)\",\n    subtitle = \"Represents uncertainty in the model parameters (Credible interval)\",\n    caption = \"95% CI in yellow, 50% CI in orange\")\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-extrapolation-1-1.png){#fig-extrapolation-1 width=672}\n:::\n:::\n\n\n\n\n## with residual error\n\nThis represents uncertainty in the model parameters plus the residual error.\n\n- `p(ynew|xnew, yobs)`\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .prediction)) +\n  geom_ribbon(\n    data = ersim_sigemax_med_qi |> filter(.width == 0.95),\n    aes(ymin = .prediction.lower, ymax = .prediction.upper),\n    fill = \"yellow3\", alpha = 0.5\n  ) +\n  geom_ribbon(\n    data = ersim_sigemax_med_qi |> filter(.width == 0.5),\n    aes(ymin = .prediction.lower, ymax = .prediction.upper),\n    fill = \"orange1\"\n  ) +\n  geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.5), col = \"darkred\") +\n  geom_point(data = d_example_emax_nocov, aes(y = Y)) +\n  coord_cartesian(ylim = c(-1, 22)) +\n  scale_y_continuous(\"Response\", breaks = 5 * 0:4) +\n  labs(x = \"Exposure\",\n    title = \"Extrapolation (incl. residual)\",\n    subtitle = \"Represents uncertainty + residual error (Prediction interval)\",\n    caption = \"95% PI in yellow, 50% PI in orange\")\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-extrapolation-2-1.png){#fig-extrapolation-2 width=672}\n:::\n:::\n\n\n\n\n## Overlay Emax and Sigmoidal Emax\n\nNo discernible difference between the two models.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_emax <-\n  sim_er_new_exp(ermod_emax, new_conc_vec)\nersim_emax_med_qi <-\n  ersim_emax |>\n  calc_ersim_med_qi(qi_width = c(0.5, 0.95))\n\nggplot(data = ersim_sigemax_med_qi, aes(x = Conc, y = .epred)) +\n  geom_ribbon(\n    data = ersim_sigemax_med_qi |> filter(.width == 0.95),\n    aes(ymin = .epred.lower, ymax = .epred.upper),\n    fill = \"orange1\", alpha = 0.5\n  ) +\n  geom_line(data = ersim_sigemax_med_qi |> filter(.width == 0.95), col = \"darkred\") +\n  geom_ribbon(\n    data = ersim_emax_med_qi |> filter(.width == 0.95),\n    aes(ymin = .epred.lower, ymax = .epred.upper),\n    fill = \"turquoise3\", alpha = 0.5\n  ) +\n  geom_line(data = ersim_emax_med_qi |> filter(.width == 0.95), col = \"steelblue3\") +\n  geom_point(data = d_example_emax_nocov, aes(y = Y)) +\n  coord_cartesian(ylim = c(-1, 22)) +\n  scale_y_continuous(\"Response\", breaks = 5 * 0:4) +\n  labs(x = \"Exposure\",\n    title = \"Extrapolation (no residual)\",\n    subtitle = \"Represents uncertainty in the model parameters (Credible interval)\",\n    caption = \"Sigmoidal Emax: Orange, Emax: Blue, 50% CI\")\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-emax-sigemax-compare-1.png){#fig-emax-sigemax-compare width=672}\n:::\n:::\n\n\n\n\n:::\n\n\n",
+    "supporting": [
+      "simulation_files"
+    ],
+    "filters": [
+      "rmarkdown/pagebreak.lua"
+    ],
+    "includes": {},
+    "engineDependencies": {},
+    "preserve": {},
+    "postProcess": true
+  }
+}