Further survival analyses
In CohortSurvival: Estimate Survival from Common Data Model Cohorts

knitr::opts_chunk$set(
  collapse = TRUE, warning = FALSE, message = FALSE,
  comment = "#>"
)

Set up

Let us first load the packages required.

library(CDMConnector)
library(CohortSurvival)
library(dplyr)
library(cmprsk)
library(survival)

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

cdm <- CohortSurvival::mockMGUS2cdm()

The CohortSurvival package does not have implemented functionality to do more complex survival analyses than Kaplar Meier curves, like Cox Proportional Hazards modelling. However, the format the data has to be in to be inputted to well-known modelling functions from packages like survival or cmprskcan be retrieved from OMOP data with some in-built functions in this package. 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.

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 ?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 their outcome variable and censor them at a particular date, for instance, the 1st of January of 1994. We see how that gives us different results:

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:

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 need to add time and status information for the outcome, then for the competing outcome. Then we leverage all those variables to get what outcome (if any) to count per individual so that we can feed the result to subsequent models.

# 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 
  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:

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))