SURVIVAL ANALIZA
dr.sc.Luka Sikic
27/04/2020
Survival analiza
Survival analiza je skup statističkih pristupa koji se koriste za istraživanje vremena potrebnog za realizaciju događaja od interesa. Primjenu ima u nizu disciplina poput biologije, sociologije, ekonomije, inžinjeringu proizvoda, a u literaturi se sreće pod više različithih naziva: duracijska analiza, time-event analiza, analiza vremena do otkaza proizvoda i sl. Ključni elementi survival analize su: vrijeme i događaj, cenzori, survival i hazard funkcije.
Vrijeme se odnosi na period do nastupanja događaja koji se analizira. Događaj se može odnositi na nastupanje smrti, zaposlenja, otkaza proizvoda, promjene političkog režima i dr.
Važan aspekt survival analize je pravilan tretman cenzoriranih opservacija. Cenzorirane opservacije se odnose na potrebu za prilagodbom analiziranih podataka usljed događaja koji nisu nastupili unutar istraživanog vremenskog perioda. Razlozi za ne-nastupanje događaja mogu biti nastupanje događaja nakon analizom razmatranog perioda, gubitak analiziranog objekta u razdoblju analize, nastupanje drugog događaja koji je u konflikutu sa istraživanim.
Statistička vjerojatnost u survival analizi se procijenjuje survival funkcijom i hazard funkcijom. Survival funkcija opisuje vjerojanost da će objekt istraživanja preživjeti u razdoblju od početka do analizom određenog trenutka nastupanja događaja. Log-rank test se može koristiti za analizu razlika između survival krivulja za analizirane grupe. Hazard funkcija određuje vjerojatnost instantnog nastupanja događaja u određenom vremenskom trenutku. Koristi se kao dijagnostički alat ili u cilju specifikacije matematičkog modela u survival analizi.
Najčešće korištena statistička procedura za procjenu survival funkcije je Kaplan-Meier metoda(KM). KM je neparametarski statistički model za procjenu survival vjerojatnosti S(t) u definiranom vremenskom intervalu. Procijenjena survival vjerojatnost S(t) je stepenasta funkcija koja odražava promjenu vrijednosti u trenutku nastupanja događaja za svaki promatrani objekt. Na osnovi survival krivulje je moguće donijeti zaključke o prosječnom survival vremenu za objekte određenih karakteristika.
U survival analizi se najčešće koriste sljedeće metode: - Kaplan-Meier krivulje za vizualizaciju - Log-rank test za usporedbu survival krivulja kod dvije ili više grupa - Cox regresija proporcionalnog hazada za analizu utjecaja varijabli na survival
Praktični primjer Survival analize
# PODATCI
# Učitaj pakete
library("survival")
library("survminer")
# Podatci
data("lung")
head(lung)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0# inst: institucija
# time: survival vrijeme u danima
# status: cenzori 1=cenzurirano, 2=nastupanje smrti
# age: Age in years
# sex: mr=1 žr=2
# ph.ecog: ECOG pokazatelj (0=good 5=dead)
# ph.karno: Karnofsky pokazatelj (bad=0-good=100) rated by physician
# pat.karno: subjektivni Karnofsky pokazatelj
# meal.cal: kalorijski unos
# wt.loss: gubitak težine u kilogramima# PROCJENA
# Procijeni survival krivulje
fit <- survfit(Surv(time, status) ~ sex, data = lung)
# Pregled podataka
print(fit)## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## n events median 0.95LCL 0.95UCL
## sex=1 138 112 270 212 310
## sex=2 90 53 426 348 550# Pregled podataka (not run)
##summary(fit)
summary(fit)$table## records n.max n.start events *rmean *se(rmean) median 0.95LCL 0.95UCL
## sex=1 138 138 138 112 325.0663 22.59845 270 212 310
## sex=2 90 90 90 53 458.2757 33.78530 426 348 550# survfit objekti
names(fit)## [1] "n" "time" "n.risk" "n.event" "n.censor" "surv"
## [7] "std.err" "cumhaz" "std.chaz" "strata" "type" "logse"
## [13] "conf.int" "conf.type" "lower" "upper" "call"d <- data.frame(time = fit$time,
n.risk = fit$n.risk,
n.event = fit$n.event,
n.censor = fit$n.censor,
surv = fit$surv,
upper = fit$upper,
lower = fit$lower
)
head(d)## time n.risk n.event n.censor surv upper lower
## 1 11 138 3 0 0.9782609 1.0000000 0.9542301
## 2 12 135 1 0 0.9710145 0.9994124 0.9434235
## 3 13 134 2 0 0.9565217 0.9911586 0.9230952
## 4 15 132 1 0 0.9492754 0.9866017 0.9133612
## 5 26 131 1 0 0.9420290 0.9818365 0.9038355
## 6 30 130 1 0 0.9347826 0.9768989 0.8944820# VIZUALIZACIJA SURVIVAL KRIVULJA
# VIZUALIZACIJA 1
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"))
# VIZUALIZACIJA 2
ggsurvplot(
fit,
pval = TRUE,
conf.int = TRUE,
conf.int.style = "step",
xlab = "Time in days",
break.time.by = 200,
ggtheme = theme_light(),
risk.table = "abs_pct",
risk.table.y.text.col = T,
risk.table.y.text = FALSE,
ncensor.plot = TRUE,
surv.median.line = "hv",
legend.labs =
c("Male", "Female"),
palette =
c("#E7B800", "#2E9FDF")
)
# VIZUALIZACIJA 3
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
xlim = c(0, 600))
# VIZUALIZACIJA KUMULATIVNIH DOGAĐAJA I HAZARDA
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_bw(),
palette = c("#E7B800", "#2E9FDF"),
fun = "event")
# VIZUALIZACIJA KUMULATIVNOG HAZARDA
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata",
ggtheme = theme_bw(),
palette = c("#E7B800", "#2E9FDF"),
fun = "cumhaz")
# KAPLAN-MEIER TABLICA
res.sum <- surv_summary(fit)
head(res.sum)## time n.risk n.event n.censor surv std.err upper lower strata
## 1 11 138 3 0 0.9782609 0.01268978 1.0000000 0.9542301 sex=1
## 2 12 135 1 0 0.9710145 0.01470747 0.9994124 0.9434235 sex=1
## 3 13 134 2 0 0.9565217 0.01814885 0.9911586 0.9230952 sex=1
## 4 15 132 1 0 0.9492754 0.01967768 0.9866017 0.9133612 sex=1
## 5 26 131 1 0 0.9420290 0.02111708 0.9818365 0.9038355 sex=1
## 6 30 130 1 0 0.9347826 0.02248469 0.9768989 0.8944820 sex=1
## sex
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1# LOG-RANK TEST ZA USPOREDBU SURVIVAL KRIVULJA
surv_razlika <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_razlika## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001# SLOŽENIJI SURVIVAL MODEL
# Provedi model
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
data = colon )
fit2## Call: survfit(formula = Surv(time, status) ~ sex + rx + adhere, data = colon)
##
## n events median 0.95LCL 0.95UCL
## sex=0, rx=Obs , adhere=0 250 126 2213 1375 NA
## sex=0, rx=Obs , adhere=1 48 32 963 594 2077
## sex=0, rx=Lev , adhere=0 222 113 2079 1191 NA
## sex=0, rx=Lev , adhere=1 44 24 1737 997 NA
## sex=0, rx=Lev+5FU, adhere=0 282 123 NA NA NA
## sex=0, rx=Lev+5FU, adhere=1 44 26 1142 649 NA
## sex=1, rx=Obs , adhere=0 286 161 1692 1216 2288
## sex=1, rx=Obs , adhere=1 46 26 1079 760 NA
## sex=1, rx=Lev , adhere=0 300 156 1976 1135 NA
## sex=1, rx=Lev , adhere=1 54 40 1000 717 1540
## sex=1, rx=Lev+5FU, adhere=0 248 80 NA NA NA
## sex=1, rx=Lev+5FU, adhere=1 34 13 NA 2318 NA# Vizualiziraj rezultate
gg_surv_complex <- ggsurvplot(fit2,
fun = "event",
conf.int = TRUE,
ggtheme = theme_bw())
gg_surv_complex$plot +theme_bw() +
theme (legend.position = "right")+
facet_grid(rx ~ adhere)
# SLOŽENIJI SURVIVAL (COX) MODEL
# Univarijatni model
res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
res.cox## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## coef exp(coef) se(coef) z p
## sex -0.5310 0.5880 0.1672 -3.176 0.00149
##
## Likelihood ratio test=10.63 on 1 df, p=0.001111
## n= 228, number of events= 165summary(res.cox)## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.5310 0.5880 0.1672 -3.176 0.00149 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.588 1.701 0.4237 0.816
##
## Concordance= 0.579 (se = 0.021 )
## Likelihood ratio test= 10.63 on 1 df, p=0.001
## Wald test = 10.09 on 1 df, p=0.001
## Score (logrank) test = 10.33 on 1 df, p=0.001# Alternativno
covariates <- c("age", "sex", "ph.karno", "ph.ecog", "wt.loss")
univ_formulas <- sapply(covariates,
function(x) {as.formula(paste('Surv(time, status)~', x))})
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = lung)})
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value <- signif(x$wald["pvalue"], digits=2)
wald.test <- signif(x$wald["test"], digits=2)
beta <- signif(x$coef[1], digits=2) #coeff beta
HR <- signif(x$coef[2], digits=2) #exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res <- c(beta, HR, wald.test, p.value)
names(res) <-c ("beta", "HR (95% CI for HR)", "wald.test",
"p.value")
return(res)
})
res <- t(as.data.frame(univ_results, check.names = FALSE))
as.data.frame(res)## beta HR (95% CI for HR) wald.test p.value
## age 0.019 1 (1-1) 4.1 0.042
## sex -0.53 0.59 (0.42-0.82) 10 0.0015
## ph.karno -0.016 0.98 (0.97-1) 7.9 0.005
## ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05
## wt.loss 0.0013 1 (0.99-1) 0.05 0.83# Multivarijatni model
res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
summary(res.cox)## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
##
## n= 227, number of events= 164
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.011067 1.011128 0.009267 1.194 0.232416
## sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
## ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0111 0.9890 0.9929 1.0297
## sex 0.5754 1.7378 0.4142 0.7994
## ph.ecog 1.5900 0.6289 1.2727 1.9864
##
## Concordance= 0.637 (se = 0.025 )
## Likelihood ratio test= 30.5 on 3 df, p=1e-06
## Wald test = 29.93 on 3 df, p=1e-06
## Score (logrank) test = 30.5 on 3 df, p=1e-06## Procjeni efekte za spol
sex_df <- with(lung, data.frame(sex = c(1, 2),
age = rep(mean(age, na.rm = TRUE), 2),
ph.ecog = c(1, 1)
)
)
sex_df## sex age ph.ecog
## 1 1 62.44737 1
## 2 2 62.44737 1# Vizualiziraj multivarijatni model (not run)
#fit3 <- survfit(res.cox, newdata = sex_df)
#ggsurvplot(fit3, newdata = sex_df, conf.int = TRUE, legend.labs=c("Sex=1", "Sex=2"),
# ggtheme = theme_minimal())