Some restructuring so that the simulations can be run without access to the raw data

Author: dincerti

01-setup.R

@@ -0,0 +1,63 @@
+# R Packages -------------------------------------------------------------------
+library("data.table")
+library("doParallel")
+library("doRNG") # Reproducible parallel foreach loops
+library("dplyr") 
+library("foreach")
+library("glmnet")
+library("ggplot2")
+library("ggpubr")
+library("impute") # From bioconductor
+library("knitr")
+library("parallel")
+library("purrr")
+library("rngtools") # Required by dorng
+library("rsample")
+library("survival")
+library("tidyr") # (>= 1.0.0) for pivot_longer() and pivot_wider()
+# library("xfun") # We will use xfun::cache_rds() below, but not attach the package here
+library("xtable")
+
+# Custom R functions -----------------------------------------------------------
+source("R/plot_patient_followup.R")
+source("R/impute_genes.R")
+source("R/make_xy.R")
+source("R/calibrate_sim.R")
+source("R/run_sim.R")
+source("R/adjust_surv.R")
+source("R/tidycoef.R")
+source("R/concordance.R")
+source("R/calibrate.R")
+
+# Settings ---------------------------------------------------------------------
+set.seed(77)
+theme_set(theme_bw())
+center_title <- function() theme(plot.title = element_text(hjust = 0.5))
+N_SIMS <- 200
+
+# Caching
+RERUN_CACHE <- TRUE # Set to TRUE to rerun all cached results
+if (!dir.exists("cache")){
+  dir.create("cache")
+}
+
+# Parallel
+PARALLEL <- TRUE
+if (PARALLEL) {
+  cl <- parallel::makeCluster(20, setup_strategy = "sequential", outfile = "simout")
+  registerDoParallel(cl)
+}
+
+# Model formula ----------------------------------------------------------------
+# This formula is used for both the simulations and the real-world data application
+f_small <- formula(
+  ~ PracticeType + index_date_year + Race + age_at_dx + RE_TP53 +
+    CN_TP53 + SV_TP53 + SV_KRAS + SV_EGFR
+)
+vars_small <- c("Practice type: community", "Index year", "African American",
+                "Other race", "White", "Age", "TP53 RE",
+                "TP53 CN", "TP53 SV", "KRAS SV", "EGFR SV")
+
+# Followup plot ----------------------------------------------------------------
+p_followup <- plot_patient_followup()
+ggsave("figs/followup.pdf", p_followup, height = 5, width = 7)

02-simulation.R

@@ -0,0 +1,89 @@
+# Calibrate simulation ---------------------------------------------------------
+if (file.exists("sim_settings.rds")) { # Interested readers can run this
+  sim_settings <- readRDS("sim_settings.rds")
+} else { # This setting is run for the paper
+  data <- readRDS("data.rds")
+  sim_settings <- calibrate_sim(f_small, data = data, store_x = TRUE)
+}
+
+# Save cumulative hazard plot
+ggsave("figs/sim_calibration_cumhaz.pdf", 
+       sim_settings$cumhaz_plot,
+       height = 5, width = 10)
+
+# Example simulated data -------------------------------------------------------
+# No predictors (intercept only model)
+params <- set_params(sim_settings = sim_settings, dist = "weibullPH")
+simdata <- sim_survdata(params = params, n_pats = 10000)
+simdata_summary_int <- summarize_simdata(simdata, km_rwd = sim_settings$os_km$rc,  
+                                         save = TRUE, name = "int")
+
+# Now include predictors (and hence create informative censoring when
+# using the Kaplan-Meier estimator). The design matrix X and parameters will 
+# be used in the simulations that follow
+x_sim <- sim_x(n_pats = 5000, sim_settings, p_bin = 0)
+params <- set_params(x_sim, sim_settings, dist = "weibullPH")
+simdata <- sim_survdata(x_sim, params)
+simdata_summary_p11 <- summarize_simdata(simdata, km_rwd = sim_settings$os_km$rc,
+                                         save = TRUE, name = "p11")
+
+# Run simulation for unpenalized Cox model -------------------------------------
+sim_coxph_p21 <- xfun::cache_rds({
+  run_sim(n_sims = N_SIMS, x = x_sim, params = params, method = "coxph")
+}, file = "sim_coxph_p21.rds", rerun = RERUN_CACHE)
+sim_coxph_21_summary <- summarize_sim(sim_coxph_p21, save = TRUE,
+                                      model_name = "coxph_p21")
+rm(sim_coxph_p21)
+
+# Run simulation for lasso model with lambda = 0 and small p -------------------
+# run_sim1(simdata, method = "coxnet", lambda = 0) # For debugging
+sim_coxlasso_lambda0_p21 <- xfun::cache_rds({
+  run_sim(n_sims = N_SIMS, x = x_sim, params = params, 
+          method = "coxnet", lambda = c(1, 0))
+}, file = "sim_coxlasso_lambda0_p21.rds", rerun = RERUN_CACHE)
+coxlasso_lambda0_p21_summary <- summarize_sim(
+  sim_coxlasso_lambda0_p21, save = TRUE, 
+  model_name = "coxlasso_lambda0_p21"
+)
+rm(sim_coxlasso_lambda0_p21)
+
+# Run simulation for lasso model with small p ----------------------------------
+sim_coxlasso_p21 <- xfun::cache_rds({
+  run_sim(n_sims = N_SIMS, x = x_sim, params = params, 
+          method = "coxnet")
+}, file = "sim_coxlasso_p21.rds", rerun = RERUN_CACHE)
+sim_coxlasso_p21_summary <- summarize_sim(sim_coxlasso_p21, save = TRUE, 
+                                          model_name = "coxlasso_p21")
+rm(sim_coxlasso_p21)
+
+# Run simulation for ridge model with small p ----------------------------------
+sim_coxridge_p21 <- xfun::cache_rds({
+  run_sim(n_sims = N_SIMS, x = x_sim, params = params, 
+          method = "coxnet", alpha = 0)
+}, file = "sim_coxridge_p21.rds", rerun = RERUN_CACHE)
+sim_coxridge_p21_summary <- summarize_sim(sim_coxridge_p21, save = TRUE, 
+                                          model_name = "coxridge_p21")
+rm(sim_coxridge_p21)
+
+# Run simulation for lasso model with large p ----------------------------------
+x_sim <- sim_x(n_pats = 5000, sim_settings, p_bin = 1000)
+params <- set_params(x_sim, sim_settings, dist = "weibullPH")
+
+sim_coxlasso_p1011 <- xfun::cache_rds({
+  run_sim(n_sims = N_SIMS, x = x_sim, params = params, 
+          method = "coxnet")
+}, file = "sim_coxlasso_p1011.rds", rerun = RERUN_CACHE)
+sim_coxlasso_p1011_summary <- summarize_sim(sim_coxlasso_p1011, save = TRUE, 
+                                            model_name = "coxlasso_p1011")
+rm(sim_coxlasso_p1011)
+
+# Calibration plot comparing small and large p simulation ----------------------
+sim_coxlasso_calplot <- ggarrange(
+  sim_coxlasso_p21_summary$calibration_plots$complete +
+    ggtitle("(A) Small p") + center_title(),
+  sim_coxlasso_p1011_summary$calibration_plots$complete +
+    ggtitle("(B) Large p") + center_title(),
+  nrow = 2, common.legend = TRUE, legend = "bottom"
+)
+ggsave("figs/sim_coxlasso_calibration_complete.pdf", sim_coxlasso_calplot, 
+       width = 7, height = 9)