inst/doc/survival.R

In DynForest: Random Forest with Multivariate Longitudinal Predictors

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----message = FALSE----------------------------------------------------------
library("DynForest")
data(pbc2)
head(pbc2)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  set.seed(1234)
#  id <- unique(pbc2$id)
#  id_sample <- sample(id, length(id)*2/3)
#  id_row <- which(pbc2$id %in% id_sample)
#  pbc2_train <- pbc2[id_row,]
#  pbc2_pred <- pbc2[-id_row,]

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  timeData_train <- pbc2_train[,c("id","time",
#                                  "serBilir","SGOT",
#                                  "albumin","alkaline")]
#  fixedData_train <- unique(pbc2_train[,c("id","age","drug","sex")])

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  timeVarModel <- list(serBilir = list(fixed = serBilir ~ time,
#                                       random = ~ time),
#                       SGOT = list(fixed = SGOT ~ time + I(time^2),
#                                   random = ~ time + I(time^2)),
#                       albumin = list(fixed = albumin ~ time,
#                                      random = ~ time),
#                       alkaline = list(fixed = alkaline ~ time,
#                                       random = ~ time))

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  Y <- list(type = "surv",
#            Y = unique(pbc2_train[,c("id","years","event")]))

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  res_dyn <- dynforest(timeData = timeData_train,
#                       fixedData = fixedData_train,
#                       timeVar = "time", idVar = "id",
#                       timeVarModel = timeVarModel, Y = Y,
#                       ntree = 200, mtry = 3, nodesize = 2, minsplit = 3,
#                       cause = 2, ncores = 7, seed = 1234)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  summary(res_dyn)
#  
#  dynforest executed for survival (competing risk) outcome
#  	Splitting rule: Fine & Gray statistic test
#  	Out-of-bag error type: Integrated Brier Score
#  	Leaf statistic: Cumulative incidence function
#  ----------------
#  Input
#  	Number of subjects: 208
#  	Longitudinal: 4 predictor(s)
#  	Numeric: 1 predictor(s)
#  	Factor: 2 predictor(s)
#  ----------------
#  Tuning parameters
#  	mtry: 3
#  	nodesize: 2
#  	minsplit: 3
#  	ntree: 200
#  ----------------
#  ----------------
#  dynforest summary
#  	Average depth per tree: 6.62
#  	Average number of leaves per tree: 27.68
#  	Average number of subjects per leaf: 4.78
#  	Average number of events of interest per leaf: 1.95
#  ----------------
#  Computation time
#  	Number of cores used: 7
#  	Time difference of 3.18191 mins
#  ----------------

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  head(get_tree(dynforest_obj = res_dyn, tree = 1))
#  
#            type id_node var_split feature  threshold   N depth
#  1 Longitudinal       1         4       2 -0.6931377 123     1
#  2      Numeric       2         2      NA -1.4960589  25     2
#  3      Numeric       3         1      NA -1.0052000  98     2
#  4 Longitudinal       4         3       1  0.3332725   3     3
#  5       Factor       5         2      NA         NA  22     3
#  6 Longitudinal       6         5       2  0.7550109  10     3

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  tail(get_tree(dynforest_obj = res_dyn, tree = 1))
#  
#              type id_node var_split feature     threshold N depth
#  156         Leaf   14522        NA      NA            NA 1    14
#  157 Longitudinal   14523         1       1  8.894535e-03 3    14
#  158 Longitudinal   29046         6       1 -1.385950e-07 2    15
#  159         Leaf   29047        NA      NA            NA 1    15
#  160         Leaf   58092        NA      NA            NA 1    16
#  161         Leaf   58093        NA      NA            NA 1    16

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  plot(res_dyn, tree = 1, nodes = 251)

## ----fig.cap = "Figure 1: Estimated cumulative incidence functions in tree 1 and node 251.", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_CIF.png")

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  plot(res_dyn, id = 104, max_tree = 9)

## ----DynForestRCIF, fig.cap = "Figure 2: Estimated cumulative incidence functions for subject 104 over 9 trees.", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_CIF9.png")

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  res_dyn_OOB <- compute_ooberror(dynforest_obj = res_dyn)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  res_dyn_OOB
#  
#  [1] 0.1265053

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  id_pred <- unique(pbc2_pred$id[which(pbc2_pred$years>4)])
#  pbc2_pred_tLM <- pbc2_pred[which(pbc2_pred$id %in% id_pred),]
#  timeData_pred <- pbc2_pred_tLM[,c("id","time",
#                                    "serBilir","SGOT",
#                                    "albumin","alkaline")]
#  fixedData_pred <- unique(pbc2_pred_tLM[,c("id","age","drug","sex")])
#  pred_dyn <- predict(object = res_dyn,
#                      timeData = timeData_pred,
#                      fixedData = fixedData_pred,
#                      idVar = "id", timeVar = "time",
#                      t0 = 4)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  plot(pred_dyn, id = c(102, 260))

## ----DynForestRpredCIF, fig.cap = "Figure 3: Predicted cumulative incidence function for subjects 102 and 260 from landmark time of 4 years (represented by the dashed vertical line)", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_CIF9.png")

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  res_dyn_VIMP <- compute_vimp(dynforest_obj = res_dyn, seed = 123)

## ----eval = FALSE, echo = TRUE, fig.show='hide'-------------------------------
#  p1 <- plot(res_dyn_VIMP, PCT = TRUE)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  group <- list(group1 = c("serBilir","SGOT"),
#                group2 = c("albumin","alkaline"))
#  res_dyn_gVIMP <- compute_gvimp(dynforest_obj = res_dyn,
#                                 group = group, seed = 123)

## ----eval = FALSE, echo = TRUE, fig.show='hide'-------------------------------
#  p2 <- plot(res_dyn_gVIMP, PCT = TRUE)

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  plot_grid(p1, p2, labels = c("A", "B"))

## ----DynForestRVIMPgVIMP, fig.cap = "Figure 4: (A) VIMP statistic and (B) grouped-VIMP statistic displayed as a percentage of loss in OOB error of prediction. group1 includes serBilir and SGOT; group2 includes albumin and alkaline.", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_VIMP_gVIMP.png")

## ----eval = FALSE, echo = TRUE------------------------------------------------
#  res_dyn_max <- dynforest(timeData = timeData_train,
#                           fixedData = fixedData_train,
#                           timeVar = "time", idVar = "id",
#                           timeVarModel = timeVarModel, Y = Y,
#                           ntree = 200, mtry = 7, nodesize = 2, minsplit = 3,
#                           cause = 2, ncores = 7, seed = 1234)

## ----echo = TRUE, eval = FALSE, fig.show='hide'-------------------------------
#  depth_dyn <- compute_vardepth(dynforest_obj = res_dyn_max)
#  p1 <- plot(depth_dyn, plot_level = "predictor")
#  p2 <- plot(depth_dyn, plot_level = "feature")

## ----echo = TRUE, eval = FALSE------------------------------------------------
#  plot_grid(p1, p2, labels = c("A", "B"))

## ----DynForestRmindepth, fig.cap = "Figure 5: Average minimal depth level by predictor (A) and feature (B).", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_mindepth.png")

## ----DynForestRmtrytuned, fig.cap = "Figure 6: OOB error according to mtry hyperparameter", eval = TRUE, echo = FALSE, out.width="70%"----
knitr::include_graphics("Figures/DynForestR_surv_mtrytuned.png")