docs: re-worked old vignette, added two new ones

re-factor: example code in main pages only being executed if suggested packages can be loaded
RobinDenz1 · Mar 29, 2024 · 48c05c7 · 48c05c7
1 parent 94dfa02
commit 48c05c7
Show file tree

Hide file tree

Showing 51 changed files with 618 additions and 269 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -27,7 +27,8 @@ Imports:
     doRNG,
     dplyr (>= 1.0.0),
     foreach,
-    rlang
+    rlang,
+	survival (>= 3.0.0)
 Suggests: 
     MASS,
     Matching (>= 4.9),
@@ -44,7 +45,6 @@ Suggests:
     prodlim (>= 2019.11.13),
     riskRegression (>= 2020.12.08),
     rmarkdown,
-    survival (>= 3.0.0),
     testthat (>= 3.2.1),
     tidyr,
 	ggpp (>= 0.4.3),

diff --git a/NEWS.md b/NEWS.md
@@ -54,3 +54,7 @@ Documentation
 * Re-factored internal code to vastly increase speed of bootstrapping related computations
 * Added risk table functionality for `plot.adjustedsurv()`
 * Renamed `adjsurv` and `adjcif` output objects of `adjustedsurv()` and `adjsutedcif()` respectively to `adj`
+* Re-worked introduction vignette
+* Added FAQ vignette
+* Added Group Comparison vignette
+* Re-factored examples to only be executed if suggested packages are installed
diff --git a/R/adjustedcif.r b/R/adjustedcif.r
@@ -176,11 +176,11 @@ adjustedcif <- function(data, variable, ev_time, event, cause, method,
     }
 
     if (force_bounds) {
-      plotdata <- force_bounds_surv(plotdata)
+      plotdata <- force_bounds_est(plotdata)
     }
 
     if (iso_reg) {
-      plotdata <- iso_reg_surv(plotdata)
+      plotdata <- iso_reg_est(plotdata)
     }
 
     # output object
@@ -346,11 +346,11 @@ adjustedcif <- function(data, variable, ev_time, event, cause, method,
     plotdata$group <- factor(plotdata$group, levels=levs)
 
     if (force_bounds) {
-      plotdata <- force_bounds_cif(plotdata)
+      plotdata <- force_bounds_est(plotdata)
     }
 
     if (iso_reg) {
-      plotdata <- iso_reg_cif(plotdata)
+      plotdata <- iso_reg_est(plotdata)
     }
 
     out <- list(adj=plotdata,
@@ -492,9 +492,9 @@ summary.adjustedcif <- function(object, ...) {
   }
 
   if (is.null(object$mids_analyses)) {
-    cat("   - Using a single dataset")
+    cat("   - Using a single dataset\n")
   } else {
-    cat("   - Using multiply imputed dataset")
+    cat("   - Using multiply imputed dataset\n")
   }
 }
 

diff --git a/R/adjustedsurv.r b/R/adjustedsurv.r
@@ -180,18 +180,17 @@ adjustedsurv <- function(data, variable, ev_time, event, method,
     }
 
     if (force_bounds) {
-      plotdata <- force_bounds_surv(plotdata)
+      plotdata <- force_bounds_est(plotdata)
     }
 
     if (iso_reg) {
-      plotdata <- iso_reg_surv(plotdata)
+      plotdata <- iso_reg_est(plotdata)
     }
 
     # output object
     out_obj <- list(mids_analyses=out,
                     adj=plotdata,
-                    data=data$data,
-                    mids=data, # TODO: change after rework to plot()
+                    data=data,
                     method=method,
                     categorical=ifelse(length(levs)>2, TRUE, FALSE),
                     conf_level=conf_level,
@@ -358,17 +357,17 @@ adjustedsurv <- function(data, variable, ev_time, event, method,
     plotdata$group <- factor(plotdata$group, levels=levs)
 
     if (force_bounds) {
-      plotdata <- force_bounds_surv(plotdata)
+      plotdata <- force_bounds_est(plotdata)
     }
 
     if (iso_reg) {
-      plotdata <- iso_reg_surv(plotdata)
+      plotdata <- iso_reg_est(plotdata)
     }
 
     out <- list(adj=plotdata,
                 data=data,
                 method=method,
-                categorical=ifelse(length(levs) > 2, TRUE, FALSE),
+                categorical=length(levs) > 2,
                 conf_level=conf_level,
                 call=match.call())
 
@@ -531,9 +530,9 @@ summary.adjustedsurv <- function(object, ...) {
   }
 
   if (is.null(object$mids_analyses)) {
-    cat("   - Using a single dataset")
+    cat("   - Using a single dataset\n")
   } else {
-    cat("   - Using multiply imputed dataset")
+    cat("   - Using multiply imputed dataset\n")
   }
 }
 

diff --git a/R/helper_functions.r b/R/helper_functions.r
@@ -499,23 +499,29 @@ add_rows_with_zero <- function(plotdata, mode="surv") {
   return(plotdata)
 }
 
-## perform isotonic regression on survival estimates
-iso_reg_surv <- function(plotdata) {
+## perform isotonic regression on survival / CIF estimates
+iso_reg_est <- function(plotdata) {
 
-  if (anyNA(plotdata$surv)) {
+  mode <- ifelse("surv" %in% colnames(plotdata), "surv", "cif")
+
+  if (anyNA(plotdata[, mode])) {
     stop("Isotonic Regression cannot be used when there are missing",
-         " values in the final survival estimates.")
+         " values in the final estimates.")
   }
 
   for (lev in levels(plotdata$group)) {
     temp <- plotdata[plotdata$group==lev, ]
-    # to surv estimates
-    new <- rev(stats::isoreg(rev(temp$surv))$yf)
-    plotdata$surv[plotdata$group==lev] <- new
+
+    if (mode=="surv") {
+      new <- rev(stats::isoreg(rev(temp$surv))$yf)
+    } else {
+      new <- stats::isoreg(temp$cif)$yf
+    }
+    plotdata[, mode][plotdata$group==lev] <- new
 
     # shift confidence intervals accordingly
     if ("ci_lower" %in% colnames(temp)) {
-      diff <- temp$surv - new
+      diff <- temp[, mode] - new
 
       plotdata$ci_lower[plotdata$group==lev] <- temp$ci_lower - diff
       plotdata$ci_upper[plotdata$group==lev] <- temp$ci_upper - diff
@@ -526,50 +532,24 @@ iso_reg_surv <- function(plotdata) {
 }
 
 ## force probabilities to be in the 0/1 range
-force_bounds_surv <- function(plotdata) {
-  plotdata <- within(plotdata, {
-    surv <- ifelse(surv < 0, 0, surv)
-    surv <- ifelse(surv > 1, 1, surv)
-  })
-
-  return(plotdata)
-}
-
-## isotonic regression for CIF estimates
-iso_reg_cif <- function(plotdata) {
-
-  if (anyNA(plotdata$cif)) {
-    stop("Isotonic Regression cannot be used when there are missing",
-         " values in the final CIF estimates.")
-  }
-
-  for (lev in levels(plotdata$group)) {
-    temp <- plotdata[plotdata$group==lev, ]
-
-    new <- stats::isoreg(temp$cif)$yf
-    plotdata$cif[plotdata$group==lev] <- new
-
-    if ("ci_lower" %in% colnames(temp)) {
-      diff <- temp$cif - new
-
-      plotdata$ci_lower[plotdata$group==lev] <- temp$ci_lower - diff
-      plotdata$ci_upper[plotdata$group==lev] <- temp$ci_upper - diff
-    }
+force_bounds_est <- function(plotdata) {
+  mode <- ifelse("surv" %in% colnames(plotdata), "surv", "cif")
+
+  if (mode=="surv") {
+    plotdata <- within(plotdata, {
+      surv <- ifelse(surv < 0, 0, surv)
+      surv <- ifelse(surv > 1, 1, surv)
+    })
+  } else {
+    plotdata <- within(plotdata, {
+      cif <- ifelse(cif < 0, 0, cif)
+      cif <- ifelse(cif > 1, 1, cif)
+    })
   }
 
   return(plotdata)
 }
 
-## forcing bounds for CIF estimates
-force_bounds_cif <- function(plotdata) {
-  plotdata <- within(plotdata, {
-    cif <- ifelse(cif < 0, 0, cif)
-    cif <- ifelse(cif > 1, 1, cif)
-  })
-
-  return(plotdata)
-}
-
 ## Given a mids object and our column names of interest, calculate the
 ## maximum observed (cause-specific) event time
 max_observed_time <- function(mids, variable, ev_time, event, levs, cause,

diff --git a/R/input_checks.r b/R/input_checks.r
@@ -1229,7 +1229,7 @@ check_inputs_auc_diff <- function(times, max_t, color, difference, ratio,
 }
 
 ## input checks when using risk tavles in plot.adjustedsurv()
-check_inputs_risk_table <- function(method, type, use_weights, warn) {
+check_inputs_risk_table <- function(method, type, use_weights, stratify, warn) {
 
   # errors
   if (!type %in% c("n_at_risk", "n_cens", "n_events")) {
@@ -1238,8 +1238,8 @@ check_inputs_risk_table <- function(method, type, use_weights, warn) {
   }
 
   # (optional) warnings
-  if (!method %in% c("km", "iptw_km", "iptw_cox") && warn) {
-    warning("Adding risk tables may produce confusing output when",
+  if (!method %in% c("km", "iptw_km", "iptw_cox") && stratify && warn) {
+    warning("Adding stratified risk tables may produce confusing output when",
             " using methods other then 'km', 'iptw_km' or 'iptw_cox',",
             " because all other methods do not use risk tables to estimate",
             " the survival curves. See details. Set risk_table_warn=FALSE",
@@ -1252,4 +1252,3 @@ check_inputs_risk_table <- function(method, type, use_weights, warn) {
             " details. Set risk_table_warn=FALSE to silence this warning.")
   }
 }
-
diff --git a/R/method_direct.r b/R/method_direct.r
@@ -118,11 +118,11 @@ surv_g_comp <- function(outcome_model, data, variable, times,
                                          "wglm", "hal9001"))) {
       requireNamespace("riskRegression")
 
-      surv_lev <- riskRegression::predictRisk(object=outcome_model,
+      surv_lev <- quiet(riskRegression::predictRisk(object=outcome_model,
                                               newdata=data_temp,
                                               times=times,
                                               cause=1,
-                                              ...)
+                                              ...))
       surv_lev <- 1 - surv_lev
 
       # using predictProb
@@ -255,11 +255,11 @@ cif_g_comp <- function(outcome_model, data, variable, times,
                                          "riskRegression", "ARR",
                                          "hal9001"))) {
 
-      surv_lev <- riskRegression::predictRisk(object=outcome_model,
+      surv_lev <- quiet(riskRegression::predictRisk(object=outcome_model,
                                               newdata=data_temp,
                                               times=times,
                                               cause=cause,
-                                              ...)
+                                              ...))
 
       # for fastCrr in fastcmprsk
     } else if (inherits(outcome_model, "fcrr")) {

diff --git a/R/plot.adjustedcif.r b/R/plot.adjustedcif.r
@@ -59,12 +59,12 @@ plot.adjustedcif <- function(x, conf_int=FALSE, max_t=Inf,
   # in some methods estimates can be outside the 0, 1 bounds,
   # if specified set those to 0 or 1 respectively
   if (force_bounds) {
-    plotdata <- force_bounds_cif(plotdata)
+    plotdata <- force_bounds_est(plotdata)
   }
 
   # apply isotonic regression if specified
   if (iso_reg) {
-    plotdata <- iso_reg_cif(plotdata)
+    plotdata <- iso_reg_est(plotdata)
   }
 
   mapping <- ggplot2::aes(x=.data$time, y=.data$cif, color=.data$group,

diff --git a/R/plot.adjustedsurv.r b/R/plot.adjustedsurv.r
@@ -84,12 +84,12 @@ plot.adjustedsurv <- function(x, conf_int=FALSE, max_t=Inf,
   # in some methods estimates can be outside the 0, 1 bounds,
   # if specified set those to 0 or 1 respectively
   if (force_bounds) {
-    plotdata <- force_bounds_surv(plotdata)
+    plotdata <- force_bounds_est(plotdata)
   }
 
   # apply isotonic regression if specified
   if (iso_reg) {
-    plotdata <- iso_reg_surv(plotdata)
+    plotdata <- iso_reg_est(plotdata)
   }
 
   # plot CIF instead of survival
@@ -231,6 +231,7 @@ plot.adjustedsurv <- function(x, conf_int=FALSE, max_t=Inf,
 
     check_inputs_risk_table(method=x$method, type=risk_table_type,
                             use_weights=risk_table_use_weights,
+                            stratify=risk_table_stratify,
                             warn=risk_table_warn)
 
     # set correct weights if specified
@@ -251,15 +252,8 @@ plot.adjustedsurv <- function(x, conf_int=FALSE, max_t=Inf,
       variable <- NULL
     }
 
-    # set correct data
-    if (!is.null(x$mids_analyses)) {
-      data <- x$mids
-    } else {
-      data <- x$data
-    }
-
     p <- add_risk_table(p_surv=p,
-                        data=data,
+                        data=x$data,
                         event=x$call$event,
                         ev_time=x$call$ev_time,
                         variable=variable,
@@ -395,16 +389,22 @@ get_censoring_ind_data <- function(x, steps, max_t, plotdata) {
 
   levs <- levels(x$adj$group)
 
+  if (is.null(x$mids_analyses)) {
+    data <- x$data
+  } else {
+    data <- x$data$data
+  }
+
   # keep only relevant data
-  x$data <- x$data[which(x$data[, x$call$ev_time] <= max_t), ]
+  data <- data[which(data[, x$call$ev_time] <= max_t), ]
 
   # create needed data.frame
   cens_dat <- vector(mode="list", length=length(levs))
   for (i in seq_len(length(levs))) {
 
     # times with censoring
-    cens_times <- sort(unique(x$data[, x$call$ev_time][
-      x$data[, x$call$event]==0 & x$data[, x$call$variable]==levs[i]]))
+    cens_times <- sort(unique(data[, x$call$ev_time][
+      data[, x$call$event]==0 & data[, x$call$variable]==levs[i]]))
     # y axis place to put them
     adjsurv_temp <- plotdata[plotdata$group==levs[i] & !is.na(plotdata$surv), ]
 

diff --git a/man/CSC_MI.Rd b/man/CSC_MI.Rd
@@ -45,6 +45,8 @@ Robin Denz
 \donttest{
 library(adjustedCurves)
 library(survival)
+
+if (requireNamespace("riskRegression") & requireNamespace("mice")) {
 library(riskRegression)
 library(mice)
 
@@ -56,11 +58,12 @@ sim_dat$group <- as.factor(sim_dat$group)
 sim_dat$x1 <- ifelse(runif(n=50) < 0.5, sim_dat$x1, NA)
 
 # perform multiple imputation
-mids <- mice::mice(data=sim_dat, method="pmm", m=5)
+mids <- mice::mice(data=sim_dat, method="pmm", m=5, printFlag=0)
 
 # use the function
 csc_mods <- CSC_MI(mids=mids,
                    formula=Hist(time, event) ~ x1 + x2 + x3 + x4 + x5 + x6 + group
                    )
 }
 }
+}
diff --git a/man/FGR_MI.Rd b/man/FGR_MI.Rd
@@ -48,6 +48,10 @@ Robin Denz
 \donttest{
 library(adjustedCurves)
 library(survival)
+
+if (requireNamespace("riskRegression") & requireNamespace("prodlim") &
+    requireNamespace("mice")) {
+
 library(riskRegression)
 library(mice)
 library(prodlim)
@@ -68,3 +72,4 @@ fgr_mods <- FGR_MI(mids=mids,
                    cause=1)
 }
 }
+}