library("subtee")

################################################################################
## Load the needed packages for saving the figures
library("ggpubr")
library("ggplot2")

################################################################################
## Load the needed dataset

## We use the dataset from Rosenkranz (2016)
## https://onlinelibrary.wiley.com/doi/abs/10.1002/bimj.201500147 to
## illustrate the methods proposed in this work.

## The data comes from a clinical trial of an prostate cancer
## treatment.

## Data is loaded from
## Royston, Patrick, and Willi Sauerbrei. Multivariable
## model-building: a pragmatic approach to regression anaylsis based
## on fractional polynomials for modelling continuous
## variables. Vol. 777. John Wiley & Sons,
## 2008. https://www.imbi.uni-freiburg.de/Royston-Sauerbrei-book

prca <- get_prca_data()

################################################################################
## Section 4.1 ------------------------------------------------------------------
## first create candidate subgroups
cand.groups <- subbuild(prca, BM == 1, PF == 1, HX == 1, STAGE == 4, 
  AGE > 65, WT > 100)
fitdat <- cbind(prca, cand.groups)

################################################################################
## Section 4.2 ------------------------------------------------------------------
subgr.names <- names(cand.groups)
prog <- paste0("`", subgr.names,"`", collapse = " + ")
prog <- as.formula(paste(" ~ ", prog))
res_unadj <- unadj(resp = "SURVTIME", trt = "RX", subgr = subgr.names,
  data = fitdat, covars = prog, event = "CENS", fitfunc = "coxph", 
  ties = "breslow")
res_unadj

confint(res_unadj, level = 0.80) # Results are not shown in the mansucript

summary(res_unadj)               # Results are not shown in the mansucript

plot(res_unadj, show.compl = TRUE)

plot(res_unadj, type = "trtEffDiff")


## Code to save figure side to side
## Caption: Figure 1. Forest plot of unadjusted treatment effects for subgroups and complements 
## (left) and treatment-subgroup interactions (right). 
## These graphics are the output of the \\code{plot}
## function applied to a \\code{subtee} object
p1 <- plot(res_unadj, show.compl = TRUE)
p2 <- plot(res_unadj, type = "trtEffDiff")
ggpubr::ggarrange(p1, p2, widths = c(.20,.15))
ggplot2::ggsave(filename = "01_unadj.pdf", width = 12, height = 6)

###############################################################################
## Section 4.3 ------------------------------------------------------------------
res_modav = modav(resp = "SURVTIME", trt = "RX", subgr = subgr.names,
  data = fitdat,  covars = prog,  event = "CENS",  fitfunc = "coxph", 
  ties = "breslow")
res_modav

plot(res_unadj, res_modav, show.compl = TRUE)

plot(res_unadj, res_modav, type = "trtEffDiff")

## Code to save figure side to side
## Caption: Figure 2. Forest plot of treatment effects for subgroups and complements (left) 
## and treatment-subgroup interactions (right). The plots compare the 
## results using unadjusted estimates and those obtained with model averaging.
## These graphics are the output of the \\code{plot} function using
## two \\code{subtee} objects resulting from the \\code{unadj} and 
## the \\code{modav} functions
p1 <- plot(res_unadj, res_modav, show.compl = TRUE)
p2 <- plot(res_unadj, res_modav, type = "trtEffDiff")
p3 <- ggpubr::ggarrange(NULL, p2, NULL, ncol = 1, nrow = 3, heights = c(1, 3, 1))
ggpubr::ggarrange(p1, p3, widths = c(1.3,1))
ggplot2::ggsave(filename = "02_modav.pdf", width = 12, height = 6)

################################################################################
## Section 4.4 ------------------------------------------------------------------
set.seed(46312)
res_bagged <- bagged(resp = "SURVTIME", trt = "RX", subgr = subgr.names,
  data = fitdat, covars = prog, event = "CENS", fitfunc = "coxph", 
  ties = "breslow", select.by = "BIC", B = 2000)
res_bagged
res_bagged$boot_results
sum(res_bagged$boot_results$percSel)
summary(res_bagged)  # Results are not shown in the mansucript

plot(res_unadj, res_modav, res_bagged, show.compl = T)

plot(res_unadj, res_modav, res_bagged, type = "trtEffDiff")

## Code to save figure side to side
## Forest plot of treatment effects for subgroups and complements (left) and 
## treatment-subgroup interactions (right) using unadjusted, model averaging 
## and bagged estimates. These graphics are the output of the \\code{plot} 
## function applied to three \\code{subtee} objects.
p1 <- plot(res_unadj, res_modav, res_bagged, show.compl = TRUE)
p2 <- plot(res_unadj, res_modav, res_bagged, type = "trtEffDiff")
p3 <- ggpubr::ggarrange(NULL, p2, NULL, ncol = 1, nrow = 3, heights = c(1, 2, 1))
ggpubr::ggarrange(p1, p3, widths = c(1.1, 1))
ggsave(filename = "03_bagged.pdf", width = 12, height = 4)