SM_PRISM.R
In StratifiedMedicine: Stratified Medicine

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=10, fig.height=8.5
)

## ----sim_ctns, warning=FALSE, message=FALSE-----------------------------------
library(ggplot2)
library(dplyr)
library(partykit)
library(StratifiedMedicine)
library(survival)
dat_ctns = generate_subgrp_data(family="gaussian")
Y = dat_ctns$Y
X = dat_ctns$X # 50 covariates, 46 are noise variables, X1 and X2 are truly predictive
A = dat_ctns$A # binary treatment, 1:1 randomized 
length(Y)
table(A)
dim(X)

## ----filter_glmnet, warning=FALSE, message=FALSE------------------------------
res_f <- filter_train(Y, A, X, filter="glmnet")
res_f$filter.vars
plot_importance(res_f)

## ----filter_glmnet2, warning=FALSE, message=FALSE, include=FALSE--------------
res_f2 <- filter_train(Y, A, X, filter="glmnet", hyper=list(interaction=T))
res_f2$filter.vars
plot_importance(res_f2)

## ----ple_train, warning=FALSE, message=FALSE----------------------------------
res_p1 <- ple_train(Y, A, X, ple="ranger", meta="T-learner")
summary(res_p1$mu_train)
plot_ple(res_p1)
plot_dependence(res_p1, X=X, vars="X1")

## ----ple_train2, warning=FALSE, message=FALSE---------------------------------
res_p2 <- ple_train(Y, A, X, ple="ranger", meta="T-learner", hyper=list(mtry=5))
plot_dependence(res_p2, X=X, vars=c("X1", "X2"))

## ----submod_train1, warning=FALSE, message=FALSE------------------------------
res_s1 <- submod_train(Y, A, X, submod="lmtree")
summary(res_s1)
plot_tree(res_s1)

## ----submod_train2, warning=FALSE, message=FALSE------------------------------
res_s2 <- submod_train(Y, A, X,  mu_train=res_p2$mu_train, submod="rpart_cate")
summary(res_s2)

## ----param1, warning=FALSE, message=FALSE-------------------------------------
param.dat1 <- param_est(Y, A, X, Subgrps = res_s1$Subgrps.train, param="lm")
param.dat1

## ----table_steps, echo=FALSE--------------------------------------------------
library(knitr)
summ.table = data.frame( `Step` = c("estimand(s)", "filter", "ple", "submod", "param"),
                        `gaussian` = c("E(Y|A=0)<br>E(Y|A=1)<br>E(Y|A=1)-E(Y|A=0)",
                                       "Elastic Net<br>(glmnet)", 
                                       "X-learner: Random Forest<br>(ranger)",
                                       "MOB(OLS)<br>(lmtree)", 
                                       "Double Robust<br>(dr)"),
                        `binomial` = c("E(Y|A=0)<br>E(Y|A=1)<br>E(Y|A=1)-E(Y|A=0)",
                                       "Elastic Net<br>(glmnet)", 
                                       "X-learner: Random Forest<br>(ranger)",
                                       "MOB(GLM)<br>(glmtree)", 
                                       "Doubly Robust<br>(dr)"),    
                        `survival` = c("HR(A=1 vs A=0)",
                                       "Elastic Net<br>(glmnet)", 
                                       "T-learner: Random Forest<br>(ranger)",
                                       "MOB(OLS)<br>(lmtree)", 
                                       "Hazard Ratios<br>(cox)") )                        
                      
kable( summ.table, caption = "Default PRISM Configurations (With Treatment)", full_width=T)

summ.table = data.frame(`Step` = c("estimand(s)", "filter", "ple", "submod", "param"),
                        `gaussian` = c("E(Y)",
                                       "Elastic Net<br>(glmnet)", 
                                       "Random Forest<br>(ranger)",
                                       "Conditional Inference Trees<br>(ctree)",
                                       "OLS<br>(lm)"),
                        `binomial` = c("Prob(Y)",
                                       "Elastic Net<br>(glmnet)", 
                                       "Random Forest<br>(ranger)",
                                       "Conditional Inference Trees<br>(ctree)", 
                                       "OLS<br>(lm)"),    
                        `survival` = c("RMST", "Elastic Net<br>(glmnet)", 
                                       "Random Forest<br>(ranger)",
                                       "Conditional Inference Trees<br>(ctree)",
                                       "RMST<br>(rmst)"))                        
                      
kable( summ.table, caption = "Default PRISM Configurations (Without Treatment, A=NULL)", full_width=T)

## ----default_ctns, warning=FALSE----------------------------------------------
# PRISM Default: filter_glmnet, ranger, lmtree, dr #
res0 = PRISM(Y=Y, A=A, X=X)
summary(res0)
plot(res0) # same as plot(res0, type="tree")

## ----default_ctns_prog, warning=FALSE-----------------------------------------
# PRISM Default: filter_glmnet, ranger, ctree, param_lm #
res_prog = PRISM(Y=Y, X=X)
# res_prog = PRISM(Y=Y, A=NULL, X=X) #also works
summary(res_prog)

## ----default_ctns_filter, include=FALSE---------------------------------------
# Elastic net model: loss by lambda #
plot(res0$filter.mod)
# Variables that remain after filtering #
res0$filter.vars
# All predictive variables (X1,X2) and prognostic variables (X3,X5, X7) remains.
plot_importance(res0)

## ----default_ctns_ple---------------------------------------------------------
summary(res0$mu_train)
plot_ple(res0)
plot_dependence(res0, vars=c("X2"))

## ----default_ctns_submod------------------------------------------------------
plot(res0$submod.fit$mod, terminal_panel = NULL)
table(res0$out.train$Subgrps)
table(res0$out.test$Subgrps)

## ----default_ctns2------------------------------------------------------------
## Overall/subgroup specific parameter estimates/inference
res0$param.dat

## ----default_hyper, eval=FALSE------------------------------------------------
#  res_new_hyper = PRISM(Y=Y, A=A, X=X, filter.hyper = list(lambda="lambda.1se"),
#                        ple.hyper = list(min.node.pct=0.05),
#                        submod.hyper = list(minsize=200, maxdepth=3), verbose=FALSE)
#  summary(res_new_hyper)

## ----default_binary-----------------------------------------------------------
dat_bin = generate_subgrp_data(family="binomial", seed = 5558)
Y = dat_bin$Y
X = dat_bin$X # 50 covariates, 46 are noise variables, X1 and X2 are truly predictive
A = dat_bin$A # binary treatment, 1:1 randomized 

res0 = PRISM(Y=Y, A=A, X=X)
summary(res0)

## ----default_surv-------------------------------------------------------------
# Load TH.data (no treatment; generate treatment randomly to simulate null effect) ##
data("GBSG2", package = "TH.data")
surv.dat = GBSG2
# Design Matrices ###
Y = with(surv.dat, Surv(time, cens))
X = surv.dat[,!(colnames(surv.dat) %in% c("time", "cens")) ]
set.seed(6345)
A = rbinom(n = dim(X)[1], size=1, prob=0.5)

# Default: glmnet ==> ranger (estimates patient-level RMST(1 vs 0) ==> mob_weib (MOB with Weibull) ==> cox (Cox regression)
res_weib = PRISM(Y=Y, A=A, X=X)
summary(res_weib)
plot(res_weib)

## ----default_boot, warning=FALSE, message=FALSE, eval=FALSE-------------------
#  res_boot = PRISM(Y=Y, A=A, X=X, resample = "Bootstrap", R=50, ple="None")
#  summary(res_boot)
#  # Plot of distributions #
#  plot(res_boot, type="resample", estimand = "HR(A=1 vs A=0)")+geom_vline(xintercept = 1)