Survival-grafer i R

library(survival) # här finns alla survivalfunktioner
library(survminer) # här finns funktioner för grafer

# Ladda exempeldata (koloncancer)
data(colon)

# Inspektera data
# View(colon)

Skapa en Cox-regression och kalla objektet för “fit1”

fit1 <- coxph(Surv(time, status) ~ age + sex + rx, data=colon)

Skapa en Cox-justerad Kaplan-Meier (efter Kalbfleish & Prentice; man prediktar alltså överlevnaden vid varje given tidpunkt)

ggadjustedcurves(fit1, variable="rx", data = colon)

Få ut koefficienterna på forestgraf

ggforest(fit1)

Inspektera Schonfelds residualer och erhåll globalt p-värde för PH (detta kan man ifrågasätta såklart)

ggcoxzph(cox.zph(fit1))

Skapa ett överlevnadsobjekt

Vi undersöker överlevnaden i relation till cancerbehandling (“rx”)

fit2 <- survfit(Surv(time, status) ~ rx, data = colon)

Kaplan-Meier

Enkel Kaplan-Meier

ggsurvplot(fit2, data = colon, pval = TRUE)

Mer invecklad Kaplan-Meier

ggsurvplot(fit2, data = colon,
           conf.int = TRUE,
           pval = TRUE,
           fun = "pct",
           risk.table = TRUE,
           size = 1,
           linetype = "strata",
           palette = c("#E7B800",
                       "#2E9FDF",
                       "#2E3FDF"),
           legend = "bottom",
           legend.title = "Exponering",
           legend.labs = c("Hej 1",
                           "Hejsan 2",
                           "Halloj 3"))

Kumulativ incidens

ggsurvplot(fit2, data = colon, fun = "event")

Kumulativ hazard

ggsurvplot(fit2, data = colon, fun = "cumhaz")

Landmark analysis - inkluderar endast personer som levt >500 dagar

library(dplyr)
fit3 <- survfit(Surv(time, status) ~ rx, data = filter(colon, time>500))
ggsurvplot(fit3, data = colon, fun = "event")

Mer invecklad landmark analysis

# Survival första 500 dagarna
km1        <- summary(survfit(Surv(time, status) ~ rx, data = colon), times = 0:500)
km1x       <- data.frame(km1[["time"]], km1[["surv"]])
n          <- nrow(km1x)/3 # siffran 3 för att det finns 3 behandlingsnivåer i data
km1x$grupp <- c(rep("Obs", n), rep("Lev", n), rep("Lev+5FU", n))

# Survival från dag 501
km2        <- summary(survfit(Surv(time, status) ~ rx, data = filter(colon, time>500)), times = 501:max(colon$time))
km2x       <- data.frame(km2[["time"]], km2[["surv"]])
n          <- nrow(km2x)/3 # siffran 3 för att det finns 3 behandlingsnivåer i data
km2x$grupp <- c(rep("Obs", n), rep("Lev", n), rep("Lev+5FU", n))

# Ändrar namn på kolumnerna
names(km1x) <- c("Time", "Survival", "Group")
names(km2x) <- c("Time", "Survival", "Group")

# Grafen
ggplot() +
  geom_line(aes(x=km1x$Time, y=km1x$Survival, color=km1x$Group)) +
  geom_line(aes(x=km2x$Time, y=km2x$Survival, color=km2x$Group)) +
  theme_minimal() + xlab("Time (days)") + ylab("Survival probability") + ggtitle("Titel på grafen") +
  geom_vline(lty=2, aes(xintercept=500)) +
  labs(color = "Treatment")