From 0201a90456977cc99fd616f5a9f5b2262c95b843 Mon Sep 17 00:00:00 2001
From: Kris Sankaran <sankaran.kris@gmail.com>
Date: Wed, 28 Aug 2024 22:58:02 -0500
Subject: [PATCH] add sensitivity analysis to IBD example

---
 R/sensitivity.R            | 20 ++++++++++++++++++++
 man/sensitivity_perturb.Rd | 21 +++++++++++++++++++++
 vignettes/IBD.Rmd          | 17 ++++++++++++++---
 3 files changed, 55 insertions(+), 3 deletions(-)

diff --git a/R/sensitivity.R b/R/sensitivity.R
index 27375d6..5ab98ba 100644
--- a/R/sensitivity.R
+++ b/R/sensitivity.R
@@ -393,6 +393,26 @@ covariance_matrix <- function(model, confound_ix = NULL, rho = 0.0) {
 #' @param progress A logical indicating whether to show a progress bar.
 #' @return A `date.frame` giving the outputs of `indirect_overall` across many
 #'   values of the correlation rho.
+#' @examples
+#' xy_data <- demo_spline()
+#' exper <- mediation_data(
+#'     xy_data, starts_with("outcome"), "treatment", "mediator"
+#' )
+#' model <- multimedia(
+#'     exper,
+#'     outcome_estimator = glmnet_model(lambda = 1e-2)
+#' ) |>
+#'     estimate(exper)
+#' rho_seq <- c(-0.2, 0.2)
+#' perturb <- matrix(
+#'     c(
+#'         0, 3, 0,
+#'         3, 0, 0,
+#'         0, 0, 0
+#'     ),
+#'     nrow = 3, byrow = TRUE
+#' )
+#' sensitivity_perturb(model, exper, perturb, n_bootstrap = 2)
 #' @export
 sensitivity_perturb <- function(
     model, exper, perturb, nu_seq = NULL, n_bootstrap = 100, progress = TRUE) {
diff --git a/man/sensitivity_perturb.Rd b/man/sensitivity_perturb.Rd
index 9dc4fb5..751f3b9 100644
--- a/man/sensitivity_perturb.Rd
+++ b/man/sensitivity_perturb.Rd
@@ -53,3 +53,24 @@ outcomes according to diag(sigma^2_mediator, sigma^2_outcome) + nu * perturb.
 The estimates sigma^2 are taken from the residuals in the original mediation
 and outcome models, and perturb and nu are provided by the user.
 }
+\examples{
+xy_data <- demo_spline()
+exper <- mediation_data(
+    xy_data, starts_with("outcome"), "treatment", "mediator"
+)
+model <- multimedia(
+    exper,
+    outcome_estimator = glmnet_model(lambda = 1e-2)
+) |>
+    estimate(exper)
+rho_seq <- c(-0.2, 0.2)
+perturb <- matrix(
+    c(
+        0, 3, 0,
+        3, 0, 0,
+        0, 0, 0
+    ),
+    nrow = 3, byrow = TRUE
+)
+sensitivity_perturb(model, exper, perturb, n_bootstrap = 2)
+}
diff --git a/vignettes/IBD.Rmd b/vignettes/IBD.Rmd
index 132c219..58e8cd8 100644
--- a/vignettes/IBD.Rmd
+++ b/vignettes/IBD.Rmd
@@ -276,6 +276,15 @@ ggplot(coords) +
     facet_wrap(~ type + reorder(outcome, -abundance))
 ```
 
+```{r}
+confound_ix <- data.frame(
+    outcome = vis_outcomes,
+    mediator = "gEnterocloster"
+)
+sensitivity_overall <- sensitivity(model, exper, confound_ix)
+saveRDS(sensitivity_overall, "~/Downloads/sensitivity_overall-IBD.rds")
+```
+
 Pathwise indirect effects are those that appear when modifying the treatment
 status for a single mediator. Unlike overall indirect effects, they give us
 effects for individual mediator-outcome pairs. Now, instead of averaging over
@@ -387,11 +396,13 @@ p1 | p2
 ```{r}
 confound_ix <- data.frame(
     mediator = c("gEnterocloster", "gEnterocloster", "gCAG103"),
-    outcome = c("m1129_chenodeoxycholate", "m1403_cholate", "m0900_ketodeoxycholate")
+    outcome = c(
+        "m1129_chenodeoxycholate", "m1403_cholate", "m0900_ketodeoxycholate"
+    )
 )
 
-rho_seq <- seq(-0.7, 0.7, length.out = 2)
-sensitivity_analysis <- sensitivity_pathwise(model, exper, confound_ix)#, rho_seq = rho_seq)
+sensitivity_pathwise <- sensitivity_pathwise(model, exper, confound_ix)
+saveRDS(sensitivity_pathwise, "~/Downloads/sensitivity_pathwise-IBD.rds")
 ```
 
 # Hurdle-Lognormal Approach