---
title: "Further survival analyses"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Further survival analyses}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Set up

Let us first load the packages required.

```{r}
library(CDMConnector)
library(CohortSurvival)
library(dplyr)
library(cmprsk)
library(survival)
```

We will create a cdm reference to use our example MGUS2 survival dataset.

```{r}
cdm <- CohortSurvival::mockMGUS2cdm()
```

The CohortSurvival package focuses on Kaplan-Meier survival estimates and cumulative incidence in a competing-risk setting. It does not fit more complex models, such as Cox proportional hazards or Fine and Gray models, directly. However, the format the data has to be in to be inputted to well-known modelling functions from packages like `survival` or `cmprsk` can be retrieved from OMOP data with `addCohortSurvival()`. Let us see how to do it in both single-event and competing-risk survival settings.

## Further analysis with single event survival

To get the `time` and `status` information we need for the `coxph` function in the package `survival`, for instance, we only need to call `addCohortSurvival()`. The stratification variables need to be columns previously added to the cohort by the user. `status` is `1` for people with the outcome event and `0` for people censored before the event. `time` is the number of days from cohort entry to the event or censoring.

```{r, fig.width=5}
input_survival_single <- cdm$mgus_diagnosis |>
       addCohortSurvival(
       cdm = cdm,
       outcomeCohortTable = "death_cohort",
       outcomeCohortId = 1
       ) 

input_survival_single |> 
  glimpse()
```

We can decide to change some of the default parameters in this function. Information on all these can be found in `?addCohortSurvival`. For instance, we can choose to exclude people with an outcome only 180 days before index date, instead of anytime, and follow them up for only one year. We can also decide to use `cohort_end_date` as the outcome date variable and censor them at a particular date, for instance, the 1st of January of 1994. We see how that gives us different results:

```{r}
cdm$mgus_diagnosis |>
       addCohortSurvival(
       cdm = cdm,
       outcomeCohortTable = "death_cohort",
       outcomeWashout = 180,
       followUpDays = 365
       ) |>
  filter(cohort_start_date > "1993-01-01") |>
  glimpse()
cdm$mgus_diagnosis |>
       addCohortSurvival(
       cdm = cdm,
       outcomeCohortTable = "death_cohort",
       outcomeDateVariable = "cohort_end_date",
       censorOnDate = as.Date("1994-01-01")
       ) |>
    filter(cohort_start_date > "1993-01-01") |>
  glimpse()
```

This table with the added `time` and `status` information should be enough to call any advanced function, like the aforementioned Cox Proportional Hazards model:

```{r}
survival::coxph(survival::Surv(time, status) ~ age + sex, data = input_survival_single)
survival::survdiff(survival::Surv(time, status) ~ sex, data = input_survival_single)
```

## Further analysis with competing risk survival

For competing-risk settings, we need to use the same function that adds `time` and `status` information, but twice. We first add time and status information for the outcome, then for the competing outcome. Then we combine those variables to identify which outcome, if any, happened first for each individual so that we can feed the result to subsequent models. In the coding below, `status = 0` means censored, `status = 1` means the event of interest, and `status = 2` means the competing outcome.

```{r}

# Add all status and time information for both outcomes
  input_survival_cr <- cdm$mgus_diagnosis |>
    addCohortSurvival(cdm, "progression") |>
    dplyr::rename(
      "outcome_time" = "time",
      "outcome_status" = "status"
    ) |>
     addCohortSurvival(cdm, "death_cohort") |>
    dplyr::rename(
      "competing_outcome_time" = "time",
      "competing_outcome_status" = "status"
    )
  
  # Collect and combine the two event processes
  input_survival_cr <- input_survival_cr |>
    dplyr::collect() |>
    dplyr::mutate(
      time = pmin(outcome_time, competing_outcome_time),
      status = factor(
        dplyr::if_else(competing_outcome_time <= outcome_time, 2 * competing_outcome_status, outcome_status))
    ) |>
    dplyr::select(-c("outcome_time", "outcome_status", "competing_outcome_time", "competing_outcome_status"))

```

We can use the package `cmprsk` to fit a Fine and Gray model to the competing risk data. We first change our `sex` covariate to numeric, and then we can run the analysis:

```{r, fig.height=6, fig.width=8}
input_survival_cr <- input_survival_cr |>
  dplyr::mutate(sex = dplyr::if_else(sex == "M", 0, 1))

covs <- data.frame(input_survival_cr$age, input_survival_cr$sex)
names(covs) <- c("age", "sex")

summary(cmprsk::crr(ftime = input_survival_cr$time,
            fstatus = input_survival_cr$status,
            cov1 = covs,
            failcode = 1,
            cencode = 0))
```

## Disconnect from the cdm database connection

We finish by disconnecting from the cdm database connection.

```{r}
cdmDisconnect(cdm)
```