---
title: "Guide to the Cox Proportional Hazards model"
author: "James Hickey"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{CoxPH}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
```{r, message=FALSE, echo=FALSE}
library("gbm3")
```

`gbm` facilitates the creation of boosted Cox proportional hazards models, a particularly useful feature when dealing with survival data.  The package can handle two types of survival object as the response, namely a right censored or counting survival object.  Both of these objects can be created using the `Surv` command found in the `survival` package. 

# Set-up of data and distribution object 
Right censored survival data consists of a time to event number and the event indicator - 0 if no event has taken place and 1 if the event has happened. On the other hand, counting survival data contains start and stop times along with an event indicator again indicating if an event has taken place in that period.  The data may be organised into strata and this should be passed to `gbm_dist` on creation of the `CoxPHGBMDist` object - see  the "Model Specific Parameters" vignette for more details.  The dataset considered here is provided by the `survival` package.

```{r, message=FALSE}
## Install package
require(survival)

# get datasets
right_cens <- cgd[cgd$enum==1, ]
start_stop <- cgd

# Set up GBMDist objects
right_cens_dist <- gbm_dist("CoxPH", strata=right_cens$hos.cat)
start_stop_dist <- gbm_dist("CoxPH", strata=start_stop$hos.cat)
```

# Creating a boosted model
Now to create the underlying boosted model the training parameters need to be defined and `gbmt` called. In this instance the data has observation ids associated with it and so it is necessary to create specific `GBMTrainParams` objects rather than relying on the defaults.

```{r, messagee=FALSE}
# Set-up training parameters
params_right_cens <- training_params(num_trees = 2000, interaction_depth = 3, 
                                     id=right_cens$id,
                                     num_train=round(0.5 * length(unique(right_cens$id))) )
params_start_stop <- training_params(num_trees = 2000, interaction_depth = 3, 
                                     id=start_stop$id,
                                     num_train=round(0.5 * length(unique(start_stop$id))) )

# Call to gbmt
fit_right_cens <- gbmt(Surv(tstop, status)~ age + sex + inherit +
                     steroids + propylac, data=right_cens, 
                     distribution = right_cens_dist,
                     train_params = params_right_cens, cv_folds=10,
                     keep_gbm_data = TRUE)
fit_start_stop <- gbmt(Surv(tstart, tstop, status)~ age + sex + inherit +
                     steroids + propylac, data=start_stop, 
                     distribution = start_stop_dist,
                     train_params = params_start_stop, cv_folds=10, 
                     keep_gbm_data = TRUE)

# Plot performance
best_iter_right <- gbmt_performance(fit_right_cens, method='test')
best_iter_stop_start <- gbmt_performance(fit_start_stop, method='test')
```
# Strata Updates
During the fitting process the original strata vector is updated in the following way.
When the data is split into a training and validation set the strata vector is also split.  The strata vector is then updated so as represent the cumulative count of the number of observations in each stratum in the training and validation sets.  The vector is padded with `NAs` so it is of the same length as the original strata vector provided and such that the validation set cumulative strata sums are separated from the training set strata counts by the appropriate amount.  

The original strata vector is stored within the `GBMFit` object and can be accessed as follows: `fit$distribution$original_strata_id`. The data in the `original_strata_id` field 
is used to recreate the correct strata when performing additional iterations using `gbm_more`.


# Role of additional parameters in GBMDist
## ties and prior_node_coeff_var
The `ties` and `prior_node_coeff_var` parameters may also be specified on construction of the `CoxPHGBMDist` object.  The former is a string specifying the method by which the algorithm deals with tied event times.  This may be set to either "breslow" or "efron" depending on your preference, with the latter being the default.  The role of the `prior_node_coeff_var` parameter is slightly more subtle and complex.  When fitting a boosted tree, the optimal predictions of the terminal nodes must be set.  These predictions determine the predictions made by the `GBMFit` object.  The role of `prior_node_coeff_var` is to ensure that the predictions are finite and it does this by acting as a regularization for the terminal node predictions.  It should be a finite positive double and is by default set to a 1000.  An exact description of its role in the underlying algorithm is described in the next section.

```{r, message=FALSE, fig.height=6, fig.width=6}
# Example using Breslow and Efron tie-breaking

# Create data
require(survival)

set.seed(1)
N <- 3000
X1 <- runif(N)
X2 <- runif(N)
X3 <- factor(sample(letters[1:4],N,replace=T))
mu <- c(-1,0,1,2)[as.numeric(X3)]

f <- 0.5*sin(3*X1 + 5*X2^2 + mu/10)
tt.surv <- rexp(N,exp(f))
tt.cens <- rexp(N,0.5)
delta <- as.numeric(tt.surv <= tt.cens)
tt <- apply(cbind(tt.surv,tt.cens),1,min)

# throw in some missing values
X1[sample(1:N,size=100)] <- NA
X3[sample(1:N,size=300)] <- NA

# random weights if you want to experiment with them
w <- rep(1,N)
data <- data.frame(tt=tt,delta=delta,X1=X1,X2=X2,X3=X3)

# Set up distribution objects
cox_breslow <- gbm_dist("CoxPH", ties="breslow", prior_node_coeff_var=100)
cox_efron <- gbm_dist("CoxPH", ties="efron", prior_node_coeff_var=100)

# Define training parameters
params <- training_params(num_trees=3000, interaction_depth=3, min_num_obs_in_node=10, 
                          shrinkage=0.001, bag_fraction=0.5, id=seq(nrow(data)), 
                          num_train=N/2, num_features=3)

# Fit gbm 
fit_breslow <- gbmt(Surv(tt, delta)~X1+X2+X3, data=data, distribution=cox_breslow, 
                    weights=w, train_params=params, var_monotone=c(0, 0, 0), 
                    keep_gbm_data=TRUE, cv_folds=5, is_verbose = FALSE)

fit_efron <- gbmt(Surv(tt, delta)~X1+X2+X3, data=data, distribution=cox_efron,
                  weights=w, train_params=params, var_monotone=c(0, 0, 0), 
                  keep_gbm_data=TRUE, cv_folds=5, is_verbose = FALSE)


# Evaluate fit 
plot(gbmt_performance(fit_breslow, method='test'))
legend("topleft", c("training error", "test error", "optimal iteration"),
       lty=c(1, 1, 2), col=c("black", "red", "blue"))
plot(gbmt_performance(fit_efron, method='test'))
legend("topleft", c("training error", "test error", "optimal iteration"),
       lty=c(1, 1, 2), col=c("black", "red", "blue"))
```


# Description of the underlying algorithm - specifically for CoxPH
The `gbm` algorithm works to estimate, via tree boosting, the function $f(\textbf{x})$, which maps covariates ($\textbf{x}$) to the response variable *y* - in this case the event indicator.
For `CoxPH`, the algorithm calculates both the partial log likelihood and martingale residuals ($\textbf{m}$) using the following approach.  The algorithm walks backwards in time until it encounters the "stop" time of an observation.  When this happens the weighted risk associated with that observation, $\omega_i e^{f(\textbf{x}_i)}$, is added to the total cumulative hazard: $S = \sum \omega_j e^{f(\textbf{x}_j)}$, which is initialized at $0$.   Continuing backwards in time when we reach a time before an observation was in the study, that is the algorithm leaves the associated time segment (start, stop], the observation's contribution to the cumulative hazard is subtracted off.  The algorithm is robust to overflow/underflows occurring in $e^{f(\textbf{x}_i)}$ by subtracting a constant off of the risk score.  This constant drifts to ensure overflow does not occur.

This algorithm deals with tied event times using either the Breslow or Efron approximations.  The method used is specified by the user but in the event of tied deaths, it defaults to the Efron approximation.  It also allows for the introduction of strata and start as well as stop times for each observation, see the previous Sections.  

As well as calculating the partial log likelihood the algorithm also calculates the martingale residuals.  The risk scores are related to the covariate matrix, $\mathbb{X}$, via:
$$ f(\textbf{x}_i) = (\mathbb{X}\boldsymbol{\beta})_i.   \qquad (1) $$
The derivative of the partial log likelihood, $l(\boldsymbol{\beta})$, with respect to the parameter vector $\boldsymbol{\beta}$ is related to the martingale residuals through:
$$\frac{\partial}{\partial \boldsymbol{\beta}} l(\boldsymbol{\beta}) = \mathbb{X}^{T} \textbf{m}.  \qquad (2) $$
Defining the loss function as the negative of the partial log likelihood then using the chain rule in combination with Equation (1) the residuals are given by:
$$ z_i = -\frac{\partial}{\partial f(\textbf{x}_{\textit{i}})}\Psi(\textit{y}_{\textit{i}},f(\textbf{x}_\textit{i})) = (\mathbb{X}\mathbb{X}^{T}\textbf{m})_i.  \qquad (3)$$

At this point the covariate matrix should only decide what splits the tree will make thus covariate matrix in Equation (3) is free to be set to the identity matrix and so:
$$ z_i = \textbf{m}_i. \qquad (4)$$

Finally, the updated implementation calculates the optimal terminal node predictions in the following way. Looping over the bagged observations  in the terminal node of interest the expected number of events is given by: $\sum_i \max(0.0, y_i - \textbf{m}_i) + 1/c$.  The constant $c$ represents the prior on the baseline number of events that occur within a given terminal node; it can be set on construction of the `CoxPHGBMDist` through the `prior_node_coeff_var` parameter.  From this the terminal node prediction is given by:
$$ \log(\frac{\sum_i y_i + 1/c}{\sum_i \max(0.0, y_i - \textbf{m}_i) + 1/c}). \qquad (5)$$

If this `prior_node_coeff_var` is set incorrectly, i. e. to not a finite positive double, then the predictions of the fitted model can become non-sensical.
