Simulation Experiment Recipe
Objectives
NOTE: The experiments presented herein are a subset ofthose described in “A FlexibleApproach for Predictive Biomarker Discovery” (P. Boileau, et al.2022). This subset was adapted by theoriginal author for a case study of the R packagesimChef. To reproduce this case studyand documentation, please see the repositorysimChef-case-study.
Introduction
An endeavor central to precision medicine is predictive biomarker discovery;they define patient sub-populations which stand to benefit most, or least, froma given treatment. The identification of these biomarkers is often the byproductof the related but fundamentally different task of treatment rule estimation.When the number of potentially predictive biomarkers is commensurate with, ormuch larger than, the number of observations, the set of biomarkers designatedas predictive by these methods will generally contain many false positives.
In synthetic and real-world data-inspired simulation settings emulatingrandomized control trials with high-dimensional covariate vectors, we study aproposed predictive biomarker discovery method, uniCATE. The objectives of thissimulation study are:
Evaluate the proposed method as an estimator of a parameter for directly assessing the importance of potentially predictive biomarkers.
Examine the empirical behavior, like type-I error rate control, of the proposed method as the sample size \(n\) grows to determine its correspondence to asymptotic results.
Study the statistical robustness of the proposed method under a number of misspecified nuisance parameter estimator scenarios.
This work is motivated by the need to identify predictive biomarkers inclinical trials. Of particular interest is their detection for drug targetdiscovery and diagnostic assay development. The former requires theidentification of biomarkers causally related to the outcome of interest,whereas the latter seeks a set of strongly predictive biomarkers. Thissimulation study is therefore representative of these applications.
Notation
Consider \(n\) identically and independently distributed random vectors \(X_i =(W_i, A_i, Y^{(1)}_i, Y^{(0)}_i) \sim P_X\), \(i = 1, \ldots, n\), correspondingto complete but unobserved data generated by participants in an idealizedrandomized control trial or observational study, where:
- \(V\): pre-treatment covariates, such as location and income, of dimension \(p\).
- \(B\): pre-treatment biomarkers, such as gene expression data, of dimension \(q\).
- \(W = (V, B)\): a \((q+p)\)-length random vector.
- \(A\): a binary random variable representing a treatment assignment.
- \(Y^{(1)}\) and \(Y^{(0)}\): partially observed random variables corresponding to the potential outcomes of clinical interest under both treatment and control conditions, respectively, for each patient.
- \(P_X\): The unknown data-generating distribution of the full data \((W, A, Y^{(0)}, Y^{(1)})\).
- \(Y = AY^{(1)} + (1-A)Y^{(0)}\): observed outcomes where \(Y_i = Y_i(1)\) if patient \(i\) received treatment, and \(Y_i = Y_i(0)\) otherwise.
- \(P_0\): The unknown data-generating distribution of the observed data \((W, A, Y)\).
- \(P_n\): The observed empirical distribution of a sample of size \(n\) drawn from \(P_0\).
Clinically predictive biomarker importance parameter
Clinically relevant predictive biomarkers are often those that have a stronginfluence on the outcome of interest on the absolute scale. As such, an idealtarget of inference when these outcomes are continuous and the number ofcovariates small is the CATE conditioning on the set of biomarkers:
\[\begin{equation*} \mathbb{E}_{P_X}\left[ Y^{(1)} - Y^{(0)} \big| B = b \right].\end{equation*}\]
Accurate and interpretable estimation of this parameter is generally challengingwhen \(p\) is large, preventing the accurate recovery of predictive biomarkers.
Indexing the biomarkers of by \(j = 1, \ldots, p\), such that \(B = (B_{1},\ldots, B_{p})\), centering them such that \(\mathbb{E}_{P_X}[B_{j}] = 0\), andassuming that \(\mathbb{E}_{P_X}[B^2_{j}] > 0\), we instead target the full-datavariable importance parameter \(\Psi^F(P_X) = (\Psi^F_1(P_X), \ldots,\Psi_p^F(P_X))\) where
\[\begin{equation*} \Psi_{j}^F(P_X) \equiv \frac{\mathbb{E}_{P_X}\left[\left(Y^{(1)} - Y^{(0)}\right)B_{j}\right]} {\mathbb{E}_{P_X}\left[B_{j}^2\right]}.\end{equation*}\]
Under the assumption that the expected difference in potential outcomes admits alinear form when conditioning on any given \(B_j\), \(\Psi^F(P_X)\) is the vectorof expected simple linear regression coefficients produced by regressing thedifference in potential outcomes against each biomarker. While the truerelationship between the difference of potential outcomes and a predictivebiomarker is almost surely nonlinear, \(\Psi^F(P_X)\) is a generally informativetarget of inference. Biomarkers with the largest absolute values in\(\Psi^F(P_X)\) generally modify the effect of treatment the most.
Unfortunately, we generally don’t observe both \(Y^{(0)}\) and \(Y^{(1)}\) for anygiven random unit. We only observe \(Y\). \(P_X\) is censored by a treatmentassignment mechanism; we can only interrogate parameters of the observeddata-generating distribution \(P_0\). Luckily, with assumptions of no unmeasuredconfounding and overlapping treatment support in which observations may beassigned to either treatment condition regardless of covariates, we have\[\begin{equation} \begin{split} \Psi_j(P_0) & \equiv \frac{\mathbb{E}_{P_0}\left[\left(\bar{Q}_0(1, W) - \bar{Q}_0(0, W)\right)B_{j}\right]} {\mathbb{E}_{P_0}\left[B_{j}^2\right]} \\ & = \Psi_j^F(P_X) \end{split}\end{equation}\]where \(\bar{Q}_0(a, w) \equiv \mathbb{E}_{P_0}[Y|A = a, W = w]\) is theconditional expected outcome given treatment and covariates.
It follows that \(\Psi(P_0) = \Psi^F(P_X)\) in a randomized control trial.\(\Psi_j(P_0)\) therefore measures the full data predictive biomarker importance.We propose a method, uniCATE, to perform inference about this parameter.
Data Generation
LM with TEM
“LM with TEM” stands for “Linear Model with Treatment Effect Modification”.Its generative model is as follows:\[\begin{equation} \begin{split} W = B & \sim N(0, I_{100 \times 100}) \\ A | W = A & \sim \text{Bernoulli}(1/2) \\ Y | A, W & \sim N\left(W^\top \left(\beta + I(A = 1)\gamma^{(1)} + I(A = 0)\gamma^{(0)}\right), 1/2\right). \end{split}\end{equation}\]Here, \(\beta = (\beta_1, \ldots, \beta_{100})^\top\) such that \(\beta_1 = \ldots= \beta_{20} = 2\) and \(\beta_{21} = \ldots = \beta_{100} = 0\), and\(\gamma^{(a)} = (\gamma_1^{(a)}, \ldots, \gamma_{100}^{(a)})^\top\) where\(\gamma_1^{(1)} = \ldots = \gamma_{50}^{(1)} = 5\), \(\gamma_1^{(0)} = \ldots =\gamma_{50}^{(0)} = 5\) and \(\gamma_{51}^{(1)} = \ldots = \gamma_{100}^{(1)} =0\) for \(a = \{0, 1\}\).
Two hundred datasets each of 125, 250 and 500 observations are generated fromthis model.
Function
## function (n = 125, cov_mat = diag(1, nrow = 100))
## {
## beta_main <- c(rep(2, 20), rep(0, 80))
## beta_0 <- rep(0, 100)
## beta_1 <- c(rep(5, 50), rep(0, 50))
## W <- MASS::mvrnorm(n = n, mu = rep(0, 100), Sigma = cov_mat)
## colnames(W) <- paste0("W", seq_len(100))
## A <- ifelse(runif(n) < 0.5, 0, 1)
## epsilon_0 <- rnorm(n = n, mean = 0, sd = 0.5)
## epsilon_1 <- rnorm(n = n, mean = 0, sd = 0.5)
## W_t <- t(W)
## main_effects <- crossprod(W_t, beta_main)
## Y_0 <- as.vector(main_effects + crossprod(W_t, beta_0) +
## epsilon_0)
## Y_1 <- as.vector(main_effects + crossprod(W_t, beta_1) +
## epsilon_1)
## Y <- ifelse(A == 0, Y_0, Y_1)
## sample_list <- list(Y = Y, A = A, W = W)
## return(sample_list)
## }Input Parameters
## list()Kinked with TEM
“Kinked with TEM” stands for “Kinked Model with Treatment Effect Modification”.Its generative model is as follows:\[\begin{equation} \begin{split} W = B & \sim N(0, I_{100 \times 100}) \\ A | W = A & \sim \text{Bernoulli}(1/2) \\ Y | A, W & \sim N\left( W^\top \left(I(A = 1)\gamma + I(A = 0)\;\text{diag}(I(W > 0))\;\gamma\right), 1/2 \right), \end{split}\end{equation}\]where \(\gamma = (\gamma_1, \ldots, \gamma_{100})\), \(\gamma_{1} = \ldots= \gamma_{50} = 10\), \(\gamma_{51} = \ldots = \gamma_{100} = 0\), and\(\text{diag}(\cdot)\) is a diagonal matrix whose diagonal equals theinput vector.
As with the LM with TEM mode, 200 datasets of 125, 250 and 500observations are generated.
Function
## function (n = 125, cov_mat = diag(1, nrow = 100))
## {
## beta_tem <- c(rep(10, 50), rep(0, 50))
## W <- MASS::mvrnorm(n = n, mu = rep(0, 100), Sigma = cov_mat)
## colnames(W) <- paste0("W", seq_len(100))
## A <- ifelse(runif(n) < 0.5, 0, 1)
## epsilon_1 <- rnorm(n = n, mean = 0, sd = 0.5)
## W_t <- t(W)
## Y_1 <- as.vector(crossprod(W_t, beta_tem) + epsilon_1)
## epsilon_0 <- rnorm(n = n, mean = 0, sd = 0.5)
## W_mod <- W[, 1:50]
## W_mod[W_mod < 0] <- 0
## W_mod_t <- t(cbind(W_mod, W[, 51:100]))
## Y_0 <- as.vector(crossprod(W_mod_t, beta_tem) + epsilon_0)
## Y <- ifelse(A == 0, Y_0, Y_1)
## sample_list <- list(Y = Y, A = A, W = W)
## return(sample_list)
## }Input Parameters
## list()NLM with TEM
“NLM with TEM” stands for “Nonlinear Model with Treatment Effect Modification”.Its generative model is as follows:\[\begin{equation} \begin{split} W = B & \sim N(0, I_{100 \times 100}) \\ A | W = A & \sim \text{Bernoulli}(1/2) \\ Y | A, W & \sim N\left( W^\top \left(I(A = 1)\gamma + I(A = 0)\;\text{diag}(I(W > 0))\;\gamma\right), 1/2 \right), \end{split}\end{equation}\]where \(\beta_1 = \ldots = \beta_{20} = 1\) and \(\beta_{21} = \ldots =\beta_{100} = 0\), and where \(\gamma_1 = \ldots = \gamma_{50} = 5\) and\(\gamma_{51} = \ldots = \gamma_{100} = 0\).
Again, 200 datasets are generated for each sample size: 125, 250 and 500.
Function
## function (n = 125, cov_mat = diag(1, nrow = 100))
## {
## beta_main <- c(rep(1, 20), rep(0, 80))
## beta_0 <- rep(0, 100)
## beta_1 <- c(rep(5, 50), rep(0, 50))
## W <- MASS::mvrnorm(n = n, mu = rep(0, 100), Sigma = cov_mat)
## colnames(W) <- paste0("W", seq_len(100))
## A <- ifelse(runif(n) < 0.5, 0, 1)
## epsilon_0 <- rnorm(n = n, mean = 0, sd = 0.5)
## epsilon_1 <- rnorm(n = n, mean = 0, sd = 0.5)
## W_t <- t(W)
## main_effects <- rowSums(exp(abs(t((W_t * beta_main)))))
## Y_0 <- as.vector(main_effects + crossprod(W_t, beta_0) +
## epsilon_0)
## Y_1 <- as.vector(main_effects + crossprod(W_t, beta_1) +
## epsilon_1)
## Y <- ifelse(A == 0, Y_0, Y_1)
## sample_list <- list(Y = Y, A = A, W = W)
## return(sample_list)
## }Input Parameters
## list()Methods and Evaluation
Methods
uniCATE (LASSO)
We propose a cross-validated estimator that uses all available data. Begin byrandomly partitioning the \(n\) observations of \(P_n\) into \(K\) independentvalidation sets \(P_{n, 1}^1, \ldots, P_{n, K}^1\) of approximately equal size.For \(k=1,\ldots, K\), define the training set as, in a slight abuse of notation,\(P_{n, k}^0 = P_n \setminus P_{n,k}^1\). Then the cross-validated estimator of\(\Psi_j(P_0)\) is defined as:
\[\begin{equation} \Psi_j^{(\text{CV})}(P_n) = \frac{1}{K} \sum_{k=1}^K \frac{\sum_{i=1}^n I(O_i \in P_{n, k}^1)\tilde{T}(O_i; P_{n,k}^0)B_{ij}} {\sum_{i=1}^n I(O_i \in P_{n, k}^1) B_{ij}^2},\end{equation}\]where \(\tilde{T}(O;P)\) is the difference in augmented inverse probabilityweight transformed outcomes. That is,\[\begin{equation} \tilde{T}(O; P) = \left(\frac{I(A=1)}{g(1, W)} - \frac{I(A=0)}{g(0,W)}\right) (Y - \bar{Q}(A, W)) + \bar{Q}(1, W) - \bar{Q}(0,W),\end{equation}\]where \(g(a,W) = \mathbb{P}_{P}[A=a|W]\). Here, \(P\) is omitted from the subscriptof \(g(a,W),\bar{Q}(A,W)\) to simplify notation.
This cross-validated estimator estimates \(\bar{Q}_0(A,W)\) on the training setusing a LASSO regression that includes main and treatment-biomarker interactionterms. The simple linear regression coefficient of the \(j\text{th}\) biomarkerregressed on the difference in potential outcomes is then computed on thevalidation set using the estimate of \(\bar{Q}_0(A,W)\). This procedure isrepeated \(K\) times, and the estimate is defined as the mean of the estimatedslopes for each validation set.
In a randomized control trial with known treatment assignment mechanism, thisestimator is asymptotically normal about the true parameter with variance givenby its efficient influence curve. Hypothesis testing is therefore possiblethrough the use of Wald-type confidence intervals. We define the nullhypothesis as \(\Psi(P_0)=0\), and use an FDR cutoff of 5% throughout thesimulation study. For more information, see (P. Boileau, et al.2022).
Function
## function (Y, A, W, use_sl = FALSE)
## {
## sample_tbl <- tibble::tibble(Y = Y, A = A)
## sample_tbl <- dplyr::bind_cols(sample_tbl, W)
## biomarker_names <- colnames(W)
## propensity_score_ls <- list(`1` = 0.5, `0` = 0.5)
## if (use_sl) {
## interactions <- lapply(biomarker_names, function(b) c("A",
## b))
## lrnr_interactions <- sl3::Lrnr_define_interactions$new(interactions)
## lrnr_lasso <- sl3::make_learner(sl3::Pipeline, lrnr_interactions,
## sl3::Lrnr_glmnet$new())
## lrnr_enet <- sl3::make_learner(sl3::Pipeline, lrnr_interactions,
## sl3::Lrnr_glmnet$new(alpha = 0.5))
## lrnr_ridge <- sl3::make_learner(sl3::Pipeline, lrnr_interactions,
## sl3::Lrnr_glmnet$new(alpha = 0))
## lrnr_spline <- sl3::make_learner(sl3::Pipeline, lrnr_interactions,
## sl3::Lrnr_polspline$new())
## lrnr_rf <- sl3::make_learner(sl3::Pipeline, lrnr_interactions,
## sl3::Lrnr_ranger$new())
## lrnr_mean <- sl3::Lrnr_mean$new()
## learner_library <- sl3::make_learner(sl3::Stack, lrnr_spline,
## lrnr_lasso, lrnr_enet, lrnr_ridge, lrnr_rf, lrnr_mean)
## super_learner <- sl3::Lrnr_sl$new(learners = learner_library,
## metalearner = sl3::make_learner(sl3::Lrnr_nnls))
## }
## else {
## super_learner <- NULL
## }
## unicate_tbl <- uniCATE::unicate(data = sample_tbl, outcome = "Y",
## treatment = "A", covariates = biomarker_names, biomarkers = biomarker_names,
## propensity_score_ls = propensity_score_ls, v_folds = 5L,
## super_learner = super_learner)
## tems <- unicate_tbl %>% dplyr::filter(p_value_bh <= 0.05) %>%
## dplyr::pull(biomarker)
## tems <- list(tems = tems)
## return(tems)
## }Input Parameters
## $use_sl
## [1] FALSEModified Covariates
uniCATE’s capacity to identify predictive biomarkers was compared to that ofpopular CATE estimation methods: the modified covariates approach and itsaugmented counterpart of Tian et al.(2012).Briefly, the former directly estimates the linear model coefficients of thetreatment-biomarker interactions, using a linear working model for these terms,without having to model or estimate the main effects. While Tian et al.’smethod is flexible since it avoids making any assumptions about the functionalform of the main biomarker effects, it can lack finite sample precision insmall-sample, high-dimensional settings. When using a penalized method like theLASSO, biomarkers are designated as predictive when their estimatedcoefficients are non-zero.
Function
## function (Y, A, W)
## {
## propensity_func <- function(x, trt) 0.5
## mod_cov <- personalized::fit.subgroup(x = W, y = Y, trt = A,
## propensity.func = propensity_func, loss = "sq_loss_lasso",
## nfolds = 10)
## coefs <- mod_cov$coefficients[paste0("W", seq_len(100)),
## ]
## tems <- names(coefs[which(coefs != 0)])
## tems <- list(tems = tems)
## return(tems)
## }Input Parameters
## list()Augmented Modified Covariates
Tian et al. (2012) proposed an “augmented” version of the modified covariatesmethod that explicitly account for the variation due to the main effects of thebiomarkers. When the outcome is continuous, it is equivalent to fitting a(penalized) multivariate linear regression with biomarker andtreatment-biomarker interaction terms. Again, biomarkers are classified aspredictive when their estimated treatment-biomarker coefficients are non-zero.For more information on these methods, please see “A Simple Method forEstimating Interactions Between a Treatment and a Large Number of Covariates”(L. Tian, et al.2012).
Function
## function (Y, A, W)
## {
## propensity_func <- function(x, trt) 0.5
## augment_func <- function(x, y) {
## df.x <- data.frame(x)
## form <- eval(paste(" ~ 1 + ", paste(colnames(df.x), collapse = " + ")))
## mm <- model.matrix(as.formula(form), data = df.x)
## cvmod <- glmnet::cv.glmnet(y = y, x = mm, nfolds = 10)
## predictions <- predict(cvmod, newx = mm, s = "lambda.min")
## predictions
## }
## aug_mod_cov <- personalized::fit.subgroup(x = W, y = Y, trt = A,
## propensity.func = propensity_func, loss = "sq_loss_lasso",
## augment.func = augment_func, nfolds = 10)
## coefs <- aug_mod_cov$coefficients[paste0("W", seq_len(100)),
## ]
## tems <- names(coefs[which(coefs != 0)])
## tems <- list(tems = tems)
## return(tems)
## }Input Parameters
## list()Evaluation
Empirical FDR
The empirical false discovery rate (FDR) for each method is computed as themean of the false discovery proportions (FDP) stratified over eachdata-generating process and sample size. For any given simulated dataset, theFDP is computed as follows:\[\begin{equation} \text{FDP} = \frac{\text{# of biomarkers erroneously classified as predictive}} {\text{# of biomarkers classified as predictive}}.\end{equation}\]When no biomarkers are classified as predictive, \(\text{FDP} = 0\).
Function
## function (fit_results)
## {
## group_vars <- c(".dgp_name", ".method_name", "n")
## eval_out <- fit_results %>% mutate(fdp = purrr::map_dbl(tems,
## fdp_fun)) %>% dplyr::group_by(dplyr::across({
## {
## group_vars
## }
## })) %>% summarize(fdr = mean(fdp), .groups = "drop")
## return(eval_out)
## }Input Parameters
## list()Empirical TPR
The empirical true positive rate (TPR) is the mean of the true positiveproportions (TPP) stratified over each data-generating process and sample size.For any given simulated dataset, the TPP is computed as\[\begin{equation} \text{TPP} = \frac{\text{# of predictive biomarkers classified as predictive}} {\text{# of predictive biomarkers}}.\end{equation}\]
Function
## function (fit_results)
## {
## group_vars <- c(".dgp_name", ".method_name", "n")
## eval_out <- fit_results %>% mutate(tpp = purrr::map_dbl(tems,
## tpp_fun)) %>% dplyr::group_by(dplyr::across({
## {
## group_vars
## }
## })) %>% summarize(tpr = mean(tpp), .groups = "drop")
## return(eval_out)
## }Input Parameters
## list()Empirical TNR
The empirical true negative rate (TNR) is the mean of the true negativeproportions (TNP) stratified over each data-generating process and sample size.For any given simulated dataset, the TNP is computed as\[\begin{equation} \text{TNP} = \frac{\text{# of non-predictive biomarkers not classified as predictive}} {\text{# of non-predictive biomarkers}}.\end{equation}\]
Function
## function (fit_results)
## {
## group_vars <- c(".dgp_name", ".method_name", "n")
## eval_out <- fit_results %>% mutate(tnp = purrr::map_dbl(tems,
## tnp_fun)) %>% dplyr::group_by(dplyr::across({
## {
## group_vars
## }
## })) %>% summarize(tnr = mean(tnp), .groups = "drop")
## return(eval_out)
## }Input Parameters
## list()Visualizations
Summary Plot
The simulation study results are summarized in a multi-panel plot comparing theempirical FDR, TPR and TNR of all methods across a range of sample sizes.
Function
## function (eval_results)
## {
## emp_fdr_tbl <- eval_results$`Empirical FDR` %>% dplyr::mutate(metric = "FDR",
## value = fdr, yintercept = 0.05) %>% dplyr::select(-fdr)
## emp_tpr_tbl <- eval_results$`Empirical TPR` %>% dplyr::mutate(metric = "TPR",
## value = tpr, yintercept = NA) %>% dplyr::select(-tpr)
## emp_tnr_tbl <- eval_results$`Empirical TNR` %>% dplyr::mutate(metric = "TNR",
## value = tnr, yintercept = NA) %>% dplyr::select(-tnr)
## comb_tbl <- bind_rows(emp_fdr_tbl, emp_tpr_tbl, emp_tnr_tbl)
## plt <- ggplot2::ggplot(comb_tbl) + ggplot2::aes(x = n, y = value,
## colour = as.factor(.method_name)) + ggplot2::geom_point(alpha = 0.7) +
## ggplot2::geom_line(alpha = 0.7) + ggplot2::geom_hline(aes(yintercept = yintercept),
## colour = "red", linetype = 2, alpha = 0.5) + ggplot2::facet_grid(cols = ggplot2::vars(.dgp_name),
## rows = ggplot2::vars(metric)) + ggplot2::xlab("Sample Size") +
## ggplot2::ylab("Value") + ggplot2::scale_colour_viridis_d(name = "Method",
## option = "E", end = 0.8) + ggplot2::theme_bw()
## return(plt)
## }
## <bytecode: 0x7f940b7966c0>Input Parameters
## list()LM with TEM-Kinked with TEM-NLM with TEM
Varying n-n-n
Empirical FDR
Empirical TPR
Empirical TNR
Summary Plot
Parameter Values
## $dgp
## $dgp$`LM with TEM`
## $dgp$`LM with TEM`$n
## [1] 125 250 500
##
##
## $dgp$`Kinked with TEM`
## $dgp$`Kinked with TEM`$n
## [1] 125 250 500
##
##
## $dgp$`NLM with TEM`
## $dgp$`NLM with TEM`$n
## [1] 125 250 500
##
##
##
## $method
## list()