CIF_specification.Rmd

---
title: "Approaches based on CIF specification"
author: "Karla Monterrubio-Gómez, Nathan Constantine-Cooke, and Catalina Vallejos"
date: "`r Sys.Date()`"
output:
  html_document:
    code_folding: show
    toc: yes
    css: style.css
    theme: simplex
    highlight: textmate
    toc_float:
      collapsed: no
    includes:
      in_header: head.html
      before_body: navbar.html
    self_contained: false
link-citations: yes
bibliography: references.bib
vignette: >
  %\VignetteIndexEntry{Vignette Title}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
    cache = FALSE
)
```


# Introduction

The dataset used here corresponds to the Hodgkin’s disease (HD) study described
in Pintilie, 2006. The data was loaded and pre-processed using the script that
is sourced below. The latter creates three data objects:

- `hd`: all the samples in the original data
- `hd_train`: a random sample of the data to be used as a training set
- `hd_test`: a random sample of the data to be used as a testing set

Details about data pre-processing steps can be found [here](https://github.com/VallejosGroup/CompRisksVignettes/blob/main/Source/DataPreparation.Rmd)

```{r}
source("../Data_prep/data_prep.R")
```

The following sections illustrate the usage of different methods defined using 
a CIF formulation.

# Fine-Gray 

In contrast to models that use a cause-specific CPH specification, this approach
permits inferences for the covariate effects on both the hazard function
and the survival function due to the fact that competing events are treated
differently.

As with cause-specific CPH models, there are several R packages to fit a
sub-distribution hazards model. Here we explore `cmprsk` and `riskRegression`.

## `cmprsk`

In order to fit the model, one should first create a design matrix containing
the covariates of interest, this can be done employing the `model.matrix()`
function, which creates dummy variables for the discrete predictors.

```{r, message=FALSE}
library(cmprsk)
predictors <- model.matrix(~ age +
                             sex +
                             trtgiven +
                             medwidsi +
                             extranod +
                             clinstg,
                           data = hd_train)
head(predictors)
```

Note that our covariates now correspond to binary indicators (1 or 0) of a 
specific level of the original factor variables. In order to use this design 
matrix, we will first drop the intercept. Once that has been done, we can employ 
the `crr()` function to fit the model.

```{r}
predictors <-  predictors[, -1] # drop first column (intercept)


FG_crr <- crr(ftime = hd_train$time,
              fstatus = hd_train$status,
              cov1 = predictors, # for time interactions use argument cov2
              failcode = 1,
              cencode = 0)
summary(FG_crr)
```

The above output includes the estimated covariate effects `exp(coef)`, as
well as significance level from which we can infer that the marginal effects for
age, treatment given, and clinical stage are significant. Note however that 
estimated hazard ratios are non-identical to those inferred in the 
cause-specific vignette. This is because the risks set differ between these two 
hazard functions.

In this setting, we can make inferences on both, the covariate effects on the
sub-distribution hazard and the directional effects of covariates on the CIF. 
For instance, the sub-distribution hazard ratio of relapse is 
`r round(exp(FG_crr$coef["clinstg2"]),3)` when a
patient is in clinical stage 2. In addition, because the hazard ratio is greater
than one, one can also infer that being in clinical stage 2 increases the
incidence of relapse. Here, is important to highlight that we cannot use the
magnitude of the hazard ratios to make inferences on the magnitude of the effect
on the probability of occurrence.

Because the sub-distribution hazard is directly related to the cumulative
incidence function, we can make risk predictions. In order to do this, we first
convert our test data (`hd_test`) into a design matrix format. Then, the function
`predict()` will return the estimated CIF at all unique time events for all
patients in the test set. This function computes $1-\exp(-H)$, where $H$ is the
cumulative hazard obtained using the Breslow estimator.

```{r}
predictors_test <- model.matrix(~ age +
                                  sex +
                                  trtgiven +
                                  medwidsi +
                                  extranod +
                                  clinstg,
                                data = hd_test)[, -1]

FG_CIF <- predict(object = FG_crr, cov1 = predictors_test)

dim(FG_CIF)
```

The first column in CIF_hat corresponds to the `r nrow(FG_CIF)` unique event 
times and the remaining `r nrow(FG_CIF)-1` columns are the estimated CIFs for 
each patient (given their covariate values) at such event times. As an 
illustration, the following code is used to plot the estimated CIF for events 
type 1 for two patients in the training dataset. The first patient is 37 years 
old, male, under radiation treatment, no mediastinum involvement, no extranodal 
disease, and in clinical stage 1. The second patient is 41 years old, male, 
under radiation, no mediastinum involvement, no extranodal disease, and clinical 
stage 2. 

```{r, fig.dim=c(6,4)}
par(mar = c(4, 4, 2, 0.1))
plot(FG_CIF[, 1],
     FG_CIF[, 2],
     type = "l",
     col = "red",
     ylim = c(0, 0.6),
     xlab = "Time (years)",
     ylab = "Cumulative incidence")
lines(FG_CIF[, 1], FG_CIF[, 4], col="blue")
legend("topright", legend = c("Patient 1", "Patient 2"),
       lty = c(1,1), col = c("red", "blue"))
```

> **_NOTE:_** Using `cmprsk` it is not straightforward to do predictions at a 
specific time point beyond the unique time-points present in the training data. 

## `riskRegression`

The function `FGR()` in this R package acts as a wrapper of the `crr()` function 
described above. However, it permits directly specifying the model through the 
formula argument, where we use `Hist()` to obtain the the response variable. 
When the cause argument is not passed, the default event of interest is 1.

Below we confirm that estimated coefficients are identical to those obtained
with `cmprsk` package and the interpretation is the same.

```{r, message=FALSE}
library(prodlim)         # for Hist()
library(riskRegression)

FG_riskRegression <- FGR(formula = Hist(time, status) ~
                           age +
                           sex +
                           trtgiven +
                           medwidsi +
                           extranod +
                           clinstg,
                          data = hd_train,
                          cause = 1)
print(FG_riskRegression)
```

Differently to `cmprsk`,`riskRegression` permits to do predictions of the CIF at
a specific time point (e.g. $t=5$ years) for a new dataset (e.g. `hd_test`). 
Below, we show results for the first 5 subjects in the test set.

```{r}
FG.pred <- predict(FG_riskRegression, newdata = hd_test,  times = 5)
FG.pred[c(1:5)]
```

# Sparse regression

**CoxBoost was deleted from CRAN**

# Pseudo-values

> **_NOTE:_** the package have been updated with respect to the code shown in
@klein2008. The user is referred to the package documentation for the
function.**

We define a time point grid of $M=7$ equally spaced points at which we will
compute the pseudo-values. Pseudo-observations are computed with the
`pseudoci()` function for both competing events. Below, we only consider
relapse (cause 1).

```{r, message=FALSE}
library(pseudo)

cutoffs <- c(5, 10, 15, 20, 25, 30, 35)

pseudo <- pseudoci(hd_train$time,
                   hd_train$status,
                   tmax = cutoffs)$pseudo$cause1     # we only consider cause 1
dim(pseudo)
```

`pseudo` contains the pseudo-values of the `r nrow(pseudo)` patients (rows) at 
the `r ncol(pseudo)` different time points (columns) for relapse. 

The following code will merge the pseudo-values with our train dataset.

```{r}
merged_dataset <- NULL
for (j in seq_len(ncol(pseudo))) {
    merged_dataset <- rbind(merged_dataset,
                            cbind(hd_train,
                                  pseudo = pseudo[, j],
                                  tpseudo = cutoffs[j],
                                  id = seq_len(nrow(hd_train))))
}

merged_dataset <- merged_dataset[order(merged_dataset$id), ]
merged_dataset$tpseudo <- as.factor(merged_dataset$tpseudo)
```

We will now fit a generalised linear model using `geese()` function from R
package `geepack`. We set `mean.link = "cloglog"` to use a complementary
log-log link in order to recover Fine-Gray approach (alternative link functions
can also be used, e.g. @klein2008 explored a logistic link function). 

```{r, message=FALSE}
library(geepack)

fit_pseudo <- geese(pseudo ~
                      tpseudo +
                      age +
                      sex +
                      trtgiven +
                      medwidsi +
                      extranod +
                      clinstg - 1,
                    data = merged_dataset,
                    id = id,
                    jack = TRUE,
                    scale.fix = TRUE,
                    family = gaussian,
                    mean.link = "cloglog",
                    corstr = "independence")
summary(fit_pseudo)
```

Covariate effects are shown below and are interpreted as logarithms of
sub-distribution hazard ratios. Because we used a complementary log-log link 
function, we observe that the estimated effects are akin to those obtained in 
the previous section. 

```{r}
exp(fit_pseudo$beta)
```

The estimated coefficients for the selected time points can be used to estimate 
the CIF, but this is only possible at the time-points used in the `cutoffs` grid. 

# Direct binomial

Binomial modelling approach can be fitted using the `timereg` R package.

To fit the model, we should first use the `Event()` function to obtain the
response variable to be used in the formula argument of the `comp.risk()`
function. The cause of interest must be explicitly defined through the cause
argument. The function supports several models. Here we focus on recovering the
Fine-Gray approach.

```{r, warning=FALSE}
library(timereg)

fit_DB <- comp.risk(Event(time, status) ~
                      const(age) + # const for time-invariant
                      const(sex) +
                      const(trtgiven) + 
                      const(medwidsi) +
                      const(extranod) +
                      const(clinstg),
                    data = hd_train,
                    cause = 1,
                    n.sim = 5000,
                    model = "prop",       # Fine-Gray
                    resample.iid = 1)
summary(fit_DB)
```

According to the summary output, the covariate effects associated to age, 
treatment and clinical stage are significant. Below, we show estimated hazard ratios.

```{r}
knitr::kable(round(exp(fit_DB$gamma), 4), col.names = c("Exp coef"))
```


In order to do predictions of the CIF, we use the `predict()` function. The
default makes predictions at all observed time points. Below we show estimated
incidence for the first (red) and second (blue) patients in the tests set:

```{r fig.dim=c(6,4)}
par(mar = c(4, 4, 2, 0.1))
pred_DB <- predict.timereg(fit_DB, 
                           newdata = as.data.frame(hd_test[1:2,1:6]), se = 0)
plot(pred_DB,
     multiple = 1,
     col = c("red", "blue"),
     ylim = c(0, 0.6),
     ylab = "Cumulative incidence")
legend("topright", legend = c("Patient 1", "Patient 2"),
       lty = c(1,1), col = c("red", "blue"))
```

In addition, one can also makes predictions at a specific time through the
argument `times`. Absolute risk predictions at time $t=5$ for the first 5
patients are shown:

```{r}
pred_DB.t <- predict.timereg(fit_DB, newdata = hd_test[, 1:6], times = 5)
pred_DB.t$P1[1:5]
```

`timereg` permits fitting the model time varying effects. For example, the 
following code includes a time-dependent effect for `trtgiven`.

```{r}
fit_DBtime <- comp.risk(Event(time, status) ~
                          const(age) +
                          const(sex) +
                          trtgiven + 
                          const(medwidsi) +
                          const(extranod) +
                          const(clinstg),
                        data = hd_train,
                        cause = 1,
                        model = "prop", # Fine-Gray
                        resample.iid = 1,
                        n.sim = 100)
summary(fit_DBtime)
```

# Dependent Dirichlet processes

> **_NOTE:_** `DPWeibull` is not available on CRAN because a dependency,
`binaryLogic`, was removed from CRAN as the maintainer's email address is no
longer accessible**

The formula definition uses `Hist()`.

An indication of a competing risk setting must be passed by setting
`comp = TRUE`. Below, we run the MCMC chain for 2,500 iterations, we discarded
500 as part of the burn-in period, and to reduce autocorrelation, a thinning of
10 was used (default). These values were chosen for illustration purposes only.
More iterations may be required for the algorithm to converge (the default in 
the `DPWeibull` package is to burn-in the first 8,000 iterations).


```{r, warning=FALSE , results='hide'}
library(DPWeibull)
library(prodlim) # for Hist() function
set.seed(1)
fit_DDP <- dpweib(Hist(time, status) ~
                    age +
                    sex +
                    trtgiven +
                    medwidsi +
                    extranod +
                    clinstg ,
                  data = hd_train,
                  comp = TRUE,
                  burnin = 500,
                  iteration = 2000,
                  predtime = seq(from = min(hd_train$time),
                                 to = max(hd_train$time), length =  40))
```
Note that in order to make valid inferences one should ensure that convergence 
was achieved. One can explore trace plots of the estimated parameters using the 
`coda` library. This will help us to determine if the MCMC parameters employed 
when fitting the model should be modified. For instance, below we show trace 
plots of $\alpha$ and $\lambda$ for patient 1, where we observe that the chains 
for $\alpha$ do not appear to have converged and therefore a larger number of 
MCMC iterations should be used.

```{r, message=FALSE, fig.dim=c(5,2)}
library(coda)
par(mar = c(4, 4, 2, 0.1))
traceplot(as.mcmc(fit_DDP$alpharec1[, 1]))
traceplot(as.mcmc(fit_DDP$lambda0rec1[, 1]))
```
The summary below shows estimates for the log-hazard ratio at specific time-points:

```{r}
summary(fit_DDP)
```

The following plots show the estimated (non-linear) covariate effects on the
scale of the log-subdistribution hazard ratio.

```{r}
plot(fit_DDP)
```

The latter plots can also be manually obtained as follows (here we add a 
horizontal line at zero as a visual reference):

```{r}
plot(fit_DDP$predtime, fit_DDP$loghr.est[1,], type = "l",
     main = "Log subdistribtion HR for covariate age, event 1",
     ylim = c(min(fit_DDP$loghrl[1,]), max(fit_DDP$loghru[1,])))
lines(fit_DDP$predtime, fit_DDP$loghrl[1,], lty = 2)
lines(fit_DDP$predtime, fit_DDP$loghru[1,], lty = 2)
abline(h = 0, col = "red")

plot(fit_DDP$predtime, fit_DDP$loghr.est[2,], type = "l",
     main = "Log subdistribtion HR for covariate sexM, event 1",
     ylim = c(min(fit_DDP$loghrl[2,]), max(fit_DDP$loghru[2,])))
lines(fit_DDP$predtime, fit_DDP$loghrl[2,], lty = 2)
lines(fit_DDP$predtime, fit_DDP$loghru[2,], lty = 2)
abline(h = 0, col = "red")
```

In the example above, age appears to have a diminishing effect over time. Note 
also that, in all cases, there seems to be high posterior uncertainty and that 
reported credible intervals contain zero.  

Below, we illustrate how to do predictions with `DPWeibull`. First, we show how
to estimate the CIF for two subjects in the test set at all observed event
times.

```{r fig.dim=c(6,4)}
pred_DDP <- predict(fit_DDP, newdata = as.data.frame(hd_test[1:2, 1:6]))
```

Note that `pred_DDP` contains posterior samples for all the two subjects at 
`r ncol(pred_DDP$tpred)` observed event times, as well as posterior estimated 
and credible intervals.

Here we show how to manually extract and plot the estimated CIF curves for 
patient 1 (red) and 2 (blue). Dashed lines show the associated credible intervals.

```{r fig.dim=c(6,4)}
plot(pred_DDP$tpred,
     pred_DDP$Fpred[1, ],
     type = "l",
     col = "red",
     xlab = "Time (years)",
     ylab = "Cumulative incidence",
     ylim = c(0, 0.6))
lines(pred_DDP$tpred, pred_DDP$Fpredl[1,], col = "red", lty = 2)
lines(pred_DDP$tpred, pred_DDP$Fpredu[1,], col = "red", lty = 2)
lines(pred_DDP$tpred, pred_DDP$Fpred[2, ], col="blue")
lines(pred_DDP$tpred, pred_DDP$Fpredl[2,], col = "blue", lty = 2)
lines(pred_DDP$tpred, pred_DDP$Fpredu[2,], col = "blue", lty = 2)
legend("topright", legend = c("Patient 1", "Patient 2"),
       lty = c(1,1), col = c("red", "blue"))
```

Note that `DPWeibull` also provides a `plot` method that can generate similar 
graphical summaries. However, the previous code may be preferrable when 
generating predictions for a large number of individuals.

```{r}
plot(pred_DDP)
```

In addition, we can also make predictions at a specific time point, e.g.
$t=5$. Below, we show the estimated risk for the first 5 subjects in the test
data set:

```{r}
pred_DDP.t <- predict(fit_DDP, 
                      newdata = as.data.frame(hd_test[, 1:6]), tpred = 5)
pred_DDP.t$Fpred[1:5]
```



# Comparison: estimates of regression coefficients

First, we compare estimated regression coefficients obtained above. Note that 
this is only for (semi-)parametric approaches that follow the Fine & Gray 
specification, with an linear effect in the log-scale of the sub-distribution
hazard. As the estimates obtained for the Fine & Gray were the same regardless 
of the choice of R library (`cmprsk` o `riskRegression`), we only store one
of them (the ones for `cmprsk`).

```{r}
my.vars <- rownames(summary(FG_crr)$coef)

# X1-3 is added so that the methods are plotted in the desired order
df <- data.frame("X1.fg_est" = summary(FG_crr)$coef[,1],
                 "X1.fg_se" = summary(FG_crr)$coef[,3],
                 "X2.pseudo_est" = summary(fit_pseudo)$mean[my.vars,1],
                 "X2.pseudo_se" = summary(fit_pseudo)$mean[my.vars,2],
                 "X3.db_est" = coef(fit_DB)[,1],
                 "X3.db_se" = coef(fit_DB)[,2])
round(df, 2)
df$variable <- rownames(df)

library(tidyverse)
df_long <- df %>%
  pivot_longer(!variable, 
               names_to = c("Method", ".value"), 
               names_sep = "_")
```

The following code was used to create Figure 1B:

```{r}
# coloured version
library(ggplot2)
library(patchwork)

p_all <- ggplot(df_long, aes(x = variable, y = est, col = Method)) +
  geom_point(cex = 2.5, position = position_dodge(width=0.4)) +
  geom_errorbar(aes(ymin = est - se, ymax = est + se), 
                width=.1, position = position_dodge(width=0.4)) +
  geom_hline(yintercept = 0, lty = 2, colour = "gray") +
  theme_classic() +
  ylab("Estimate") + xlab("Coefficient") +
  scale_color_manual(values = c("#999999", "#1A5276", "#B9770E"),
                     labels = c('Fine & Gray', 'Pseudo values', 'Direct binomial')) 

p_age <- ggplot(df_long[df_long$variable == "age",], 
                aes(x = variable, y = est, col = Method)) +
  geom_point(cex = 2.5, position = position_dodge(width=0.6)) +
  geom_errorbar(aes(ymin = est - se, ymax = est + se), 
                width=.1, position = position_dodge(width=0.6)) +
  theme_classic() +
  ylab("Estimate") + xlab("Coefficient") +
  scale_color_manual(values = c("#999999", "#1A5276", "#B9770E"),
                     labels = c('Fine & Gray', 'Pseudo values', 'Direct binomial')) +
  theme(legend.position="none", 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        plot.background = element_rect(colour = "black", fill=NA, size=0.6))

p_all + inset_element(p_age, left = 0.004, bottom = 0.6, right = 0.19, top = 1)

if (file.exists("/.dockerenv")){ # running in docker
  ggsave("/Outputs/Comparison_estimation_CIF.pdf", device = "pdf")
} else {
  ggsave("../Outputs/Comparison_estimation_CIF.pdf", device = "pdf")
}
```

```{r}
# black and white

p_all <- ggplot(df_long, aes(x = variable, y = est, pch = Method)) +
  geom_point(cex = 2.5, position = position_dodge(width=0.4)) +
  geom_errorbar(aes(ymin = est - se, ymax = est + se), 
                width=.1, position = position_dodge(width=0.4)) +
  geom_hline(yintercept = 0, lty = 2, colour = "gray") +
  theme_classic() +
  ylab("Estimate") + xlab("Coefficient") +
  scale_shape_manual(values = c(16, 15, 17), 
                     labels = c('Fine & Gray', 'Pseudo values', 'Direct binomial')) 

p_age <- ggplot(df_long[df_long$variable == "age",], 
                aes(x = variable, y = est, pch = Method)) +
  geom_point(cex = 2.5, position = position_dodge(width=0.6)) +
  geom_errorbar(aes(ymin = est - se, ymax = est + se), 
                width=.1, position = position_dodge(width=0.6)) +
  theme_classic() +
  ylab("Estimate") + xlab("Coefficient") +
  scale_shape_manual(values = c(16, 15, 17), 
                     labels = c('Fine & Gray', 'Pseudo values', 'Direct binomial'))  +
  theme(legend.position="none", 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        plot.background = element_rect(colour = "black", fill=NA, size=0.6))

p_all + inset_element(p_age, left = 0.004, bottom = 0.6, right = 0.19, top = 1)

if (file.exists("/.dockerenv")){ # running in docker
  ggsave("/Outputs/Comparison_estimation_CIF_BW.pdf", device = "pdf")
} else {
  ggsave("../Outputs/Comparison_estimation_CIF_BW.pdf", device = "pdf")
}
```

# Storing predictions

In order to allow comparison with the predictions generated by other methods,
we save the predictions obtained in this vignette.  

```{r save_pred}
pred_CIF <- data.frame("testID" = seq_len(nrow(hd_test)),
                      "FG_riskRegression" = FG.pred,
                      "DirectBinomial" = pred_DB.t$P1,
                      "DPWeibull" = pred_DDP.t$Fpred)
if (file.exists("/.dockerenv")){ # running in docker
  write.csv(pred_CIF, "/Outputs/pred_CIF.csv", row.names = FALSE)
} else {
  write.csv(pred_CIF, "../Outputs/pred_CIF.csv", row.names = FALSE)
}
```

# References

<div id="refs"></div>

# Session Info
```{r}
sessionInfo()
```