#' Moderated Statistical Tests for Influence Functions

#'

#' Performs variance shrinkage via application of an empirical Bayes procedure

#' (of LIMMA) on the observed data after a transformation moving the data to

#' influence function space, based on the average treatment effect parameter.

#'

#' @param biotmle \code{biotmle} object as generated by \code{biomarkertmle}

#' @param adjust the multiple testing correction to be applied to p-values that

#'  are generated from the moderated tests. The recommended (default) method

#'  is that of Benjamini and Hochberg. See \link[limma]{topTable} for a list of

#'  appropriate methods.

#' @param pval_type The reference distribution to be used for computing the

#'  p-value. Those based on the normal approximation tend to provide misleading

#'  inference when working with moderately sized (finite) samples. Use of the

#'  logistic distribution has been found to empirically improve performance in

#'  settings where multiple hypothesis testing is a concern.

#' @param ... Other arguments passed to \code{\link[limma]{topTable}}.

#'

#' @importFrom methods is

#' @importFrom dplyr "%>%" mutate select

#' @importFrom stats plogis p.adjust

#' @importFrom limma lmFit eBayes topTable

#' @importFrom tibble as_tibble rownames_to_column

#' @importFrom assertthat assert_that

#'

#' @return \code{biotmle} object containing the results of applying both

#'  \code{\link[limma]{lmFit}} and \code{\link[limma]{topTable}}.

#'

#' @export modtest_ic

#'

#' @examples

#' library(dplyr)

#' library(biotmleData)

#' library(SuperLearner)

#' library(SummarizedExperiment)

#' data(illuminaData)

#'

#' colData(illuminaData) <- colData(illuminaData) %>%

#'   data.frame() %>%

#'   dplyr::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")

#' )

#'

#' limmaTMLEout <- modtest_ic(biotmle = biomarkerTMLEout)

modtest_ic <- function(biotmle,

                       adjust = "BH",

                       pval_type = c("normal", "logistic"),

                       ...) {

  # check for input type and set argument defaults

  assertthat::assert_that(is(biotmle, "bioTMLE"))

  pval_type <- match.arg(pval_type)

  # extract influence function estimates and number of observations

  biomarkertmle_out <- as.matrix(biotmle@tmleOut)

  n_obs <- nrow(colData(biotmle))

  # build design matrix for forcing limma to only perform shrinkage

  fit <- limma::lmFit(

    object = biomarkertmle_out,

    design = rep(1, n_obs)

  )

  fit <- limma::eBayes(fit = fit)

  # compute inference

  tt <- limma::topTable(

    fit = fit,

    coef = 1,

    number = Inf,

    adjust.method = adjust,

    sort.by = "none",

    ...

  ) %>%

    tibble::as_tibble(.name_repair = "minimal") %>%

    dplyr::mutate(

      # remove log-fold change and overwrite average expression with ATE

      ID = biotmle@NAMES,

      logFC = NULL,

      AveExpr = biotmle@ateOut,

      var_bayes = fit$s2.post / n_obs,

    ) %>%

    dplyr::select(ID, AveExpr, t, P.Value, adj.P.Val, B, var_bayes)

  # use logistic distribution as reference for p-values by default

  if (pval_type == "logistic") {

    p_val <- 2 * stats::plogis(-abs(tt$t), location = 0, scale = sqrt(3) / pi)

    p_val_adj <- stats::p.adjust(p_val, method = adjust)

    tt$P.Value <- p_val

    tt$adj.P.Val <- p_val_adj

  }

  # store in slot of output object

  biotmle@topTable <- tt

  return(biotmle)

}