############################# Analyza preziti #############################

### nacteni dat:
aml <- read.table(file = "alert_data_2010.txt", sep = "\t", dec = ".", header = TRUE)

### pocet pacientu:
nrow(aml)

### pocet zemrelych:
sum(aml$OS_cenz==1)

### graf Kaplan-Meierova (K-M) odhadu pro celkove preziti (definovano predem jako datum posledniho sledovani/umrti - datum dg):
# OS_m: cas od data dg do data posledniho sledovani/umrti 
# OS_cenz: 1..zemreli pacienti, 0..cenzorovani pacienti
library(survival) # nacteni knihovny "survival"
km.surv <- survfit(Surv(OS_m, OS_cenz) ~ 1, data = aml) # vypocet K-M odhadu
plot(km.surv, xlab="Mesice", ylab="Podil zijicich pacientu") # vykresleni K-M krivky preziti
summary(km.surv) # vypis jednotlivych casu udalosti a poctu udalosti

### K-M krivka pro vice skupin pacientu - celkove preziti (dle K-M) pro zeny a muze:
km.pohlavi <- survfit(Surv(OS_m, OS_cenz) ~ pohlavi, data = aml) # pohlavi: 1..zeny, 2..muzi
plot(km.pohlavi[pohlavi=1]) # vykresleni krivky preziti pro zeny
plot(km.pohlavi[pohlavi=2]) # vykresleni krivky preziti pro muze
plot(km.pohlavi) # vykresleni obou krivek preziti do 1 grafu (krivky jsou bez intervalu spolehlivosti)

# vykresleni obou krivek preziti do 1 grafu i s intervaly spolehlivosti:
plot(km.pohlavi[pohlavi=1])
lines(km.pohlavi[pohlavi=2]$time, km.pohlavi[pohlavi=2]$surv, type="s", col=2)
lines(km.pohlavi[pohlavi=2]$time, km.pohlavi[pohlavi=2]$lower, type="s", col=2, lty=2)
lines(km.pohlavi[pohlavi=2]$time, km.pohlavi[pohlavi=2]$upper, type="s", col=2, lty=2)

# logrank test pro rozdil v preziti:
km.pohlavi.test <- survdiff(Surv(OS_m, OS_cenz) ~ pohlavi, data = aml)
km.pohlavi.test