Add regression plot

ndunnewind · Apr 26, 2024 · fb3ef46 · fb3ef46
1 parent 0b094dd
commit fb3ef46
Show file tree

Hide file tree

Showing 2 changed files with 108 additions and 24 deletions.
diff --git a/tests/testthat/test-example_certolizumab.R b/tests/testthat/test-example_certolizumab.R
@@ -8,44 +8,44 @@ skip_on_cran()
 params <-
 list(run_tests = FALSE)
 
-## ----code=readLines("children/knitr_setup.R"), include=FALSE-------------------------------------------------------------------------------------
+## ----code=readLines("children/knitr_setup.R"), include=FALSE--------------------------------------
 
-## ----include=FALSE-------------------------------------------------------------------------------------------------------------------------------
+## ----include=FALSE--------------------------------------------------------------------------------
 set.seed(76439)
 
 
-## ----eval = FALSE--------------------------------------------------------------------------------------------------------------------------------
+## ----eval = FALSE---------------------------------------------------------------------------------
 ## library(multinma)
 ## options(mc.cores = parallel::detectCores())
 
-## ----setup, echo = FALSE-------------------------------------------------------------------------------------------------------------------------
+## ----setup, echo = FALSE--------------------------------------------------------------------------
 library(multinma)
-nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")),
-             "true" =, "warn" = 2,
+nc <- switch(tolower(Sys.getenv("_R_CHECK_LIMIT_CORES_")), 
+             "true" =, "warn" = 2, 
              parallel::detectCores())
 options(mc.cores = nc)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 head(certolizumab)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 cert_net <- set_agd_arm(certolizumab,
                         study = study, trt = trt, n = n, r = r,
                         trt_class = dplyr::if_else(trt == "Placebo", "Placebo", "Treatment"))
 cert_net
 
 
-## ----eval=FALSE----------------------------------------------------------------------------------------------------------------------------------
+## ----eval=FALSE-----------------------------------------------------------------------------------
 ## plot(cert_net, weight_edges = TRUE, weight_nodes = TRUE)
 
-## ----certolizumab_network_plot, echo=FALSE-------------------------------------------------------------------------------------------------------
+## ----certolizumab_network_plot, echo=FALSE--------------------------------------------------------
 plot(cert_net, weight_edges = TRUE, weight_nodes = TRUE) +
   ggplot2::theme(legend.box.margin = ggplot2::unit(c(0, 0, 0, 4), "lines"))
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 cert_fit_FE <- nma(cert_net,
                    trt_effects = "fixed",
                    regression = ~.mu:.trt,
@@ -56,11 +56,11 @@ cert_fit_FE <- nma(cert_net,
                    adapt_delta = 0.95)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 cert_fit_FE
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 cert_fit_RE <- nma(cert_net,
                    trt_effects = "random",
                    regression = ~.mu:.trt,
@@ -70,50 +70,92 @@ cert_fit_RE <- nma(cert_net,
                    iter = 4000)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 cert_fit_RE
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 (dic_FE <- dic(cert_fit_FE))
 (dic_RE <- dic(cert_fit_RE))
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 plot(dic_FE)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 plot(dic_RE)
 
 
-## ------------------------------------------------------------------------------------------------------------------------------------------------
+## -------------------------------------------------------------------------------------------------
 newdata <- data.frame(.mu = cert_fit_FE$xbar[[".mu"]])
 
 
-## ----certolizumab_releff_FE, fig.height=3--------------------------------------------------------------------------------------------------------
+## ----certolizumab_releff_FE, fig.height=3---------------------------------------------------------
 (cert_releff_FE <- relative_effects(cert_fit_FE, newdata = newdata))
 plot(cert_releff_FE, ref_line = 0)
 
-## ----certolizumab_releff_RE, fig.height=3--------------------------------------------------------------------------------------------------------
+## ----certolizumab_releff_RE, fig.height=3---------------------------------------------------------
 (cert_releff_RE <- relative_effects(cert_fit_RE, newdata = newdata))
 plot(cert_releff_RE, ref_line = 0)
 
 
-## ----certolizumab_ranks--------------------------------------------------------------------------------------------------------------------------
+## ----certolizumab_ranks---------------------------------------------------------------------------
 (cert_ranks <- posterior_ranks(cert_fit_RE, newdata = newdata))
 plot(cert_ranks)
 
-## ----certolizumab_rankprobs----------------------------------------------------------------------------------------------------------------------
+## ----certolizumab_rankprobs-----------------------------------------------------------------------
 (cert_rankprobs <- posterior_rank_probs(cert_fit_RE, newdata = newdata))
 plot(cert_rankprobs)
 
-## ----certolizumab_cumrankprobs-------------------------------------------------------------------------------------------------------------------
+## ----certolizumab_cumrankprobs--------------------------------------------------------------------
 (cert_cumrankprobs <- posterior_rank_probs(cert_fit_RE, cumulative = TRUE, newdata = newdata))
 plot(cert_cumrankprobs)
 
 
-## ----certolizumab_tests, include=FALSE, eval=params$run_tests------------------------------------------------------------------------------------
+## ----certolizumab_reg_plot------------------------------------------------------------------------
+library(dplyr)
+library(ggplot2)
+
+mu_dat <- seq(log(0.02), log(0.5), length.out = 20)
+
+cert_mu_reg <- cert_fit_FE %>%
+  relative_effects(newdata = tibble::tibble(.mu = mu_dat), study = .mu) %>%
+  as_tibble() %>%
+  mutate(
+    .trt = .trtb,
+    .mu = as.numeric(as.character(.study))
+  )
+
+cert_lor <-
+  cert_net$agd_arm %>%
+  mutate(
+    probability = .r / .n,
+    log_odds = log(probability / (1 - probability))
+  ) %>%
+  group_by(.study) %>%
+  mutate(
+    .mu = log_odds[.trt == "Placebo"],
+    log_odds_ratio = log_odds - .mu,
+    sample_size = sum(n),
+  ) %>%
+  ungroup() %>%
+  filter(.trt != "Placebo")
+
+cert_lor %>%
+  ggplot(aes(x = .mu)) +
+  geom_hline(yintercept = 0, colour = "grey60") +
+  geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), data = cert_mu_reg,
+              fill = "darkred", alpha = 0.3) +
+  geom_line(aes(y = mean), data = cert_mu_reg,
+            colour = "darkred") +
+  geom_point(aes(y = log_odds_ratio, size = sample_size), alpha = 0.6) +
+  facet_wrap(vars(.trt)) +
+  labs(x = "Placebo log odds", y = "log Odds Ratio", size = "Sample Size") +
+  theme_multinma()
+
+
+## ----certolizumab_tests, include=FALSE, eval=params$run_tests-------------------------------------
 #--- Test against TSD 3 results ---
 library(testthat)
 library(dplyr)

diff --git a/vignettes/example_certolizumab.Rmd b/vignettes/example_certolizumab.Rmd
@@ -130,6 +130,48 @@ plot(cert_rankprobs)
 plot(cert_cumrankprobs)
 ```
 
+```{r certolizumab_reg_plot}
+library(dplyr)
+library(ggplot2)
+
+mu_dat <- seq(log(0.02), log(0.5), length.out = 20)
+
+cert_mu_reg <- cert_fit_FE %>%
+  relative_effects(newdata = tibble::tibble(.mu = mu_dat), study = .mu) %>%
+  as_tibble() %>%
+  mutate(
+    .trt = .trtb,
+    .mu = as.numeric(as.character(.study))
+  )
+
+cert_lor <-
+  cert_net$agd_arm %>%
+  mutate(
+    probability = .r / .n,
+    log_odds = log(probability / (1 - probability))
+  ) %>%
+  group_by(.study) %>%
+  mutate(
+    .mu = log_odds[.trt == "Placebo"],
+    log_odds_ratio = log_odds - .mu,
+    sample_size = sum(n),
+  ) %>%
+  ungroup() %>%
+  filter(.trt != "Placebo")
+
+cert_lor %>%
+  ggplot(aes(x = .mu)) +
+  geom_hline(yintercept = 0, colour = "grey60") +
+  geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), data = cert_mu_reg,
+              fill = "darkred", alpha = 0.3) +
+  geom_line(aes(y = mean), data = cert_mu_reg,
+            colour = "darkred") +
+  geom_point(aes(y = log_odds_ratio, size = sample_size), alpha = 0.6) +
+  facet_wrap(vars(.trt)) +
+  labs(x = "Placebo log odds", y = "log Odds Ratio", size = "Sample Size") +
+  theme_multinma()
+```
+
 ## References
 
 ```{r certolizumab_tests, include=FALSE, eval=params$run_tests}