I am studying survival analysis and am trying to see if there's a way to probabilistically forecast future outcomes, using simulation or other means.
In the first example below, I fit a Cox model to the complete "lung" data from the
survival package, showing 1000 months of outcomes. In the second example, I adjust the "lung" data as-if I only had 500 months of survival data, creating object "lung1".
Using survival analysis, how could I probabilistically forecast events for months 501-1000 for lung1, assuming I only had data for months 1-500? I've used time-series forecasting models (ETS, ARIMA, etc.) but I wonder if there's a better solution using survival analysis? A problem with these time-series models is generating negative survival outcomes which obviously is impossible. Nevertheless, I post an image below of an ETS forecast model I've used before with log adjustments to eliminate negative-value outcomes.
I post simple code for the Cox survival models at the bottom. Images for "lung" and truncated "lung1" data:
Example of ETS time-series model forecast (using other data):
# Example from http://www.sthda.com/english/wiki/cox-proportional-hazards-model library(survival) library(survminer) # status 1 = censored # status 2 = dead ### Full data set ### # Cox regression of time to death on the time-constant covariates cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung) # Plot the baseline survival function ggsurvplot(survfit(cox, data = lung), palette = "#2E9FDF", ggtheme = theme_minimal()) ### Truncate the full data set "as if" we only had the first half of the time series available # lung1 reduces study time to 500 months (from 1000) and adjusts status (via status1) at month 500 cut-off lung1 <- lung %>% mutate(time1 = pmin(time,500)) %>% mutate(status1 = if_else(time > time1,as.integer(1),as.integer(status))) # Cox regression of time to death on the time-constant covariates cox1 <- coxph(Surv(time1, status1) ~ age + sex + ph.ecog, data = lung1) # Plot the truncated survival data myplot <- ggsurvplot(survfit(cox1, data = lung1), palette = "#2E9FDF", ggtheme = theme_minimal(),xlim = c(0, 1000)) myplot$plot <- myplot$plot + scale_x_continuous(breaks = sort(c(seq(0, 1000, 250)))) myplot