R/biotmle.R

In nhejazi/biotmle: Targeted Learning with Moderated Statistics for Biomarker Discovery

utils::globalVariables(c("assay<-"))

#' Biomarker Evaluation with Targeted Minimum Loss Estimation of the ATE
#'
#' Computes the causal target parameter defined as the difference between the
#' biomarker expression values under treatment and those same values under no
#' treatment, using Targeted Minimum Loss Estimation.
#'
#' @param se A \code{SummarizedExperiment} containing microarray expression
#'  or next-generation sequencing data in the \code{assays} slot and a matrix
#'  of phenotype-level data in the \code{colData} slot.
#' @param varInt A \code{numeric} indicating the column of the design matrix
#'  corresponding to the treatment or outcome of interest (in the
#'  \code{colData} slot of the \code{SummarizedExperiment} argument "se").
#' @param normalized A \code{logical} indicating whether the data included in
#'  the \code{assay} slot of the input \code{SummarizedExperiment} object has
#'  been normalized externally. The default is set to \code{TRUE} with the
#'  expectation that an appropriate normalization method has been applied. If
#'  set to \code{FALSE}, median normalization is performed for microarray data.
#' @param ngscounts A \code{logical} indicating whether the data are counts
#'  generated from a next-generation sequencing experiment (e.g., RNA-seq). The
#'  default setting assumes continuous expression measures as generated by
#'  microarray platforms.
#' @param bppar_type A parallelization option specified by \code{BiocParallel}.
#'  Consult the manual page for \code{\link[BiocParallel]{BiocParallelParam}}
#'  for possible types and their descriptions. The default for this argument is
#'  \code{\link[BiocParallel]{MulticoreParam}}, for multicore evaluation.
#' @param bppar_debug A \code{logical} indicating whether or not to rely upon
#'  pkg{BiocParallel}. Setting this argument to \code{TRUE}, replaces the call
#'  to \code{\link[BiocParallel]{bplapply}} by a call to \code{lapply}, which
#'  significantly reduces the overhead of debugging. Note that invoking this
#'  option overrides all other parallelization arguments.
#' @param cv_folds A \code{numeric} scalar indicating how many folds to use in
#'  performing targeted minimum loss estimation. Cross-validated estimates have
#'  been demonstrated to allow relaxation of certain theoretical conditions and
#'  and accommodate the construction of more conservative variance estimates.
#' @param g_lib A \code{character} vector specifying the library of machine
#'  learning algorithms for use in fitting the propensity score P(A = a | W).
#' @param Q_lib A \code{character} vector specifying the library of machine
#'  learning algorithms for use in fitting the outcome regression E[Y | A,W].
#' @param ... Additional arguments to be passed to \code{\link[drtmle]{drtmle}}
#'  in computing the targeted minimum loss estimator of the average treatment
#'  effect.
#'
#' @importFrom SummarizedExperiment assay colData rowData SummarizedExperiment
#' @importFrom BiocParallel register bplapply bpprogressbar MulticoreParam
#' @importFrom tibble as_tibble
#'
#' @return S4 object of class \code{biotmle}, inheriting from
#'  \code{SummarizedExperiment}, with additional slots \code{tmleOut} and
#'  \code{call}, among others, containing TML estimates of the ATE of exposure
#'  on biomarker expression.
#'
#' @export biomarkertmle
#'
#' @examples
#' library(dplyr)
#' library(biotmleData)
#' library(SuperLearner)
#' library(SummarizedExperiment)
#' data(illuminaData)
#'
#' colData(illuminaData) <- colData(illuminaData) %>%
#'   data.frame() %>%
#'   mutate(age = as.numeric(age > median(age))) %>%
#'   DataFrame()
#' benz_idx <- which(names(colData(illuminaData)) %in% "benzene")
#'
#' biomarkerTMLEout <- biomarkertmle(
#'   se = illuminaData[1:2, ],
#'   varInt = benz_idx,
#'   bppar_type = BiocParallel::SerialParam(),
#'   g_lib = c("SL.mean", "SL.glm"),
#'   Q_lib = c("SL.mean", "SL.glm")
#' )
biomarkertmle <- function(se,
                          varInt,
                          normalized = TRUE,
                          ngscounts = FALSE,
                          bppar_type = BiocParallel::MulticoreParam(),
                          bppar_debug = FALSE,
                          cv_folds = 1,
                          g_lib = c(
                            "SL.mean", "SL.glm", "SL.bayesglm"
                          ),
                          Q_lib = c(
                            "SL.mean", "SL.bayesglm", "SL.earth", "SL.ranger"
                          ),
                          ...) {

  # catch input and invoke S4 class constructor for "bioTMLE" object
  call <- match.call(expand.dots = TRUE)
  biotmle <- .biotmle(
    SummarizedExperiment(
      assays = list(expMeasures = assay(se)),
      rowData = rowData(se),
      colData = colData(se)
    ),
    call = call,
    tmleOut = tibble::as_tibble(matrix(NA, 10, 10), .name_repair = "minimal"),
    topTable = tibble::as_tibble(matrix(NA, 10, 10), .name_repair = "minimal")
  )

  # invoke the voom transform from LIMMA if next-generation sequencing data)
  if (ngscounts) {
    voom_out <- rnaseq_ic(biotmle)
    voom_exp <- 2^(voom_out$E)
    assay(se) <- voom_exp
  }

  # set up parallelization based on input
  BiocParallel::bpprogressbar(bppar_type) <- TRUE
  BiocParallel::register(bppar_type, default = TRUE)

  # TMLE procedure to identify biomarkers based on an EXPOSURE
  if (!ngscounts && !normalized) {
    # median normalization
    exp_normed <- limma::normalizeBetweenArrays(as.matrix(assay(se)),
      method = "scale"
    )
    Y <- tibble::as_tibble(t(exp_normed), .name_repair = "minimal")
  } else {
    Y <- tibble::as_tibble(t(as.matrix(assay(se))), .name_repair = "minimal")
  }
  # simple sanity check of whether Y includes array values
  if (!all(apply(Y, 2, class) == "numeric")) {
    stop("Warning - values in Y do not appear to be numeric.")
  }

  # exposure / treatment
  A <- as.numeric(SummarizedExperiment::colData(se)[, varInt])

  # baseline covariates
  W <- tibble::as_tibble(SummarizedExperiment::colData(se)[, -varInt],
                         .name_repair = "minimal")
  if (is.null(dim(W)[2])) {
    W <- as.numeric(rep(1, length(A)))
  }

  # coerce matrix of baseline covariates to numeric
  if (!all(is.numeric(apply(W, 2, class)))) {
    W <- tibble::as_tibble(apply(W, 2, as.numeric), .name_repair = "minimal")
  }

  # perform multi-level TMLE (of the ATE) for genes as Y
  if (!bppar_debug) {
    biomarkertmle_out <- BiocParallel::bplapply(Y[, seq_along(Y)],
      exp_biomarkertmle,
      W = W,
      A = A,
      g_lib = g_lib,
      Q_lib = Q_lib,
      cv_folds = cv_folds,
      ...
    )
  } else {
    biomarkertmle_out <- lapply(Y[, seq_along(Y)],
      exp_biomarkertmle,
      W = W,
      A = A,
      g_lib = g_lib,
      Q_lib = Q_lib,
      cv_folds = cv_folds,
      ...
    )
  }
  biomarkertmle_params <- do.call(c, lapply(biomarkertmle_out, `[[`, "param"))
  biomarkertmle_eifs <- do.call(
    cbind.data.frame,
    lapply(biomarkertmle_out, `[[`, "eif")
  )

  biotmle@ateOut <- as.numeric(biomarkertmle_params)
  if (!ngscounts) {
    biomarker_eifs <- t(as.matrix(biomarkertmle_eifs))
    colnames(biomarker_eifs) <- colnames(se)
    biotmle@tmleOut <- tibble::as_tibble(
      biomarker_eifs,
      .name_repair = "minimal"
    )
  } else {
    voom_out$E <- t(as.matrix(biomarkertmle_eifs))
    biotmle@tmleOut <- voom_out
  }
  return(biotmle)
}

###############################################################################

#' TMLE procedure using ATE for Biomarker Identication from Exposure
#'
#' This function performs influence curve-based estimation of the effect of an
#' exposure on biological expression values associated with a given biomarker,
#' controlling for a user-specified set of baseline covariates.
#'
#' @param Y A \code{numeric} vector of expression values for a given biomarker.
#' @param A A \code{numeric} vector of discretized exposure vector (e.g., from
#'  a design matrix whose effect on expression values is of interest.
#' @param W A \code{Matrix} of \code{numeric} values corresponding to baseline
#'  covariates to be marginalized over in the estimation process.
#' @param g_lib A \code{character} vector identifying the library of learning
#'  algorithms to be used in fitting the propensity score P[A = a | W].
#' @param Q_lib A \code{character} vector identifying the library of learning
#'  algorithms to be used in fitting the outcome regression E[Y | A, W].
#' @param cv_folds A \code{numeric} scalar indicating how many folds to use in
#'  performing targeted minimum loss estimation. Cross-validated estimates are
#'  more robust, allowing relaxing of theoretical conditions and construction
#'  of conservative variance estimates.
#' @param ... Additional arguments passed to \code{\link[drtmle]{drtmle}} in
#'  computing the targeted minimum loss estimator of the average treatment
#'  effect.
#'
#' @importFrom assertthat assert_that
#' @importFrom drtmle drtmle
#'
#' @return TMLE-based estimate of the relationship between biomarker expression
#'  and changes in an exposure variable, computed iteratively and saved in the
#'  \code{tmleOut} slot in a \code{biotmle} object.
exp_biomarkertmle <- function(Y,
                              A,
                              W,
                              g_lib,
                              Q_lib,
                              cv_folds,
                              ...) {
  # check the case that Y is passed in as a column of a data.frame
  if (any(class(Y) == "data.frame")) Y <- as.numeric(unlist(Y[, 1]))
  if (any(class(A) == "data.frame")) A <- as.numeric(unlist(A[, 1]))
  assertthat::assert_that(length(unique(A)) > 1)

  # fit standard (possibly CV) TML estimator (n.b., guard = NULL)
  a_0 <- sort(unique(A[!is.na(A)]))
  suppressWarnings(
    tmle_fit <- drtmle::drtmle(
      Y = Y,
      A = A,
      W = W,
      a_0 = a_0,
      SL_g = g_lib,
      SL_Q = Q_lib,
      cvFolds = cv_folds,
      stratify = TRUE,
      guard = NULL,
      parallel = FALSE,
      use_future = FALSE,
      ...
    )
  )

  # compute ATE and estimated EIF by delta method
  ate_tmle <- tmle_fit$tmle$est[seq_along(a_0)[-1]] - tmle_fit$tmle$est[1]
  eif_tmle_delta <- tmle_fit$ic$ic[, seq_along(a_0)[-1]] - tmle_fit$ic$ic[, 1]

  # return only highest contrast (e.g., a[1] v a[5]) if many contrasts
  if (!is.vector(eif_tmle_delta)) {
    param_out <- ate_tmle[length(ate_tmle)]
    eif_out <- eif_tmle_delta[, ncol(eif_tmle_delta)] +
      ate_tmle[length(ate_tmle)]
  } else {
    param_out <- ate_tmle
    eif_out <- eif_tmle_delta + ate_tmle
  }
  assertthat::assert_that(is.vector(eif_out))

  # output
  out <- list(param = param_out, eif = eif_out)
  return(out)
}