## ----setup_data----------------------------------------------------------
library(dplyr)
library(biotmle)
library(biotmleData)
data(illuminaData)
library(SummarizedExperiment)
"%ni%" = Negate("%in%")

## ----clean_data----------------------------------------------------------
# discretize "age" in the phenotype-level data
colData(illuminaData) <- colData(illuminaData) %>%
  data.frame %>%
  dplyr::mutate(age = as.numeric(age > median(age))) %>%
  DataFrame

# specify column index of treatment/exposure variable of interest
varInt_index <- which(names(colData(illuminaData)) %in% "benzene")

## ----biomarkerTMLE_eval, eval=FALSE--------------------------------------
#  biomarkerTMLEout <- biomarkertmle(se = illuminaData,
#                                    varInt = varInt_index,
#                                    type = "exposure",
#                                    family = "gaussian",
#                                    g_lib = c("SL.glmnet", "SL.randomForest",
#                                              "SL.polymars", "SL.mean"),
#                                    Q_lib = c("SL.glmnet", "SL.randomForest",
#                                              "SL.nnet", "SL.mean")
#                                   )

## ----load_biomarkerTMLE_result, echo=FALSE-------------------------------
data(biomarkertmleOut)

## ----limmaTMLE_eval------------------------------------------------------
varInt_index <- which(names(colData(illuminaData)) %in% "benzene")
designVar <- as.data.frame(colData(illuminaData))[, varInt_index]
design <- as.numeric(designVar == max(designVar))

limmaTMLEout <- modtest_ic(biotmle = biomarkerTMLEout, design = design)

## ----pval_hist_limma_adjp------------------------------------------------
plot(x = limmaTMLEout, type = "pvals_adj")

## ----pval_hist_limma_rawp------------------------------------------------
plot(x = limmaTMLEout, type = "pvals_raw")

## ----heatmap_limma_results-----------------------------------------------
heatmap_ic(x = limmaTMLEout, design = design, FDRcutoff = 0.05, top = 25)

## ----volcano_plot_limma_results------------------------------------------
volcano_ic(biotmle = limmaTMLEout)

## ----sessionInfo, echo=FALSE---------------------------------------------
sessionInfo()