Further survival analyses
a03_Further_survival_analyses.RmdSet up
Let us first load the packages required.
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()
#> Rows: ??
#> Columns: 13
#> Database: DuckDB v1.1.3 [unknown@Linux 6.8.0-1021-azure:R 4.4.2/:memory:]
#> $ cohort_definition_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ subject_id <int> 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 14, 15, 16, 19, 2…
#> $ cohort_start_date <date> 1981-01-01, 1968-01-01, 1980-01-01, 1977-01-01, …
#> $ cohort_end_date <date> 1981-01-01, 1968-01-01, 1980-01-01, 1977-01-01, …
#> $ age <dbl> 88, 78, 94, 68, 90, 90, 89, 87, 79, 86, 80, 85, 9…
#> $ sex <fct> F, F, M, M, F, M, F, F, F, M, F, M, F, M, M, F, F…
#> $ hgb <dbl> 13.1, 11.5, 10.5, 15.2, 10.7, 12.9, 10.5, 12.3, 9…
#> $ creat <dbl> 1.3, 1.2, 1.5, 1.2, 0.8, 1.0, 0.9, 1.2, 1.1, 1.0,…
#> $ mspike <dbl> 0.5, 2.0, 2.6, 1.2, 1.0, 0.5, 1.3, 1.6, 2.3, 2.3,…
#> $ age_group <chr> ">=70", ">=70", ">=70", "<70", ">=70", ">=70", ">…
#> $ days_to_exit <int> 30, 25, 46, 92, 8, 4, 151, 2, 136, 2, 14, 18, 43,…
#> $ status <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ time <dbl> 30, 25, 46, 92, 8, 4, 151, 2, 136, 2, 14, 18, 43,…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()
#> Rows: ??
#> Columns: 13
#> Database: DuckDB v1.1.3 [unknown@Linux 6.8.0-1021-azure:R 4.4.2/:memory:]
#> $ cohort_definition_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ subject_id <int> 1340, 1376, 1380, 1382, 1248, 213, 1245, 1374, 67…
#> $ cohort_start_date <date> 1994-01-01, 1994-01-01, 1994-01-01, 1994-01-01, …
#> $ cohort_end_date <date> 1994-01-01, 1994-01-01, 1994-01-01, 1994-01-01, …
#> $ age <dbl> 67, 67, 69, 66, 87, 93, 82, 68, 81, 79, 85, 86, 8…
#> $ sex <fct> F, F, M, M, M, F, M, F, M, M, M, F, F, F, F, F, M…
#> $ hgb <dbl> 12.2, 13.7, 15.0, 12.1, 12.7, 12.8, 11.4, 9.2, 13…
#> $ creat <dbl> 1.4, 1.1, 0.8, 2.0, 1.5, 1.1, 1.5, 1.8, 1.4, 1.1,…
#> $ mspike <dbl> 1.2, 1.5, 0.0, 0.0, 0.5, 0.8, 1.4, 0.5, 1.3, 1.7,…
#> $ age_group <chr> "<70", "<70", "<70", "<70", ">=70", ">=70", ">=70…
#> $ days_to_exit <dbl> 46, 41, 22, 31, 10, 19, 1, 40, 43, 6, 12, 57, 52,…
#> $ status <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0…
#> $ time <dbl> 46, 41, 22, 31, 10, 19, 1, 40, 43, 6, 12, 57, 52,…
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()
#> Rows: ??
#> Columns: 13
#> Database: DuckDB v1.1.3 [unknown@Linux 6.8.0-1021-azure:R 4.4.2/:memory:]
#> $ cohort_definition_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ subject_id <int> 1340, 1376, 1380, 1382, 1248, 213, 1245, 1374, 67…
#> $ cohort_start_date <date> 1994-01-01, 1994-01-01, 1994-01-01, 1994-01-01, …
#> $ cohort_end_date <date> 1994-01-01, 1994-01-01, 1994-01-01, 1994-01-01, …
#> $ age <dbl> 67, 67, 69, 66, 87, 93, 82, 68, 81, 79, 85, 86, 8…
#> $ sex <fct> F, F, M, M, M, F, M, F, M, M, M, F, F, F, F, F, M…
#> $ hgb <dbl> 12.2, 13.7, 15.0, 12.1, 12.7, 12.8, 11.4, 9.2, 13…
#> $ creat <dbl> 1.4, 1.1, 0.8, 2.0, 1.5, 1.1, 1.5, 1.8, 1.4, 1.1,…
#> $ mspike <dbl> 1.2, 1.5, 0.0, 0.0, 0.5, 0.8, 1.4, 0.5, 1.3, 1.7,…
#> $ age_group <chr> "<70", "<70", "<70", "<70", ">=70", ">=70", ">=70…
#> $ days_to_exit <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ status <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
#> $ time <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…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)
#> Call:
#> survival::coxph(formula = survival::Surv(time, status) ~ age +
#> sex, data = input_survival_single)
#>
#> coef exp(coef) se(coef) z p
#> age 0.061622 1.063561 0.003402 18.114 < 2e-16
#> sexM 0.358258 1.430835 0.065693 5.454 4.94e-08
#>
#> Likelihood ratio test=391.2 on 2 df, p=< 2.2e-16
#> n= 1384, number of events= 963
survival::survdiff(survival::Surv(time, status) ~ sex, data = input_survival_single)
#> Call:
#> survival::survdiff(formula = survival::Surv(time, status) ~ sex,
#> data = input_survival_single)
#>
#> N Observed Expected (O-E)^2/E (O-E)^2/V
#> sex=F 631 423 471 4.88 9.67
#> sex=M 753 540 492 4.67 9.67
#>
#> Chisq= 9.7 on 1 degrees of freedom, p= 0.002Further analysis with competing risk survival
For competing risk settings, there is a similar function that adds
time and status information to the cohort of
interest. We only need to specify which are the outcome and competing
outcome of interest. We can also choose other options such as varying
the follow up time, censoring on a specific date, washout periods, or
others. You can check all the options this function offers by running
?addCompetingRiskCohortSurvival. With the default
parameters:
input_survival_cr <- cdm$mgus_diagnosis %>%
addCompetingRiskCohortSurvival(
cdm = cdm,
outcomeCohortTable = "progression",
outcomeCohortId = 1,
competingOutcomeCohortTable = "death_cohort",
competingOutcomeCohortId = 1
) %>%
glimpse()
#> Rows: 1,384
#> Columns: 13
#> $ cohort_definition_id <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ subject_id <int> 56, 81, 124, 127, 147, 163, 165, 180, 186, 190, 1…
#> $ cohort_start_date <date> 1978-01-01, 1985-01-01, 1974-01-01, 1978-01-01, …
#> $ cohort_end_date <date> 1978-01-01, 1985-01-01, 1974-01-01, 1978-01-01, …
#> $ age <dbl> 78, 91, 73, 73, 58, 57, 80, 70, 76, 78, 54, 79, 6…
#> $ sex <fct> M, F, M, M, M, F, M, F, F, M, F, M, M, M, F, M, F…
#> $ hgb <dbl> 10.3, 5.9, 15.3, 12.4, 13.1, 12.2, 11.0, 14.3, 12…
#> $ creat <dbl> 3.0, 1.0, 1.2, 1.6, 1.1, 0.8, 1.4, 1.2, 0.9, 1.1,…
#> $ mspike <dbl> 1.9, 0.0, 1.7, 1.4, 0.8, 1.9, 2.0, 1.6, 1.9, 0.9,…
#> $ age_group <chr> ">=70", ">=70", ">=70", ">=70", "<70", "<70", ">=…
#> $ days_to_exit <int> 44, 21, 82, 60, 189, 260, 85, 107, 104, 101, 171,…
#> $ time <dbl> 29, 14, 80, 30, 188, 201, 76, 81, 51, 101, 153, 8…
#> $ status <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1…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))
#> Competing Risks Regression
#>
#> Call:
#> cmprsk::crr(ftime = input_survival_cr$time, fstatus = input_survival_cr$status,
#> cov1 = covs, failcode = 1, cencode = 0)
#>
#> coef exp(coef) se(coef) z p-value
#> age -0.0192 0.981 0.00585 -3.28 0.001
#> sex 0.2871 1.333 0.19309 1.49 0.140
#>
#> exp(coef) exp(-coef) 2.5% 97.5%
#> age 0.981 1.02 0.970 0.992
#> sex 1.333 0.75 0.913 1.945
#>
#> Num. cases = 1384
#> Pseudo Log-likelihood = -726
#> Pseudo likelihood ratio test = 8.32 on 2 df,