diff --git a/.github/workflows/publish.yml b/.github/workflows/publish.yml new file mode 100644 index 0000000..e82eb7b --- /dev/null +++ b/.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 }} diff --git a/README.Rmd b/README.Rmd index 1961488..1408537 100644 --- a/README.Rmd +++ b/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 diff --git a/README.md b/README.md index 771f6c1..721a45f 100644 --- a/README.md +++ b/README.md @@ -95,7 +95,7 @@ install.packages('rstanemax') ## Important links -- Rendered book: +- Rendered book: ## Note for developer diff --git a/_freeze/index/execute-results/html.json b/_freeze/index/execute-results/html.json new file mode 100644 index 0000000..0e21801 --- /dev/null +++ b/_freeze/index/execute-results/html.json @@ -0,0 +1,15 @@ +{ + "hash": "ded2f03a0ffb872db228ca1ac3c691d2", + "result": { + "engine": "knitr", + "markdown": "---\nformat: html\n---\n\n\n\n\n# Welcome {.unnumbered}\n\n\n\n\n\nThis book provides examples of exposure-response analysis with Bayesian\nmethods.\n\n## Install necessary packages \"BayesERtools\n\nThe examples utilizes\n[`BayesERtools`](https://genentech.github.io/BayesERtools/) package.\n\n- Tutorial (This page): \n- Documentation: \n- GitHub repo: \n\nYou can install the package as follows:\n\n``` r\n# install.packages('BayesERtools') # Once on CRAN\ndevtools::install_github(\"genentech/BayesERtools\") # development version\n```\n\nYou also need latest version of `rstanemax` (\\>= 0.1.7) to use\nEmax model.\n\n``` r\ninstall.packages('rstanemax')\n```\n\n## Model types supported by `BayesERtools`\n\n
\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n \n\n\n\n\n \n\n\n\n\n \n\n\n\n\n \n\n\n\n\n \n\n\n\n\n \n\n\n\n\n \n \n \n \n \n \n
\n
Binary endpoint
\n 		\n
Continuous endpoint
\n
Linear (logit)Emax (logit)LinearEmax
backendrstanarmrstanemaxrstanarmrstanemax
referenceπŸ”—πŸ”—πŸ”—πŸ”—
develop modelβœ…βœ…βœ…βœ…
simulate & plot ERβœ…βœ…βœ…βœ…
exposure metrics selectionβœ…βœ…βœ…βœ…
covariate selectionβœ…βŒβœ…βŒ
covariate forest plotβœ…βŒπŸŸ‘βŒ
βœ… Available, 🟑 In plan/under development, ❌ Not in a current plan
\n
\n\n## Important links\n\n- Rendered book: \n\n## Note for developer\n\nRun `usethis::use_tidy_style(strict = FALSE)` before committing to\nensure that the code is formatted appropriately.\n\n\n\n\n## Session info\n\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndevtools::session_info()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n─ Session info ───────────────────────────────────────────────────────────────\n setting value\n version R version 4.4.1 (2024-06-14)\n os Ubuntu 22.04.5 LTS\n system x86_64, linux-gnu\n ui X11\n language (EN)\n collate en_US.UTF-8\n ctype en_US.UTF-8\n tz Etc/UTC\n date 2025-02-07\n pandoc 3.1.11 @ /usr/lib/rstudio-server/bin/quarto/bin/tools/x86_64/ (via rmarkdown)\n\n─ Packages ───────────────────────────────────────────────────────────────────\n package * version date (UTC) lib source\n abind 1.4-8 2024-09-12 [2] RSPM (R 4.4.0)\n arrayhelpers 1.1-0 2020-02-04 [1] RSPM (R 4.4.0)\n backports 1.5.0 2024-05-23 [2] RSPM (R 4.4.0)\n BayesERtools * 0.1.11 2025-01-30 [1] local\n bayesplot * 1.11.1 2024-02-15 [2] RSPM (R 4.4.0)\n bayestestR 0.14.0 2024-07-24 [1] RSPM (R 4.4.0)\n cachem 1.1.0 2024-05-16 [2] RSPM (R 4.4.0)\n checkmate 2.3.2 2024-07-29 [2] RSPM (R 4.4.0)\n cli 3.6.3 2024-06-21 [2] RSPM (R 4.4.0)\n coda 0.19-4.1 2024-01-31 [2] RSPM (R 4.4.0)\n colorspace 2.1-1 2024-07-26 [2] RSPM (R 4.4.0)\n devtools 2.4.5 2022-10-11 [1] RSPM (R 4.4.0)\n digest 0.6.37 2024-08-19 [2] RSPM (R 4.4.0)\n distributional 0.5.0 2024-09-17 [2] RSPM (R 4.4.0)\n dplyr * 1.1.4 2023-11-17 [2] RSPM (R 4.4.0)\n ellipsis 0.3.2 2021-04-29 [1] RSPM (R 4.4.0)\n evaluate 1.0.1 2024-10-10 [2] RSPM (R 4.4.0)\n fansi 1.0.6 2023-12-08 [2] RSPM (R 4.4.0)\n fastmap 1.2.0 2024-05-15 [2] RSPM (R 4.4.0)\n forcats * 1.0.0 2023-01-29 [2] RSPM (R 4.4.0)\n fs 1.6.4 2024-04-25 [2] RSPM (R 4.4.0)\n generics 0.1.3 2022-07-05 [2] RSPM (R 4.4.0)\n ggdist 3.3.2 2024-03-05 [1] RSPM (R 4.4.0)\n ggplot2 * 3.5.1 2024-04-23 [2] RSPM (R 4.4.0)\n glue 1.8.0 2024-09-30 [2] RSPM (R 4.4.0)\n gt * 0.11.1 2024-10-04 [1] RSPM (R 4.4.0)\n gtable 0.3.5 2024-04-22 [2] RSPM (R 4.4.0)\n here * 1.0.1 2020-12-13 [1] RSPM (R 4.4.0)\n hms 1.1.3 2023-03-21 [2] RSPM (R 4.4.0)\n htmltools 0.5.8.1 2024-04-04 [2] RSPM (R 4.4.0)\n htmlwidgets 1.6.4 2023-12-06 [2] RSPM (R 4.4.0)\n httpuv 1.6.15 2024-03-26 [2] RSPM (R 4.4.0)\n insight 0.20.5 2024-10-02 [1] RSPM (R 4.4.0)\n jsonlite 1.8.9 2024-09-20 [2] RSPM (R 4.4.0)\n knitr 1.48 2024-07-07 [2] RSPM (R 4.4.0)\n later 1.3.2 2023-12-06 [2] RSPM (R 4.4.0)\n lattice 0.22-6 2024-03-20 [4] RSPM (R 4.4.1)\n lifecycle 1.0.4 2023-11-07 [2] RSPM (R 4.4.0)\n loo * 2.8.0 2024-07-03 [2] RSPM (R 4.4.0)\n lubridate * 1.9.3 2023-09-27 [2] RSPM (R 4.4.0)\n magrittr 2.0.3 2022-03-30 [2] RSPM (R 4.4.0)\n matrixStats 1.4.1 2024-09-08 [2] RSPM (R 4.4.0)\n memoise 2.0.1 2021-11-26 [2] RSPM (R 4.4.0)\n mime 0.12 2021-09-28 [2] RSPM (R 4.4.0)\n miniUI 0.1.1.1 2018-05-18 [2] RSPM (R 4.4.0)\n munsell 0.5.1 2024-04-01 [2] RSPM (R 4.4.0)\n pillar 1.9.0 2023-03-22 [2] RSPM (R 4.4.0)\n pkgbuild 1.4.4 2024-03-17 [2] RSPM (R 4.4.0)\n pkgconfig 2.0.3 2019-09-22 [2] RSPM (R 4.4.0)\n pkgload 1.4.0 2024-06-28 [1] RSPM (R 4.4.0)\n posterior * 1.6.0 2024-07-03 [2] RSPM (R 4.4.0)\n profvis 0.4.0 2024-09-20 [1] RSPM (R 4.4.0)\n promises 1.3.0 2024-04-05 [2] RSPM (R 4.4.0)\n purrr * 1.0.2 2023-08-10 [2] RSPM (R 4.4.0)\n R6 2.5.1 2021-08-19 [2] RSPM (R 4.4.0)\n Rcpp 1.0.13 2024-07-17 [2] RSPM (R 4.4.0)\n readr * 2.1.5 2024-01-10 [2] RSPM (R 4.4.0)\n remotes 2.5.0 2024-03-17 [2] RSPM (R 4.4.0)\n rlang 1.1.4 2024-06-04 [2] RSPM (R 4.4.0)\n rmarkdown 2.28 2024-08-17 [2] RSPM (R 4.4.0)\n rprojroot 2.0.4 2023-11-05 [1] RSPM (R 4.4.0)\n rstudioapi 0.16.0 2024-03-24 [2] RSPM (R 4.4.0)\n scales 1.3.0 2023-11-28 [2] RSPM (R 4.4.0)\n sessioninfo 1.2.2 2021-12-06 [1] RSPM (R 4.4.0)\n shiny 1.9.1 2024-08-01 [2] RSPM (R 4.4.0)\n stringi 1.8.4 2024-05-06 [2] RSPM (R 4.4.0)\n stringr * 1.5.1 2023-11-14 [2] RSPM (R 4.4.0)\n svUnit 1.0.6 2021-04-19 [1] RSPM (R 4.4.0)\n tensorA 0.36.2.1 2023-12-13 [2] RSPM (R 4.4.0)\n tibble * 3.2.1 2023-03-20 [2] RSPM (R 4.4.0)\n tidybayes * 3.0.7 2024-09-15 [1] RSPM (R 4.4.0)\n tidyr * 1.3.1 2024-01-24 [2] RSPM (R 4.4.0)\n tidyselect 1.2.1 2024-03-11 [2] RSPM (R 4.4.0)\n tidyverse * 2.0.0 2023-02-22 [2] RSPM (R 4.4.0)\n timechange 0.3.0 2024-01-18 [2] RSPM (R 4.4.0)\n tzdb 0.4.0 2023-05-12 [2] RSPM (R 4.4.0)\n urlchecker 1.0.1 2021-11-30 [1] RSPM (R 4.4.0)\n usethis 3.0.0 2024-07-29 [1] RSPM (R 4.4.0)\n utf8 1.2.4 2023-10-22 [2] RSPM (R 4.4.0)\n vctrs 0.6.5 2023-12-01 [2] RSPM (R 4.4.0)\n withr 3.0.1 2024-07-31 [2] RSPM (R 4.4.0)\n xfun 0.48 2024-10-03 [4] RSPM (R 4.4.1)\n xml2 1.3.6 2023-12-04 [2] RSPM (R 4.4.0)\n xtable 1.8-4 2019-04-21 [2] RSPM (R 4.4.0)\n yaml 2.3.10 2024-07-26 [2] RSPM (R 4.4.0)\n\n [1] /home/yoshidk6/R/x86_64-pc-linux-gnu-library/4.4\n [2] /usr/local/lib/R/site-library\n [3] /usr/lib/R/site-library\n [4] /usr/lib/R/library\n\n──────────────────────────────────────────────────────────────────────────────\n```\n\n\n:::\n:::\n", + "supporting": [], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/basic_workflow/execute-results/html.json b/_freeze/notebook/binary/basic_workflow/execute-results/html.json new file mode 100644 index 0000000..56c4097 --- /dev/null +++ b/_freeze/notebook/binary/basic_workflow/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "c22b5c3299b6c5e4c5f46b13efbcaab5", + "result": { + "engine": "knitr", + "markdown": "# Basic workflow\n\nIn this section, we will show a basic workflow of performing ER analysis\nfor binary endpoint using the logistic regression model.\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(posterior)\nlibrary(here)\nlibrary(gt)\n\ntheme_set(theme_bw(base_size = 12))\n```\n:::\n\n\n\n\n\n\n\n## Data\n\nWe will use an example simulated dataset included in `BayesERtools` package \n(`d_sim_binom_cov`) for the analysis.\nIn this document we use hypoglycemia (hgly2) as an example AE. Another example\nAE is diarrhea (dr2), where you would see fairly flat ER curve.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndata(d_sim_binom_cov)\n\nhead(d_sim_binom_cov) |> \n gt() |>\n fmt_number(n_sigfig = 3) |>\n fmt_integer(columns = c(\"ID\", \"AEFLAG\", \"Dose_mg\"))\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# Grade 2+ diarrhea\ndf_er_ae_dr2 <- d_sim_binom_cov_2 |> filter(AETYPE == \"dr2\")\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n \n \n \n
IDAETYPEAEFLAGDose_mgAUCssCmaxssCminssBAGEBWTBGLUCBHBA1CRACEVISC
1hgly2020086664.310.184.474.14.6531.5WhiteNo
1dr2020086664.310.184.474.14.6531.5WhiteNo
1ae_covsel_test020086664.310.184.474.14.6531.5WhiteNo
2hgly202001,71016627.359.188.27.2441.9WhiteNo
2dr202001,71016627.359.188.27.2441.9WhiteNo
2ae_covsel_test12001,71016627.359.188.27.2441.9WhiteNo
\n
\n```\n\n:::\n:::\n\n\n\n\nWe also defines variables that is used in the analysis.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nvar_resp <- \"AEFLAG\"\n# HbA1c & glucose are only relevant for hyperglycemia\nvar_cov_ae_hgly2 <-\n c(\"BAGE_10\", \"BWT_10\", \"RACE\", \"VISC\", \"BHBA1C_5\", \"BGLUC\")\nvar_cov_ae_dr2 <-\n c(\"BAGE_10\", \"BWT_10\", \"RACE\", \"VISC\")\n```\n:::\n\n\n\n\n\n## Basic model development\n\n`dev_ermod_bin()` function can be used to develop basic ER model. (Note that\nthis function can also be used for models with covariates, if you already know\nthe covariates to be included.)\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin <- dev_ermod_bin(\n data = df_er_ae_hgly2,\n var_resp = var_resp,\n var_exposure = \"AUCss_1000\"\n)\nermod_bin\n```\n\n
\n── Binary ER model ─────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\n:::\n\n\n\n\nYou can compare the observed data with the model fit using `plot_er_gof()` \nfunction.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Using `*` instead of `+` so that scale can be\n# applied for both panels (main plot and boxplot)\nplot_er_gof(ermod_bin, var_group = \"Dose\", show_coef_exp = TRUE) *\n xgxr::xgx_scale_x_log10()\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-plot-er-1.png){#fig-plot-er width=672}\n:::\n:::\n\n\n\n\nMCMC samples can be obtained with `as_draws()` family of functions, such as\n`as_draws_df()`.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndraws_df <- as_draws_df(ermod_bin)\n\ndraws_df_summary <-\n draws_df |> \n summarize_draws(mean, median, sd, ~quantile2(.x, probs = c(0.025, 0.975)), \n default_convergence_measures()) \n\ndraws_df_summary |>\n gt() |>\n fmt_number(n_sigfig = 3)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdq2.5q97.5rhatess_bulkess_tail
(Intercept)βˆ’2.05βˆ’2.040.234βˆ’2.51βˆ’1.591.002,2301,870
AUCss_10000.4120.4120.07610.2650.5611.002,1502,140
\n
\n```\n\n:::\n:::\n\n\n\n\n\nYou can predict the probability of events for a given exposure level with\n`sim_er_new_exp()` function.\n\nHere, the prediction is done for AUCss_1000 of 1, 1.5, 2, 3 (AUCss of 1000, \n1500, 2000, 3000), and the output is the median and 95% CI of the predicted \nprobability.\nYou can set `output_type = \"draw\"` to get the raw posterior draws.\n\nThere are two types of outputs here, `.epred` and `.linpred`, as follows:\n\n- `.epred`: Expected response on the probability scale (% of event)\n- `.linpred`: Expected response on the linear predictor scale (logit scale, \n i.e. log-odds)\n \nSee `?BayesERtools::sim_er` for more details.\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_med_qi <- sim_er_new_exp(\n ermod_bin,\n exposure_to_sim_vec = c(1, 1.5, 2, 3),\n output_type = \"median_qi\"\n)\nersim_med_qi |> \n gt() |>\n fmt_number(n_sigfig = 3) |>\n fmt_integer(columns = c(\".row\"))\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n \n \n
AUCss_1000.row.epred.epred.lower.epred.upper.linpred.linpred.lower.linpred.upper.width.point.interval
1.0010.1640.1220.213βˆ’1.63βˆ’1.97βˆ’1.310.950medianqi
1.5020.1940.1530.240βˆ’1.43βˆ’1.71βˆ’1.150.950medianqi
2.0030.2280.1890.270βˆ’1.22βˆ’1.46βˆ’0.9930.950medianqi
3.0040.3080.2660.354βˆ’0.811βˆ’1.02βˆ’0.6030.950medianqi
\n
\n```\n\n:::\n:::\n\n\n\n\n## Selection of exposure metrics\n\n`dev_ermod_bin_exp_sel()` function can be used to select the best exposure\nmetric(s) from a list of candidate exposure metrics. In this case, AUCss_1000\nis selected as the best exposure metric, as it had the highest elpd (expected \nlog predictive density).[^1]\n\n[^1]: Some references about elpd: \\code{?loo::`loo-glossary`}, [What is the interpretation of ELPD / elpd_loo / elpd_diff?](https://mc-stan.org/loo/articles/online-only/faq.html#elpd_interpretation)\n\nNote that whether you want to select exposure metrics using the statistical\ncriteria (e.g. elpd) or pre-specify the exposure metric(s) depends on the\ncontexts. Should you choose to pre-specify the exposure metric(s), you can\nskip this step.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_exp_sel <-\n dev_ermod_bin_exp_sel(\n # Use reduced N to make the example run faster\n data = slice_sample(df_er_ae_hgly2, n = 100),\n var_resp = var_resp,\n var_exp_candidates = c(\"AUCss_1000\", \"Cmaxss\", \"Cminss\"),\n # Use reduced iter to make the example run faster\n iter = 1000\n )\n```\n\n
β„Ή The exposure metric selected was: AUCss_1000\n
\n\n```{.r .cell-code}\nermod_bin_exp_sel\n```\n\n
\n── Binary ER model & exposure metric selection ─────────────────────────────────\nβ„Ή Use `plot_er_exp_sel()` for ER curve of all exposure metrics\nβ„Ή Use `plot_er()` with `show_orig_data = TRUE` for ER curve of the selected exposure metric\n\n── Exposure metrics comparison ──\n\n           elpd_diff se_diff\nAUCss_1000  0.00      0.00  \nCminss     -0.76      0.70  \nCmaxss     -0.91      0.97  \n\n── Selected model ──\n\nstan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000\n observations: 100\n predictors:   2\n------\n            Median MAD_SD\n(Intercept) -1.59   0.45 \nAUCss_1000   0.20   0.15 \n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n
\n:::\n\n\n\n\nThe ER curve for all the evaluated exposure metrics can be generated with \n`plot_er_exp_sel()` function.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot_er_exp_sel(ermod_bin_exp_sel) +\n xgxr::xgx_scale_x_log10()\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-exp-sel-1.png){#fig-exp-sel width=672}\n:::\n:::\n\n\n\n\n\n## Selection of covariates\n\n`dev_ermod_bin_cov_sel()` function can be used to select the best covariates\nfrom a list of candidate covariates. In this case, HbA1c (BHBA1C_5) and glucose\n(BGLUC) are selected, in addition to the exposure metric AUCss_1000 as \npredictors.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_cov_sel <-\n dev_ermod_bin_cov_sel(\n data = df_er_ae_hgly2,\n var_resp = var_resp,\n var_exposure = \"AUCss_1000\",\n var_cov_candidate = var_cov_ae_hgly2,\n verbosity_level = 2\n )\n```\n:::\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nermod_bin_cov_sel\n```\n\n
\n── Binary ER model & covariate selection ───────────────────────────────────────\nβ„Ή Use `plot_submod_performance()` to see variable selection performance\nβ„Ή Use `plot_er()` with `marginal = TRUE` to visualize marginal ER curve\n\n── Selected model ──\n\nstan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000 + BHBA1C_5 + BGLUC\n observations: 500\n predictors:   4\n------\n            Median MAD_SD\n(Intercept) -11.00   1.12\nAUCss_1000    0.49   0.08\nBHBA1C_5      0.50   0.09\nBGLUC         0.74   0.13\n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n
\n:::\n\n\n\n\nThe plot below shows that AUCss_1000, BHBA1C_5, and BGLUC contributes to \nimproving the model performance, and after then the inclusion of no other\ncovariates improves the model performance.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot_submod_performance(ermod_bin_cov_sel)\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-cov-sel-1.png){#fig-cov-sel width=672}\n:::\n:::\n\n\n\n\nIn some cases, you might see a warning message like below. This indicates that\napproximation of leave-one-out cross-validation performance (PSIS-LOO) is not\nreliable.\n\n```\nWarning: In the recalculation of the reference model's PSIS-LOO CV weights for \nthe performance evaluation, ... Pareto k-values are in the interval...`. \n```\n\nAlternatively to the default `cv_method = \"LOO\"`, you can use k-fold \ncross-validation by setting`cv_method = \"kfold\"` in `dev_ermod_bin_cov_sel()` \nfunction. This can take longer time to run, but it can be more reliable in\nthe cases where LOO is not reliable. You can also set `validate_search = TRUE`\nto let the function perform the variable selection for each fold separately,\nrather than using the selected variable sequence from the full dataset\nevaluation. \n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_cov_sel_kfold <-\n dev_ermod_bin_cov_sel(\n data = df_er_ae_hgly2,\n var_resp = var_resp,\n var_exposure = \"AUCss_1000\",\n var_cov_candidate = var_cov_ae_hgly2,\n cv_method = \"kfold\",\n validate_search = TRUE,\n verbosity_level = 2\n )\n```\n:::\n\n\n\n\n\n\n\n\nAdded bonus of using k-fold cv is that you can visualize how often each \nvariable is selected in the model. Here, as you can see (and as expected),\nHbA1c (BHBA1C_5) and glucose (BGLUC) are highly related and they are almost\ninterchangeably selected in the 2nd and 3rd positions. Note that the function\nenforces the exposure metric to be included first in the model.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nermod_bin_cov_sel_kfold\n```\n\n
\n── Binary ER model & covariate selection ───────────────────────────────────────\nβ„Ή Use `plot_submod_performance()` to see variable selection performance\nβ„Ή Use `plot_var_ranking()` to see variable ranking\nβ„Ή Use `plot_er()` with `marginal = TRUE` to visualize marginal ER curve\n\n── Selected model ──\n\nstan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000 + BHBA1C_5 + BGLUC\n observations: 500\n predictors:   4\n------\n            Median MAD_SD\n(Intercept) -10.97   1.20\nAUCss_1000    0.49   0.09\nBHBA1C_5      0.50   0.09\nBGLUC         0.73   0.14\n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n
\n\n```{.r .cell-code}\nplot_submod_performance(ermod_bin_cov_sel_kfold)\nplot_var_ranking(ermod_bin_cov_sel_kfold)\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-cov-sel-kfold-1.png){#fig-cov-sel-kfold-1 width=672}\n:::\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-cov-sel-kfold-2.png){#fig-cov-sel-kfold-2 width=672}\n:::\n:::\n\n\n\n\n## Marginal ER prediction\n\nYou can simulate the marginal ER relationship, i.e. ER relationships for\n\"marginalized\", or averaged, response for the population of interest,\nusing `sim_er_new_exp_marg()` function.\nBy default, the covariate distribution is from the original data, but you can\nalso supply other distribution with `data_cov` argument.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_new_exp_marg_med_qi <- sim_er_new_exp_marg(\n ermod_bin_cov_sel,\n exposure_to_sim_vec = c(2:6),\n output_type = \"median_qi\"\n)\nersim_new_exp_marg_med_qi\n\nplot_er(ersim_new_exp_marg_med_qi, marginal = TRUE)\n```\n\n
# A tibble: 5 Γ— 14\n  .id_exposure AUCss_1000 .epred .epred.lower .epred.upper .linpred\n                                     \n1            1          2  0.226        0.190        0.263  -1.50  \n2            2          3  0.306        0.266        0.348  -1.01  \n3            3          4  0.399        0.338        0.457  -0.518 \n4            4          5  0.496        0.412        0.583  -0.0288\n5            5          6  0.595        0.482        0.701   0.470 \n# β„Ή 8 more variables: .linpred.lower , .linpred.upper ,\n#   .prediction , .prediction.lower , .prediction.upper ,\n#   .width , .point , .interval \n
::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-marginal-er-1.png){#fig-marginal-er width=672}\n:::\n:::\n\n\n\n\n## Evaluation of covariate effects\n\nYou can visualize the effect of the covariates with `sim_coveff()` and\n`plot_coveff()` functions. You can see that all three predictors have \nfairly strong effects on the odds ratio of hypoglycemia.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncoveffsim <- sim_coveff(ermod_bin_cov_sel)\nplot_coveff(coveffsim)\nprint_coveff(coveffsim)\n```\n\n
# A tibble: 9 Γ— 5\n  var_label  value_label value_annot `Odds ratio` `95% CI`        \n                                         \n1 AUCss_1000 0.868       5th         0.514        \"[0.427, 0.613]\"\n2 AUCss_1000 2.21        median      1            \" \"             \n3 AUCss_1000 5.30        95th        4.60         \"[3.08, 7.04]\"  \n4 BHBA1C_5   5.75        5th         0.330        \"[0.237, 0.455]\"\n5 BHBA1C_5   7.97        median      1            \" \"             \n6 BHBA1C_5   10.4        95th        3.44         \"[2.41, 4.97]\"  \n7 BGLUC      4.43        5th         0.293        \"[0.198, 0.420]\"\n8 BGLUC      6.10        median      1            \" \"             \n9 BGLUC      7.59        95th        3.00         \"[2.17, 4.25]\"  \n
::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-cov-eff-1.png){#fig-cov-eff width=576}\n:::\n:::\n", + "supporting": [ + "basic_workflow_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-eff-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-eff-1.png new file mode 100644 index 0000000..79051f9 Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-eff-1.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-1.png new file mode 100644 index 0000000..9c61ba5 Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-1.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-1.png new file mode 100644 index 0000000..71013e3 Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-1.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-2.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-2.png new file mode 100644 index 0000000..53fdd26 Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-cov-sel-kfold-2.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-exp-sel-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-exp-sel-1.png new file mode 100644 index 0000000..a8bb15f Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-exp-sel-1.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-marginal-er-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-marginal-er-1.png new file mode 100644 index 0000000..b6195d4 Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-marginal-er-1.png differ diff --git a/_freeze/notebook/binary/basic_workflow/figure-html/fig-plot-er-1.png b/_freeze/notebook/binary/basic_workflow/figure-html/fig-plot-er-1.png new file mode 100644 index 0000000..300d48c Binary files /dev/null and b/_freeze/notebook/binary/basic_workflow/figure-html/fig-plot-er-1.png differ diff --git a/_freeze/notebook/binary/coveff_customize/execute-results/html.json b/_freeze/notebook/binary/coveff_customize/execute-results/html.json new file mode 100644 index 0000000..70da528 --- /dev/null +++ b/_freeze/notebook/binary/coveff_customize/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "2b65258c5886183cec6404c45d4d7c4b", + "result": { + "engine": "knitr", + "markdown": "# Customize covariate effect plots\n\nThis page showcase how to customize the covariate effect plots.\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\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}\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 )\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\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_hgly2_cov <- dev_ermod_bin(\n data = df_er_ae_hgly2,\n var_resp = var_resp,\n var_cov = c(\"RACE\", \"BGLUC\", \"BHBA1C_5\"),\n var_exposure = var_exposure\n)\n```\n:::\n\n\n\n\n## Evaluation of covariate effects\n\nBy default, the covariate effect plots are generated for 5th and 95th percentiles\nof the continuous covariates and for each level of the categorical covariates\n(in the order of frequency in data).\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncoveffsim <- sim_coveff(ermod_bin_hgly2_cov)\nplot_coveff(coveffsim)\n```\n\n::: {.cell-output-display}\n![](coveff_customize_files/figure-html/fig-cov-eff-default-1.png){#fig-cov-eff-default width=576}\n:::\n:::\n\n\n\n\nYou can modify the width of the quantile interval using the `qi_width_cov` argument.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncoveffsim_qicov_08 <- sim_coveff(ermod_bin_hgly2_cov, qi_width_cov = 0.8)\nplot_coveff(coveffsim_qicov_08)\n```\n\n::: {.cell-output-display}\n![](coveff_customize_files/figure-html/fig-cov-eff-qicov_08-1.png){#fig-cov-eff-qicov_08 width=576}\n:::\n:::\n\n\n\n\n## Further customization\n\nThese covariate effect simulation (and plotting) can be customized by providing\nspecifications to `sim_coveff()` function.\n\nFirst, we build the specification for the covariate effects with\n`build_spec_coveff()` function.\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nspec_coveff <- build_spec_coveff(ermod_bin_hgly2_cov)\nspec_coveff |> \n gt() |> \n fmt_number(n_sigfig = 3)|>\n fmt_integer(columns = c(\"value_order\", \"var_order\"))\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n \n \n
var_ordervar_namevar_labelvalue_ordervalue_annotvalue_labelvalue_contvalue_catis_ref_valueshow_ref_valueis_covariate
1AUCss_1000AUCss_100015th0.8680.868NAFALSENAFALSE
1AUCss_1000AUCss_10002median2.212.21NATRUETRUEFALSE
1AUCss_1000AUCss_1000395th5.305.30NAFALSENAFALSE
2RACERACE11st freqWhiteNAWhiteTRUETRUETRUE
2RACERACE22nd freqAsianNAAsianFALSENATRUE
2RACERACE33rd freqBlackNABlackFALSENATRUE
3BGLUCBGLUC15th4.434.43NAFALSENATRUE
3BGLUCBGLUC2median6.106.10NATRUETRUETRUE
3BGLUCBGLUC395th7.597.59NAFALSENATRUE
4BHBA1C_5BHBA1C_515th5.755.75NAFALSENATRUE
4BHBA1C_5BHBA1C_52median7.977.97NATRUETRUETRUE
4BHBA1C_5BHBA1C_5395th10.410.4NAFALSENATRUE
\n
\n```\n\n:::\n:::\n\n\n\n\nLet's say we want to customize in a follwoing way:\n\n- Re-calculate percentiles of BGLUC from new distribution (example below uses\n uniform distribution between 4 and 8) and change the width of the quantile\n interval to 0.8.\n- Use specific values for BHBA1C (30, 35, 45, 50) and label with the original \n scale rather than those devided by 5 (used in the model).\n- Use Asian as the reference level\n- Show the plots in the order of BGLUC, BHBA1C, and then RACE\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nspec_new_bgluc <- build_spec_coveff_one_variable(\n \"BGLUC\", seq(4, 8, by = 0.1),\n qi_width_cov = 0.8, show_ref_value = TRUE\n)\n\nspec_new_bhba1c <- \n tibble(value_cont = c(30, 35, 45, 50) / 5) |> \n mutate(\n value_order = row_number(),\n value_label = as.character(value_cont * 5),\n var_name = \"BHBA1C_5\",\n var_label = \"BHBA1C\",\n value_annot = glue::glue(\"{value_label}mmol/mol\"),\n is_ref_value = FALSE,\n show_ref_value = FALSE)\n\nspec_new_race <- \n spec_coveff |> \n filter(var_name == \"RACE\") |> \n mutate(\n is_ref_value = c(FALSE, TRUE, FALSE),\n show_ref_value = c(FALSE, TRUE, FALSE)\n ) |> \n select(-var_order, -is_covariate)\n\nspec_coveff_new_1 <- \n replace_spec_coveff(spec_coveff, bind_rows(spec_new_bgluc, spec_new_bhba1c)) |>\n # spec_new_race is separately provided as we want to change the reference level\n replace_spec_coveff(spec_new_race, replace_ref_value = TRUE)\n\nd_new_var_order <- \n tibble(var_name = c(\"AUCss_1000\", \"BGLUC\", \"BHBA1C_5\", \"RACE\")) |> \n mutate(var_order = row_number())\n\nspec_coveff_new <- \n spec_coveff_new_1 |> \n select(-var_order) |>\n left_join(d_new_var_order, by = \"var_name\")\n\nspec_coveff_new |> \n gt() |> \n fmt_number(n_sigfig = 3)|>\n fmt_integer(columns = c(\"value_order\", \"var_order\"))\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n \n \n \n
var_namevalue_ordervalue_annotvalue_labelvalue_contvalue_catis_ref_valueshow_ref_valuevar_labelis_covariatevar_order
AUCss_100015th0.8680.868NAFALSENAAUCss_1000FALSE1
AUCss_10002median2.212.21NATRUETRUEAUCss_1000FALSE1
AUCss_1000395th5.305.30NAFALSENAAUCss_1000FALSE1
RACE11st freqWhiteNAWhiteFALSEFALSERACETRUE4
RACE22nd freqAsianNAAsianTRUETRUERACETRUE4
RACE33rd freqBlackNABlackFALSEFALSERACETRUE4
BGLUC1NA6.106.10NATRUEFALSEBGLUCTRUE2
BGLUC210th4.404.40NAFALSENABGLUCTRUE2
BGLUC3median6.006.00NAFALSENABGLUCTRUE2
BGLUC490th7.607.60NAFALSENABGLUCTRUE2
BHBA1C_51NA7.977.97NATRUEFALSEBHBA1CTRUE3
BHBA1C_5230mmol/mol306.00NAFALSENABHBA1CTRUE3
BHBA1C_5335mmol/mol357.00NAFALSENABHBA1CTRUE3
BHBA1C_5445mmol/mol459.00NAFALSENABHBA1CTRUE3
BHBA1C_5550mmol/mol5010.0NAFALSENABHBA1CTRUE3
\n
\n```\n\n:::\n:::\n\n\n\n\nThe customized covariate effect plots can be generated using the `plot_coveff()`\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ncoveffsim_spec <- sim_coveff(ermod_bin_hgly2_cov, spec_coveff = spec_coveff_new)\nplot_coveff(coveffsim_spec)\n```\n\n::: {.cell-output-display}\n![](coveff_customize_files/figure-html/fig-cov-eff-custom-1.png){#fig-cov-eff-custom width=576}\n:::\n:::\n", + "supporting": [ + "coveff_customize_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-custom-1.png b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-custom-1.png new file mode 100644 index 0000000..932060c Binary files /dev/null and b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-custom-1.png differ diff --git a/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-default-1.png b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-default-1.png new file mode 100644 index 0000000..66f4a24 Binary files /dev/null and b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-default-1.png differ diff --git a/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-qicov_08-1.png b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-qicov_08-1.png new file mode 100644 index 0000000..cdc49c0 Binary files /dev/null and b/_freeze/notebook/binary/coveff_customize/figure-html/fig-cov-eff-qicov_08-1.png differ diff --git a/_freeze/notebook/binary/mod_structure_comparison/execute-results/html.json b/_freeze/notebook/binary/mod_structure_comparison/execute-results/html.json new file mode 100644 index 0000000..601fb6b --- /dev/null +++ b/_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
\n── Binary ER model ─────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\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
\n── Binary Emax model ───────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\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
               elpd_diff se_diff\nbin_hgly2       0.0       0.0   \nbin_emax_hgly2 -1.7       1.7   \n
\n:::\n", + "supporting": [ + "mod_structure_comparison_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-emax-1.png b/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-emax-1.png new file mode 100644 index 0000000..cc7de07 Binary files /dev/null and b/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-emax-1.png differ diff --git a/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-lin-1.png b/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-lin-1.png new file mode 100644 index 0000000..04809d6 Binary files /dev/null and b/_freeze/notebook/binary/mod_structure_comparison/figure-html/fig-bin-lin-1.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/execute-results/html.json b/_freeze/notebook/binary/model_diagnostics/execute-results/html.json new file mode 100644 index 0000000..196a0cc --- /dev/null +++ b/_freeze/notebook/binary/model_diagnostics/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "b681bf31d4f6b3823db20602900e574b", + "result": { + "engine": "knitr", + "markdown": "# Model diagnostics and performance evaluation\n\nThis page showcase the model diagnosis and performance evaluation on the \nER model for binary endpoint.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(here)\nlibrary(posterior)\nlibrary(tidybayes)\nlibrary(bayesplot)\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}\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# Grade 2+ diarrhea\ndf_er_ae_dr2 <- d_sim_binom_cov_2 |> filter(AETYPE == \"dr2\")\n\nvar_resp <- \"AEFLAG\"\nvar_exposure <- \"AUCss_1000\"\n```\n:::\n\n\n\n\n## Fit models\n\nThere is clear trend of E-R for hyperglycemia (95% CI doesn't include 0) \nwhile the evidence of E-R is not seen for diarrhea (95% CI includes 0).\n\n::: {.panel-tabset}\n\n### Hyperglycemia\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)\nplot_er_gof(ermod_bin_hgly2, var_group = \"Dose\", show_coef_exp = TRUE)\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-data-fit-hgly2-1.png){#fig-data-fit-hgly2 width=672}\n:::\n:::\n\n\n\n\n### Diarrhea\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\nermod_bin_dr2 <- dev_ermod_bin(\n data = df_er_ae_dr2,\n var_resp = var_resp,\n var_exposure = var_exposure\n)\nplot_er_gof(ermod_bin_dr2, var_group = \"Dose\", show_coef_exp = TRUE)\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-data-fit-dr2-1.png){#fig-data-fit-dr2 width=672}\n:::\n:::\n\n\n\n\n:::\n\n## Parameter summary\n\nYou can see that AUCss effect is much stronger for hyperglycemia than diarrhea.\n\n::: {.panel-tabset}\n\n### Hyperglycemia\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nermod_bin_hgly2 |> \n summarize_draws(mean, median, sd, ~quantile2(.x, probs = c(0.025, 0.975)),\n default_convergence_measures()) |>\n gt() |>\n fmt_number(n_sigfig = 3)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdq2.5q97.5rhatess_bulkess_tail
(Intercept)βˆ’2.05βˆ’2.040.234βˆ’2.51βˆ’1.591.002,2301,870
AUCss_10000.4120.4120.07610.2650.5611.002,1502,140
\n
\n```\n\n:::\n:::\n\n\n\n\n### Diarrhea\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nermod_bin_dr2 |> \n summarize_draws(mean, median, sd, ~quantile2(.x, probs = c(0.025, 0.975)),\n default_convergence_measures()) |>\n gt() |>\n fmt_number(n_sigfig = 3)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdq2.5q97.5rhatess_bulkess_tail
(Intercept)βˆ’2.15βˆ’2.150.255βˆ’2.65βˆ’1.671.002,8302,320
AUCss_10000.1350.1340.0832βˆ’0.02930.3001.002,8102,570
\n
\n```\n\n:::\n:::\n\n\n\n\n:::\n\n## Predictive performance evaluation\n\nWe can calculate predictive performance metrics such as AUC-ROC with \n`eval_ermod()` function. Options for evaluation data are:\n\n- `eval_type = \"training\"`: training data\n- `eval_type = \"test\"`: test data (supply to `newdata` argument)\n- `eval_type = \"kfold\"`: k-fold cross-validation\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmetrics_hgly2_train <- eval_ermod(ermod_bin_hgly2, eval_type = \"training\")\nmetrics_hgly2_kfold <- eval_ermod(ermod_bin_hgly2, eval_type = \"kfold\")\n```\n:::\n\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmetrics_hgly2_train |>\n gt() |>\n fmt_number(n_sigfig = 3)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n\n\n \n\n\n \n \n \n
.metric.estimator.estimate
roc_aucbinary0.650
mn_log_lossbinary0.555
\n
\n```\n\n:::\n\n```{.r .cell-code}\nmetrics_hgly2_kfold |>\n gt() |>\n fmt_number(n_sigfig = 3) |>\n fmt_integer(columns = c(\"fold_id\"))\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n\n\n\n \n \n \n
fold_id.metric.estimator.estimate
1roc_aucbinary0.631
2roc_aucbinary0.607
3roc_aucbinary0.684
4roc_aucbinary0.670
5roc_aucbinary0.663
1mn_log_lossbinary0.602
2mn_log_lossbinary0.572
3mn_log_lossbinary0.534
4mn_log_lossbinary0.581
5mn_log_lossbinary0.500
\n
\n```\n\n:::\n:::\n\n\n\n\n## Probability of direction\n\nAlthough credible intervals are preferred, there is a concept called\nthe probability of direction which is somewhat similar to the p-value, in which\nthe probability of the effect being far from NULL (usually set to 0) is \ncalculated.\n\nSee `?p_direction` and [vignette](https://easystats.github.io/bayestestR/articles/probability_of_direction.html)\nfor detail.\n\n::: {.panel-tabset}\n\n### Hyperglycemia\n\nThe exposure effect is so clear that none of the MCMC sample is below 0, \nleading to a \"p-value\" of 0. \nSince there were 4000 MCMC samples (`nrow(as_draws_df(ermod_bin_hgly2))`),\nit is expected that the p-value is less than 1/4000 * 2, i.e. < 0.0005 \n(multiplication with 2 corresponds to two-sided test).\n\nkernel\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\np_direction(ermod_bin_hgly2, as_p = TRUE, as_num = TRUE)\n```\n\n
[1] 0\n
\n\n```{.r .cell-code}\n1 / length(as_draws(ermod_bin_hgly2)$AUCss_1000) * 2\n```\n\n
[1] 5e-04\n
\n:::\n\n\n\n\n### Diarrhea\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\np_direction(ermod_bin_dr2, as_p = TRUE, as_num = TRUE)\n```\n\n
[1] 0.108\n
\n\n```{.r .cell-code}\n# Below is a direct calculation of this value\nmean(as_draws(ermod_bin_dr2)$AUCss_1000 < 0) * 2\n```\n\n
[1] 0.108\n
\n:::\n\n\n\n:::\n\n\n## Diagnostic plots\n\nWe use the [`bayesplot` package](https://mc-stan.org/bayesplot/) \n([Cheat sheet](https://rstudio.github.io/cheatsheets/bayesplot.pdf)) \nto visualize the model fit.\n\n### Convergence\n\nGood fit results in:\n\n- Parameter distributions from MCMC chains should overlap\n- Trace plots should not show any trend\n- `Rhat` close to 1 (e.g. < 1.1)\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_draws_bin_hgly2 <- as_draws_df(ermod_bin_hgly2)\nmcmc_dens_overlay(d_draws_bin_hgly2)\nmcmc_trace(d_draws_bin_hgly2)\nmcmc_rhat(rhat(ermod_bin_hgly2$mod$stanfit))\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-1.png){#fig-convergence-1 width=672}\n:::\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-2.png){#fig-convergence-2 width=672}\n:::\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-3.png){#fig-convergence-3 width=672}\n:::\n:::\n\n\n\n\n### Parameter estimates distribution\n\n::: {.panel-tabset}\n\n#### `mcmc_hist`\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmcmc_hist(d_draws_bin_hgly2)\n```\n\n
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n
::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-mcmc_hist-1.png){#fig-mcmc_hist width=672}\n:::\n:::\n\n\n\n\n#### Parameter covariance\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmcmc_pairs(d_draws_bin_hgly2,\n off_diag_args = list(size = 0.5, alpha = 0.25))\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-mcmc_pairs-1.png){#fig-mcmc_pairs width=672}\n:::\n:::\n\n\n\n\n:::\n\n", + "supporting": [ + "model_diagnostics_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-1.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-1.png new file mode 100644 index 0000000..6da2824 Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-1.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-2.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-2.png new file mode 100644 index 0000000..737eeae Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-2.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-3.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-3.png new file mode 100644 index 0000000..5031d8e Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-convergence-3.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-dr2-1.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-dr2-1.png new file mode 100644 index 0000000..d7e64fb Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-dr2-1.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-hgly2-1.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-hgly2-1.png new file mode 100644 index 0000000..c629030 Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-data-fit-hgly2-1.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_hist-1.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_hist-1.png new file mode 100644 index 0000000..a3e3e71 Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_hist-1.png differ diff --git a/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_pairs-1.png b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_pairs-1.png new file mode 100644 index 0000000..4cff854 Binary files /dev/null and b/_freeze/notebook/binary/model_diagnostics/figure-html/fig-mcmc_pairs-1.png differ diff --git a/_freeze/notebook/binary/simulation/execute-results/html.json b/_freeze/notebook/binary/simulation/execute-results/html.json new file mode 100644 index 0000000..a0b2bda --- /dev/null +++ b/_freeze/notebook/binary/simulation/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "3f267b9fdda0626bfeb19f2f9983f353", + "result": { + "engine": "knitr", + "markdown": "# Simulation\n\nThis page showcase the model simulation using the ER model for binary endpoint.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(posterior)\nlibrary(tidybayes)\nlibrary(bayesplot)\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}\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 )\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\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
\n── Binary ER model ─────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\n:::\n\n\n\n\n\n## Prediction at specific conc\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nnew_conc_vec <- c(1, 3, 9)\n\n# Sim at specific conc\nd_sim_new_conc <-\n sim_er_new_exp(ermod_bin_hgly2,\n exposure_to_sim_vec = new_conc_vec,\n output_type = c(\"median_qi\"))\n\nd_sim_new_conc |>\n select(-starts_with(\".linpred\"), -c(.row, .width, .point, .interval)) |>\n gt() |>\n fmt_number(decimals = 2) |>\n tab_header(\n title = md(\"Predicted probability of events at specific concentrations\")\n )\n```\n\n::: {#fig-spec-conc-1 .cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n \n\n\n\n \n\n\n\n \n \n \n
Predicted probability of events at specific concentrations
AUCss_1000.epred.epred.lower.epred.upper
1.000.160.120.21
3.000.310.270.35
9.000.840.670.93
\n
\n```\n\n:::\n\n```{.r .cell-code}\n# Sim to draw ER curve\nd_sim_curve <-\n sim_er_curve(ermod_bin_hgly2, output_type = c(\"median_qi\"))\n\nd_sim_curve |>\n ggplot(aes(x = AUCss_1000, y = .epred)) +\n # Emax model curve\n geom_ribbon(aes(y = .epred, ymin = .epred.lower, ymax = .epred.upper),\n alpha = 0.3, fill = \"deepskyblue\") +\n geom_line(aes(y = .epred), color = \"deepskyblue\") +\n # Prediction at the specified doses\n geom_point(data = d_sim_new_conc, aes(y = .epred), color = \"tomato\", size = 3) +\n geom_errorbar(data = d_sim_new_conc,\n aes(y = .epred, ymin = .epred.lower, ymax = .epred.upper),\n width = 0.03, color = \"tomato\") +\n coord_cartesian(ylim = c(0, 1)) +\n scale_y_continuous(\n breaks = c(0, .5, 1),\n labels = scales::percent\n ) +\n labs(\n x = \"AUC~ss~ / 1000\", y = \"Probability of event\",\n title = \"ER model predictions at new exposure levels\",\n caption = \"Area: 95% credible interval\"\n ) +\n theme(axis.title.x = ggtext::element_markdown()) +\n xgxr::xgx_scale_x_log10()\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-spec-conc-1.png){#fig-spec-conc-2 width=672}\n:::\n:::\n\n\n\n\n## Prediction for new dose\n\n### Prep PK data for simulation\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\nd_new_dose_pk <- \n tibble(Dose_mg = rep(c(100, 200, 400), each = 30)) |> \n mutate(CL = rlnorm(n(), meanlog = log(100), sdlog = 0.3),\n AUCss_1000 = Dose_mg / CL,\n Dose = glue::glue(\"{Dose_mg} mg\"))\n\nd_median_auc <- \n d_new_dose_pk |> \n group_by(Dose) |> \n summarize(AUCss_1000 = median(AUCss_1000))\n\nggplot(d_new_dose_pk, aes(x = AUCss_1000, fill = Dose)) +\n geom_histogram(position = \"identity\", alpha = 0.5) +\n labs(\n x = \"AUC~ss~ / 1000\", y = \"Count\",\n title = \"Distribution of AUC~ss~ for new doses\"\n ) +\n theme(plot.title = ggtext::element_markdown(),\n axis.title.x = ggtext::element_markdown()) +\n xgxr::xgx_scale_x_log10()\n```\n\n
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n
::: {.cell-output-display}\n![](simulation_files/figure-html/fig-new-dose-data-1.png){#fig-new-dose-data width=672}\n:::\n:::\n\n\n\n\n### Sim and plot\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_sim_new_dose_raw <- \n sim_er(ermod_bin_hgly2,\n newdata = d_new_dose_pk)\n\nd_sim_new_dose_per_dose <- \n d_sim_new_dose_raw |> \n # Calc per-dose summary probability for each MCMC draws\n ungroup() |> \n summarize(prob = mean(.epred), .by = c(Dose, Dose_mg, .draw)) |> \n # Summarize across MCMC draws\n group_by(Dose, Dose_mg) |>\n median_qi() |> \n ungroup() |> \n # Add median AUCss\n left_join(d_median_auc, by = \"Dose\")\n\n\nd_sim_curve |>\n ggplot(aes(x = AUCss_1000, y = .epred)) +\n # Emax model curve\n geom_ribbon(aes(y = .epred, ymin = .epred.lower, ymax = .epred.upper),\n alpha = 0.3, fill = \"grey\") +\n geom_line(aes(y = .epred), color = \"grey\") +\n # Prediction at the specified doses\n geom_point(data = d_sim_new_dose_per_dose, aes(y = prob, color = Dose), size = 3) +\n geom_errorbar(data = d_sim_new_dose_per_dose,\n aes(y = prob, ymin = .lower, ymax = .upper, color = Dose), width = 0.03) +\n geom_boxplot(data = d_new_dose_pk, \n aes(x = AUCss_1000, y = -0.1, fill = Dose, color = Dose), width = 0.1, alpha = 0.5,\n inherit.aes = FALSE) +\n geom_hline(yintercept = 0, linetype = \"solid\", linewidth = 0.2) +\n coord_cartesian(ylim = c(-0.15, 1)) +\n scale_y_continuous(\n breaks = c(0, .5, 1),\n labels = scales::percent\n ) +\n labs(\n x = \"AUC~ss~ / 1000\", y = \"Probability of event\",\n title = \"ER model predictions at new dose levels\",\n caption = \"Area: 95% credible interval\n Boxplot: Observed exposure levels\n Symbols: Predicted mean probability for each dose and 95% CI\"\n ) +\n guides(\n fill = guide_legend(reverse = TRUE),\n color = guide_legend(reverse = TRUE)\n ) +\n theme(axis.title.x = ggtext::element_markdown()) +\n xgxr::xgx_scale_x_log10()\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-new-dose-sim-1.png){#fig-new-dose-sim width=672}\n:::\n:::\n\n\n\n\n## PK boxplot on ER plot\n\nWhen you use `plot_er_gof()` function, you can only add boxplot for the data\nthat you used for the model fit.\n\nSee below the example of adding exposure boxplot to any ER plot.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot_er(ermod_bin_hgly2, show_orig_data = TRUE) +\n geom_boxplot(data = d_new_dose_pk, \n aes(x = AUCss_1000, y = -0.2, fill = Dose, color = Dose), width = 0.1, alpha = 0.5,\n inherit.aes = FALSE) +\n geom_hline(yintercept = -0.1, linetype = \"solid\", linewidth = 0.2) + \n coord_cartesian(ylim = c(-0.25, 1)) +\n xgxr::xgx_scale_x_log10() +\n guides(\n fill = guide_legend(reverse = TRUE),\n color = guide_legend(reverse = TRUE)\n )\n```\n\n::: {.cell-output-display}\n![](simulation_files/figure-html/fig-er-and-obs-pk-1.png){#fig-er-and-obs-pk width=672}\n:::\n:::\n", + "supporting": [ + "simulation_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/simulation/figure-html/fig-er-and-obs-pk-1.png b/_freeze/notebook/binary/simulation/figure-html/fig-er-and-obs-pk-1.png new file mode 100644 index 0000000..6a5ee5e Binary files /dev/null and b/_freeze/notebook/binary/simulation/figure-html/fig-er-and-obs-pk-1.png differ diff --git a/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-data-1.png b/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-data-1.png new file mode 100644 index 0000000..76b5b08 Binary files /dev/null and b/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-data-1.png differ diff --git a/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-sim-1.png b/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-sim-1.png new file mode 100644 index 0000000..37b9092 Binary files /dev/null and b/_freeze/notebook/binary/simulation/figure-html/fig-new-dose-sim-1.png differ diff --git a/_freeze/notebook/binary/simulation/figure-html/fig-spec-conc-1.png b/_freeze/notebook/binary/simulation/figure-html/fig-spec-conc-1.png new file mode 100644 index 0000000..9c7b40e Binary files /dev/null and b/_freeze/notebook/binary/simulation/figure-html/fig-spec-conc-1.png differ diff --git a/_freeze/notebook/binary/workflow_wo_package/execute-results/html.json b/_freeze/notebook/binary/workflow_wo_package/execute-results/html.json new file mode 100644 index 0000000..b32ce0e --- /dev/null +++ b/_freeze/notebook/binary/workflow_wo_package/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "e32c76e0aa44e73501176e2d654324f0", + "result": { + "engine": "knitr", + "markdown": "# Workflow without using `BayesERtools`\n\nThis page shows how to perform ER analysis without using `BayesERtools` package\nto help:\n\n- Understand the internal workings of the package\n- Flexibility to customize the analysis\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(here)\nlibrary(gt)\n\ntheme_set(theme_bw(base_size = 12))\n```\n:::\n\n\n\n\n\n\n\n## Data\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndata(d_sim_binom_cov, package = \"BayesERtools\")\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 )\n\n# Grade 2+ hypoglycemia\ndf_er_ae_hgly2 <- d_sim_binom_cov_2 |> filter(AETYPE == \"hgly2\")\n\nvar_exposure <- \"AUCss_1000\"\nvar_resp <- \"AEFLAG\"\nvar_cov_ae_hgly2 <- c(\"BAGE_10\", \"BWT_10\", \"RACE\", \"VISC\", \"BHBA1C_5\", \"BGLUC\")\n```\n:::\n\n\n\n\n\n## Basic model development\n\n`dev_ermod_bin()` function can be used to develop basic ER model. (Note that\nthis function can also be used for models with covariates, if you already know\nthe covariates to be included.)\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\nvar_all <- c(var_exposure) # If you have covariates, you can add here\n\nformula_all <-\n stats::formula(glue::glue(\n \"{var_resp} ~ {paste(var_all, collapse = ' + ')}\"\n ))\n\nermod_bin <- rstanarm::stan_glm(\n formula_all,\n family = stats::binomial(),\n data = df_er_ae_hgly2,\n QR = dplyr::if_else(length(var_all) > 1, TRUE, FALSE),\n refresh = 0, # Suppress output\n)\n\nermod_bin\n```\n\n
stan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000\n observations: 500\n predictors:   2\n------\n            Median MAD_SD\n(Intercept) -2.0    0.2  \nAUCss_1000   0.4    0.1  \n\n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n
\n:::\n\n\n\n\nPerform simulation for plotting purpose\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nexposure_range <-\n range(df_er_ae_hgly2[[var_exposure]])\n\nexposure_to_sim_vec <-\n seq(exposure_range[1], exposure_range[2], length.out = 51)\n\ndata_for_sim <- \n tibble(!!var_exposure := exposure_to_sim_vec)\n\nsim_epred_med_qi <-\n tidybayes::add_epred_draws(data_for_sim, ermod_bin) |>\n tidybayes::median_qi() |>\n dplyr::as_tibble()\n```\n:::\n\n\n\n\nObserved vs model predicted plot:\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(data = sim_epred_med_qi, aes(x = .data[[var_exposure]], y = .epred)) +\n geom_ribbon(aes(ymin = .lower, ymax = .upper), alpha = 0.3) +\n geom_line() +\n # Observed data plots\n geom_jitter(\n data = df_er_ae_hgly2,\n aes(x = .data[[var_exposure]], y = .data[[var_resp]]),\n width = 0, height = 0.05, alpha = 0.5\n ) +\n xgxr::xgx_stat_ci(\n data = df_er_ae_hgly2,\n aes(x = .data[[var_exposure]], y = .data[[var_resp]]),\n bins = 4, conf_level = 0.95, distribution = \"binomial\",\n geom = c(\"point\"), shape = 0, size = 4\n ) +\n xgxr::xgx_stat_ci(\n data = df_er_ae_hgly2,\n aes(x = .data[[var_exposure]], y = .data[[var_resp]]),\n bins = 4, conf_level = 0.95, distribution = \"binomial\",\n geom = c(\"errorbar\"), linewidth = 0.5\n ) +\n # Figure settings\n coord_cartesian(ylim = c(-0.05, 1.05)) +\n scale_y_continuous(\n breaks = c(0, .5, 1),\n labels = scales::percent\n ) +\n labs(x = var_exposure, y = \"Probability of event\")\n```\n\n::: {.cell-output-display}\n![](workflow_wo_package_files/figure-html/fig-plot-er-1.png){#fig-plot-er width=672}\n:::\n:::\n\n\n\n\nMCMC samples can be obtained with `as_draws()` family of functions, such as\n`as_draws_df()`.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndraws_df <- posterior::as_draws_df(ermod_bin)\n\ndraws_df_summary <-\n posterior::summarize_draws(draws_df)\n\ndraws_df_summary |>\n gt::gt() |>\n gt::fmt_number(n_sigfig = 3)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdmadq5q95rhatess_bulkess_tail
(Intercept)βˆ’2.05βˆ’2.040.2340.228βˆ’2.43βˆ’1.671.002,2301,870
AUCss_10000.4120.4120.07610.07650.2880.5371.002,1502,140
\n
\n```\n\n:::\n:::\n\n\n\n\n\n## Selection of exposure metrics\n\nFirst you fit models with all the candidate exposure metrics and then compare\nthe models using leave-one-out cross-validation (LOO).\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\n# Run models with all the candidate exposure metrics\nl_mod_exposures <-\n c(\"AUCss_1000\", \"Cmaxss\", \"Cminss\") |>\n purrr::set_names() |>\n purrr::map(\n function(.x) {\n formula <- stats::formula(glue::glue(\"{var_resp} ~ {.x}\"))\n \n mod <- rstanarm::stan_glm(\n formula,\n family = stats::binomial(),\n data = df_er_ae_hgly2,\n refresh = 0 # Suppress output\n )\n },\n .progress = TRUE\n )\n\n# Calculate loo (leave-one-out cross-validation) for each model\n# and compare the models\nl_loo_exposures <- purrr::map(l_mod_exposures, loo::loo)\nloo::loo_compare(l_loo_exposures)\n```\n\n
           elpd_diff se_diff\nAUCss_1000  0.0       0.0   \nCminss     -4.4       3.1   \nCmaxss     -5.1       2.9   \n
\n:::\n\n\n\n\n\n## Selection of covariates\n\nSelection of covariates are be done with `projpred` package in `BayesERtools`.\n\n### Step 1: Full reference model fit\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nvarnames <- paste(c(var_exposure, var_cov_ae_hgly2), collapse = \" + \")\nformula_full <-\n stats::formula(\n glue::glue(\n \"{var_resp} ~ {varnames}\"\n )\n )\n\n# Need to construct call and then evaluate. Directly calling\n# rstanarm::stan_glm with formula_full does not work for the cross-validation\ncall_fit_ref <-\n rlang::call2(rstanarm::stan_glm,\n formula = formula_full,\n family = quote(stats::binomial()), data = quote(df_er_ae_hgly2), QR = TRUE,\n refresh = 0)\nfit_ref <- eval(call_fit_ref)\n\nrefm_obj <- projpred::get_refmodel(fit_ref)\n```\n:::\n\n\n\n\n### Step 2: Variable selection\n\nThe code below shows example of variable selection with K-fold cross-validation\napproach.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Force exposure metric to be always included first\nsearch_terms <- projpred::force_search_terms(\n forced_terms = var_exposure,\n optional_terms = var_cov_ae_hgly2\n)\n\ncvvs <- projpred::cv_varsel(\n refm_obj,\n search_terms = search_terms,\n cv_method = \"kfold\",\n method = \"forward\",\n validate_search = TRUE,\n refit_prj = TRUE # Evaluation often look strange without refit\n)\n```\n\n
-----\nRunning the search using the full dataset ...\n10% of terms selected\n20% of terms selected\n40% of terms selected\n50% of terms selected\n70% of terms selected\n80% of terms selected\n100% of terms selected\n-----\n-----\nRefitting the reference model K = 5 times (using the fold-wise training data) ...\n
Fitting model 1 out of 5\n
Fitting model 2 out of 5\n
Fitting model 3 out of 5\n
Fitting model 4 out of 5\n
Fitting model 5 out of 5\n
-----\n-----\nRunning the search and the performance evaluation with `refit_prj = TRUE` for each of the K = 5 CV folds separately ...\n\n  |                                                                            \n  |                                                                      |   0%\n  |                                                                            \n  |==============                                                        |  20%\n  |                                                                            \n  |============================                                          |  40%\n  |                                                                            \n  |==========================================                            |  60%\n  |                                                                            \n  |========================================================              |  80%\n  |                                                                            \n  |======================================================================| 100%\n-----\n
\n\n```{.r .cell-code}\nrk <- projpred::ranking(cvvs)\n\nn_var_select <- projpred::suggest_size(cvvs)\nn_var_select <- max(1, n_var_select) # At least exposure metric should be included\n\nvar_selected <- head(rk[[\"fulldata\"]], n_var_select)\n```\n:::\n\n\n\n\n#### Output\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nvar_selected\n```\n\n
[1] \"AUCss_1000\" \"BHBA1C_5\"   \"BGLUC\"     \n
\n\n```{.r .cell-code}\nplot(cvvs, text_angle = 45, show_cv_proportions = FALSE, deltas = TRUE)\nplot(rk) # This only works when cv_method = \"kfold\" and validate_search = TRUE\n```\n\n::: {.cell-output-display}\n![](workflow_wo_package_files/figure-html/fig-covsel-1.png){#fig-covsel-1 width=672}\n:::\n\n::: {.cell-output-display}\n![](workflow_wo_package_files/figure-html/fig-covsel-2.png){#fig-covsel-2 width=672}\n:::\n:::\n\n\n\n\n\n### Step 3: Final model fit\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nset.seed(1234)\n\nermod_bin_cov <- rstanarm::stan_glm(\n stats::formula(glue::glue(\n \"{var_resp} ~ {paste(var_selected, collapse = ' + ')}\"\n )),\n family = stats::binomial(),\n data = df_er_ae_hgly2,\n QR = dplyr::if_else(length(var_selected) > 1, TRUE, FALSE),\n refresh = 0, # Suppress output\n)\n\nermod_bin_cov\n```\n\n
stan_glm\n family:       binomial [logit]\n formula:      AEFLAG ~ AUCss_1000 + BHBA1C_5 + BGLUC\n observations: 500\n predictors:   4\n------\n            Median MAD_SD\n(Intercept) -11.0    1.2 \nAUCss_1000    0.5    0.1 \nBHBA1C_5      0.5    0.1 \nBGLUC         0.7    0.1 \n\n------\n* For help interpreting the printed output see ?print.stanreg\n* For info on the priors used see ?prior_summary.stanreg\n
\n:::\n\n\n\n\n## Marginal ER prediction\n\nThe example below simulate the marginal ER relationship, i.e. \nER relationships for \"marginalized\", or averaged, response for the population \nof interest, using the covariate distribution is from the original data.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nexposure_to_sim <- c(2:6)\n\ndata_cov <- df_er_ae_hgly2 |> select(-!!var_exposure)\n\ndata_for_sim <- \n tibble(!!var_exposure := exposure_to_sim) |>\n mutate(.id_exposure = row_number()) |>\n expand_grid(data_cov)\n\nsim_epred_raw <-\n tidybayes::add_epred_draws(data_for_sim, ermod_bin_cov)\n\n# Calculate marginal expected response for each exposure value and draw\nsim_epred_marg <-\n sim_epred_raw |>\n dplyr::ungroup() |>\n dplyr::summarize(\n .epred = mean(.epred),\n .by = c(.id_exposure, !!var_exposure, .draw)\n )\n\nsim_epred_marg_med_qi <-\n sim_epred_marg |>\n dplyr::group_by(.id_exposure, !!sym(var_exposure)) |>\n tidybayes::median_qi() |>\n dplyr::as_tibble()\n\nsim_epred_marg_med_qi |>\n gt::gt() |>\n gt::fmt_number(n_sigfig = 3) |>\n gt::fmt_integer(columns = c(\".id_exposure\"))\n```\n\n::: {#fig-marginal-er .cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n \n \n \n
.id_exposureAUCss_1000.epred.lower.upper.width.point.interval
12.000.2280.1920.2670.950medianqi
23.000.3070.2680.3470.950medianqi
34.000.3960.3410.4550.950medianqi
45.000.4910.4090.5760.950medianqi
56.000.5880.4770.6920.950medianqi
\n
\n```\n\n:::\n:::\n", + "supporting": [ + "workflow_wo_package_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-1.png b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-1.png new file mode 100644 index 0000000..eee1359 Binary files /dev/null and b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-1.png differ diff --git a/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-2.png b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-2.png new file mode 100644 index 0000000..ca61eca Binary files /dev/null and b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-covsel-2.png differ diff --git a/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-plot-er-1.png b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-plot-er-1.png new file mode 100644 index 0000000..6756507 Binary files /dev/null and b/_freeze/notebook/binary/workflow_wo_package/figure-html/fig-plot-er-1.png differ diff --git a/_freeze/notebook/emax/basic_workflow/execute-results/html.json b/_freeze/notebook/emax/basic_workflow/execute-results/html.json new file mode 100644 index 0000000..ef0c984 --- /dev/null +++ b/_freeze/notebook/emax/basic_workflow/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "702afada72c3b9fe128d1858fb2c2920", + "result": { + "engine": "knitr", + "markdown": "# Basic workflow\n\nIn this section, we will show a basic workflow for performing an E~max~ model \nanalysis for continuous endpoint.\n\n## Setup and load\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(tidyverse)\nlibrary(BayesERtools)\nlibrary(posterior)\nlibrary(tidybayes)\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}\nd_example_emax_nocov <- read_csv(here(\"data\", \"d_example_emax_nocov.csv\"))\n\nd_example_emax_nocov |>\n head() |>\n gt() |>\n fmt_number(decimals = 2)\n```\n\n::: {#fig-data-look-1 .cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n\n \n\n \n\n \n\n \n\n \n\n \n \n \n
ConcY
47.5115.14
12.447.26
14.908.78
13.627.05
29.1610.82
45.1613.27
\n
\n```\n\n:::\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```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-data-look-1.png){#fig-data-look-2 width=672}\n:::\n:::\n\n\n\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
\n── Emax model ──────────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\n:::\n\n\n\n\n## Observation vs model fit\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_sim_ermod_sigemax <-\n sim_er(ermod_sigemax, output_type = c(\"median_qi\"))\n\nplot_er_gof(d_sim_ermod_sigemax)\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-plot-er-1.png){#fig-plot-er width=672}\n:::\n:::\n\n\n\n\n## Parameter estimates\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_draws_sigemax <- as_draws_df(ermod_sigemax)\n\nd_draws_sigemax_summary <-\n summarize_draws(d_draws_sigemax)\n\nec50_mean <-\n d_draws_sigemax_summary |>\n filter(variable == \"ec50\") |>\n pull(mean)\n\nd_draws_sigemax_summary |>\n gt() |>\n fmt_number(decimals = 2)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdmadq5q95rhatess_bulkess_tail
ec5024.9422.518.635.6516.1742.091.001,111.481,157.42
sigma0.860.840.160.150.641.151.001,434.701,777.17
gamma1.771.620.770.680.843.121.00593.05951.82
e03.133.582.001.62βˆ’0.845.381.00675.84809.63
emax15.9214.715.495.139.3526.511.00583.13809.81
\n
\n```\n\n:::\n:::\n\n\n\n\n\n## Prediction at a certain concentrations\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_sim_new_conc <-\n sim_er_new_exp(ermod_sigemax,\n exposure_to_sim_vec = c(10, 20, 30, 50),\n output_type = c(\"median_qi\"))\n\nd_sim_new_conc |>\n select(-starts_with(\".linpred\"), -c(.row, .width, .point, .interval)) |>\n gt() |>\n fmt_number(decimals = 2) |>\n tab_header(\n title = md(\"Prediction at specific concentrations\")\n )\n```\n\n::: {#fig-plot-er-new-conc-1 .cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n \n\n\n\n\n\n\n \n\n\n\n\n\n\n \n\n\n\n\n\n\n \n \n \n
Prediction at specific concentrations
Conc.epred.epred.lower.epred.upper.prediction.prediction.lower.prediction.upper
10.006.575.957.166.564.778.40
20.0010.049.3910.7710.048.2411.97
30.0012.3811.7813.0212.3810.4714.34
50.0014.7613.5515.9314.7512.5016.78
\n
\n```\n\n:::\n\n```{.r .cell-code}\nd_sim_ermod_sigemax |>\n ggplot(aes(x = Conc, y = Y)) +\n # Emax model curve\n geom_vline(xintercept = ec50_mean, linetype = \"dashed\", color = \"deepskyblue\") +\n geom_ribbon(aes(y = .epred, ymin = .epred.lower, ymax = .epred.upper),\n alpha = 0.3, fill = \"deepskyblue\") +\n geom_line(aes(y = .epred), color = \"deepskyblue\") +\n # Observed data\n geom_point(data = d_example_emax_nocov, color = \"grey\") +\n # Prediction at the specified doses\n geom_point(data = d_sim_new_conc, aes(y = .epred), color = \"tomato\", size = 3) +\n geom_errorbar(data = d_sim_new_conc,\n aes(y = .epred, ymin = .epred.lower, ymax = .epred.upper),\n width = 1, color = \"tomato\") +\n coord_cartesian(ylim = c(-1, 17)) +\n scale_fill_brewer(palette = \"Greys\") +\n labs(\n title = \"Sigmoidal E~max~ model predictions at new exposure levels\",\n caption =\n \"vertical dashed line: estimated EC~50~ value
area: 95% credible interval\") +\n theme(plot.title = ggtext::element_markdown(),\n plot.caption = ggtext::element_markdown())\n```\n\n::: {.cell-output-display}\n![](basic_workflow_files/figure-html/fig-plot-er-new-conc-1.png){#fig-plot-er-new-conc-2 width=672}\n:::\n:::\n", + "supporting": [ + "basic_workflow_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/emax/basic_workflow/figure-html/fig-data-look-1.png b/_freeze/notebook/emax/basic_workflow/figure-html/fig-data-look-1.png new file mode 100644 index 0000000..dabcb85 Binary files /dev/null and b/_freeze/notebook/emax/basic_workflow/figure-html/fig-data-look-1.png differ diff --git a/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-1.png b/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-1.png new file mode 100644 index 0000000..3e879e6 Binary files /dev/null and b/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-1.png differ diff --git a/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-new-conc-1.png b/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-new-conc-1.png new file mode 100644 index 0000000..29da186 Binary files /dev/null and b/_freeze/notebook/emax/basic_workflow/figure-html/fig-plot-er-new-conc-1.png differ diff --git a/_freeze/notebook/emax/data_generation/execute-results/html.json b/_freeze/notebook/emax/data_generation/execute-results/html.json new file mode 100644 index 0000000..a9da12b --- /dev/null +++ b/_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 + } +} \ No newline at end of file diff --git a/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-cov-1.png b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-cov-1.png new file mode 100644 index 0000000..c011b5c Binary files /dev/null and b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-cov-1.png differ diff --git a/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-1.png b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-1.png new file mode 100644 index 0000000..2a0e5db Binary files /dev/null and b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-1.png differ diff --git a/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-2.png b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-2.png new file mode 100644 index 0000000..4fa3a5f Binary files /dev/null and b/_freeze/notebook/emax/data_generation/figure-html/fig-check-data-no-cov-2.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/execute-results/html.json b/_freeze/notebook/emax/model_diagnostics/execute-results/html.json new file mode 100644 index 0000000..f22dd59 --- /dev/null +++ b/_freeze/notebook/emax/model_diagnostics/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "2772926bb500d62c529cb1e78d9cc7c7", + "result": { + "engine": "knitr", + "markdown": "# Model diagnostics & comparison\n\nThis page showcase the model diagnosis and comparison for the E~max~ model \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
\n── Emax model ──────────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\n:::\n\n\n\n\n### Regular E~max~ (Ξ³ 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
\n── Emax model ──────────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\n:::\n\n\n\n\n:::\n\n## Parameter summary table\n\n::: {.panel-tabset}\n\n### `posterior` package\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_draws_sigemax_summary <-\n summarize_draws(ermod_sigemax)\n\nd_draws_sigemax_summary |>\n gt() |>\n fmt_number(decimals = 2)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n \n \n \n
variablemeanmediansdmadq5q95rhatess_bulkess_tail
ec5024.9422.518.635.6516.1742.091.001,111.481,157.42
sigma0.860.840.160.150.641.151.001,434.701,777.17
gamma1.771.620.770.680.843.121.00593.05951.82
e03.133.582.001.62βˆ’0.845.381.00675.84809.63
emax15.9214.715.495.139.3526.511.00583.13809.81
\n
\n```\n\n:::\n:::\n\n\n\n\n### Highest density interval\n\nHere is the example of highest-density continuous interval (HDCI) for the median of `ED50`.\nSee [here](https://mjskay.github.io/ggdist/reference/point_interval.html) for more details.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# HDCI of median.ED50\nas_draws_df(ermod_sigemax) |>\n tidybayes::spread_rvars(ec50) |>\n tidybayes::median_hdci()\n```\n\n
# A tibble: 1 Γ— 6\n   ec50 .lower .upper .width .point .interval\n               \n1  22.5   13.3   43.3   0.95 median hdci     \n
\n:::\n\n\n\n\n:::\n\n\n## Fitted values\n\nFitted values without residual errors (i.e. PRED in NONMEM term) can be\nextracted with `sim_er()` function. `.epred` is the expected value prediction.\nSee `?sim_er` for detail.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsim_er(ermod_sigemax) |> head()\n```\n\n
# A tibble: 6 Γ— 7\n# Groups:   Conc, Y, .row [1]\n   Conc     Y  .row .draw    .epred .prediction  .linpred\n              \n1  47.5  15.1     1     1      14.5        14.7      14.5\n2  47.5  15.1     1     2      14.6        16.8      14.6\n3  47.5  15.1     1     3      14.7        13.5      14.7\n4  47.5  15.1     1     4      15.0        15.2      15.0\n5  47.5  15.1     1     5      14.4        13.9      14.4\n6  47.5  15.1     1     6      14.5        14.7      14.5\n
\n:::\n\n\n\n\nYou can specify `output_type = \"median_qi\"` to get median and quantile\nintervals of the prediction.\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_sigemax_med_qi <-\n sim_er(ermod_sigemax, output_type = \"median_qi\")\n\nersim_sigemax_med_qi |>\n arrange(.row) |>\n head() |>\n gt(rownames_to_stub = TRUE) |>\n fmt_number(decimals = 2, columns = -.row)\n```\n\n::: {.cell-output-display}\n\n```{=html}\n
\n\n\n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n \n \n \n
ConcY.row.epred.epred.lower.epred.upper.prediction.prediction.lower.prediction.upper.linpred.linpred.lower.linpred.upper.width.point.interval
147.5115.14114.5613.4715.6114.5412.5516.5714.5613.4715.610.95medianqi
212.447.2627.506.818.117.495.629.307.506.818.110.95medianqi
314.908.7838.407.699.028.366.5310.208.407.699.020.95medianqi
413.627.0547.947.248.567.936.109.837.947.248.560.95medianqi
529.1610.82512.2211.6212.8712.2010.2814.0412.2211.6212.870.95medianqi
645.1613.27614.3613.3615.2714.3412.3616.2914.3613.3615.270.95medianqi
\n
\n```\n\n:::\n:::\n\n\n\n\n\n## Diagnostic plots\n\nWe use the [`bayesplot` package](https://mc-stan.org/bayesplot/) \n([Cheat sheet](https://rstudio.github.io/cheatsheets/bayesplot.pdf)) \nto visualize the model fit.\n\n### Convergence\n\nGood fit results in:\n\n- Parameter distributions from MCMC chains should overlap\n- Trace plots should not show any trend\n- `Rhat` close to 1 (e.g. < 1.1)\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nd_draws_sigemax <- as_draws_df(ermod_sigemax)\nmcmc_dens_overlay(d_draws_sigemax)\nmcmc_trace(d_draws_sigemax)\nmcmc_rhat(rhat(ermod_sigemax$mod$stanfit))\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-1.png){#fig-convergence-1 width=672}\n:::\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-2.png){#fig-convergence-2 width=672}\n:::\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-convergence-3.png){#fig-convergence-3 width=672}\n:::\n:::\n\n\n\n\n### Parameter estimates distribution\n\n::: {.panel-tabset}\n\n#### `mcmc_hist`\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmcmc_hist(d_draws_sigemax)\n```\n\n
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.\n
::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-mcmc_hist-1.png){#fig-mcmc_hist width=672}\n:::\n:::\n\n\n\n\n#### Parameter covariance\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nmcmc_pairs(d_draws_sigemax,\n off_diag_args = list(size = 0.5, alpha = 0.25))\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-mcmc_pairs-1.png){#fig-mcmc_pairs width=672}\n:::\n:::\n\n\n\n\n:::\n\n### GOF plots\n\n::: {.panel-tabset}\n\n#### Obs & Pred vs predictor\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot_er(ermod_sigemax, show_orig_data = TRUE)\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-idv_vs_pred-1.png){#fig-idv_vs_pred width=672}\n:::\n:::\n\n\n\n\n#### Obs vs pred\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\n# Preparation for diagnostic plots\nersim_sigemax_med_qi <- sim_er(ermod_sigemax, output_type = \"median_qi\")\n\nersim_sigemax_med_qi |>\n ggplot(aes(.epred, Y)) +\n geom_abline(linetype = 2, color = \"grey\") +\n geom_point() +\n geom_errorbarh(aes(xmin = .epred.lower, xmax = .epred.upper), height = 0) +\n labs(title = \"Observed vs Predicted\",\n x = \"Predicted\",\n y = \"Observed\",\n caption = \"Symbol: median and 95% credible interval\")\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-obs_vs_pred-1.png){#fig-obs_vs_pred width=576}\n:::\n:::\n\n\n\n\n#### Residuals\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_sigemax_w_resid <-\n sim_er(ermod_sigemax) |>\n mutate(.residual = Y - .epred) # Add residuals for plotting\nersim_sigemax_w_resid_med_qi <- median_qi(ersim_sigemax_w_resid)\n\nersim_sigemax_w_resid_med_qi |>\n ggplot(aes(x = .epred, y = .residual)) +\n xlab(\"Predicted (linear)\") +\n ylab(\"Residuals\") +\n geom_point() +\n geom_errorbar(aes(ymin = .residual.lower, ymax = .residual.upper), width = 0) +\n geom_hline(aes(yintercept = 2), lty = 2, colour = \"grey70\") +\n geom_hline(aes(yintercept = -2), lty = 2, colour = \"grey70\")\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-residuals-1.png){#fig-residuals width=672}\n:::\n:::\n\n\n\n\n#### Q-Q plot of residuals\n\n\n\n\n::: {.cell}\n\n```{.r .cell-code}\nersim_sigemax_w_resid_med_qi |>\n ggplot(aes(sample = .residual)) +\n geom_qq() +\n geom_qq_line(colour = \"steelblue\", lty = 2, alpha = 0.4) +\n coord_equal() +\n xlab(\"Theoretical\") +\n ylab(\"Sample\")\n```\n\n::: {.cell-output-display}\n![](model_diagnostics_files/figure-html/fig-qq-1.png){#fig-qq 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 E~max~ model with `Ξ³` fixed to be 1 appears better.\nHowever, `elpd_diff` is tiny and smaller than `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\n::: {.cell}\n\n```{.r .cell-code}\nloo_sigemax <- loo(ermod_sigemax)\nloo_emax <- loo(ermod_emax)\n\nloo_compare(list(sigemax = loo_sigemax, emax = loo_emax))\n```\n\n
        elpd_diff se_diff\nemax     0.0       0.0   \nsigemax -0.3       0.8   \n
\n:::\n", + "supporting": [ + "model_diagnostics_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-1.png new file mode 100644 index 0000000..e8ae22d Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-2.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-2.png new file mode 100644 index 0000000..ddf9ed1 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-2.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-3.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-3.png new file mode 100644 index 0000000..be19345 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-convergence-3.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-idv_vs_pred-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-idv_vs_pred-1.png new file mode 100644 index 0000000..4dbfa8d Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-idv_vs_pred-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_hist-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_hist-1.png new file mode 100644 index 0000000..8ce3d19 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_hist-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_pairs-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_pairs-1.png new file mode 100644 index 0000000..fe2c310 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-mcmc_pairs-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-obs_vs_pred-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-obs_vs_pred-1.png new file mode 100644 index 0000000..06db95c Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-obs_vs_pred-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-qq-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-qq-1.png new file mode 100644 index 0000000..f3a0bf9 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-qq-1.png differ diff --git a/_freeze/notebook/emax/model_diagnostics/figure-html/fig-residuals-1.png b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-residuals-1.png new file mode 100644 index 0000000..13df889 Binary files /dev/null and b/_freeze/notebook/emax/model_diagnostics/figure-html/fig-residuals-1.png differ diff --git a/_freeze/notebook/emax/simulation/execute-results/html.json b/_freeze/notebook/emax/simulation/execute-results/html.json new file mode 100644 index 0000000..a5aa75f --- /dev/null +++ b/_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
\n── Emax model ──────────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\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
\n── Emax model ──────────────────────────────────────────────────────────────────\nβ„Ή Use `plot_er()` to visualize ER curve\n\n── Developed model ──\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
\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 + } +} \ No newline at end of file diff --git a/_freeze/notebook/emax/simulation/figure-html/fig-emax-sigemax-compare-1.png b/_freeze/notebook/emax/simulation/figure-html/fig-emax-sigemax-compare-1.png new file mode 100644 index 0000000..1f26a0c Binary files /dev/null and b/_freeze/notebook/emax/simulation/figure-html/fig-emax-sigemax-compare-1.png differ diff --git a/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-1-1.png b/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-1-1.png new file mode 100644 index 0000000..66f2408 Binary files /dev/null and b/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-1-1.png differ diff --git a/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-2-1.png b/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-2-1.png new file mode 100644 index 0000000..4e19633 Binary files /dev/null and b/_freeze/notebook/emax/simulation/figure-html/fig-extrapolation-2-1.png differ diff --git a/_freeze/site_libs/clipboard/clipboard.min.js b/_freeze/site_libs/clipboard/clipboard.min.js new file mode 100644 index 0000000..1103f81 --- /dev/null +++ b/_freeze/site_libs/clipboard/clipboard.min.js @@ -0,0 +1,7 @@ +/*! + * clipboard.js v2.0.11 + * https://clipboardjs.com/ + * + * Licensed MIT Β© Zeno Rocha + */ +!function(t,e){"object"==typeof exports&&"object"==typeof module?module.exports=e():"function"==typeof define&&define.amd?define([],e):"object"==typeof exports?exports.ClipboardJS=e():t.ClipboardJS=e()}(this,function(){return n={686:function(t,e,n){"use strict";n.d(e,{default:function(){return b}});var e=n(279),i=n.n(e),e=n(370),u=n.n(e),e=n(817),r=n.n(e);function c(t){try{return document.execCommand(t)}catch(t){return}}var a=function(t){t=r()(t);return c("cut"),t};function o(t,e){var n,o,t=(n=t,o="rtl"===document.documentElement.getAttribute("dir"),(t=document.createElement("textarea")).style.fontSize="12pt",t.style.border="0",t.style.padding="0",t.style.margin="0",t.style.position="absolute",t.style[o?"right":"left"]="-9999px",o=window.pageYOffset||document.documentElement.scrollTop,t.style.top="".concat(o,"px"),t.setAttribute("readonly",""),t.value=n,t);return e.container.appendChild(t),e=r()(t),c("copy"),t.remove(),e}var f=function(t){var e=1