########## SURVIVAL ANALYSIS ###################

#### 1. Brain Cancer Data ####

install.packages("ISLR2")
library(ISLR2)
fix(BrainCancer)
attach(BrainCancer)
table(sex)           #frekvencije za kategorijska obelezja sex, diagnosis, status
table(diagnosis)
table(status)     # 1 znaci da je pacijen umro pre kraja studije

#Crtanje Kaplan-Majerove krive

library(survival)  #ili stikliraj sa desne strane medju paketima
KMcurve=survfit(Surv(time,status)~1) #nakon znaka ~ uvek ide oznaka za ishod SMRT
plot(KMcurve,xlab="Months",ylab="Estimated Probability of Survival")

#Razdvajanje krive po polu

KMcurve1=survfit(Surv(time,status)~sex)
plot(KMcurve1,xlab="Months",ylab="Estimated Probability of Survival", col=c(2,4))
legend("bottomleft", levels(sex), col=c(2,4),lty=1)  #lty je line type (probaj da menjas od 1-5)

#Log-rank test

logrank=survdiff(Surv(time,status)~sex) 
logrank

#Cox-ova regresija

Cox=coxph(Surv(time,status)~sex) #samo sex kao prediktor
summary(Cox) #p=0,05 odgovara z vrednosti od 1,96

#Komentar: Skor test jednak Log-rank testu (uporedi prakticno). Dalje, sve p-vrednosti su zaokruzene
#na jednu decimalu. Pune vrednosti se mogu videti

summary(Cox)$sctest
summary(Cox)$logtest
summary(Cox)$waldtest

Cox1=coxph(Surv(time,status)~sex+diagnosis+loc+ki+gtv+stereo) #ukljuceni svi prediktori
Cox1 # ovde istakni razliku izmedju Cox1 i summary(Cox1) (intervali poverenja i dodatni testovi)

#Komentar (kontrolne grupe). diagnosis~meningioma, sex~female, loc~Infrantentorial 


#### 2. Publication #####

fix(Publication)
attach(Publication)
table(status) # oznaka 1 znaci da je rad objavljen podaci tokom ove studije
table(posres)
Pub=survfit(Surv(time, status)~posres, data=Publication) #ovde je potrebno navesti koje podatke koristimo, jer se obelezja Time i Status javljaju u dve baze
plot(Pub , xlab = "Months",
     ylab = "Probability of Not Being Published", col = c(3,2))
legend("topright", c("Negative Result", "Positive Result"),
       col = c(3,2), lty = 1)
#Komentar. Za razliku od prethodnog slucaja, ovde ne moze stajati komanda 'levels(posres)', 
          #jer su kodovi za posres numericke vrednosti 0 i1. Zato unosimo stringove koji se se ispisati
          #u legendi.

#Cox-ova regresija

CR=coxph(Surv(time,status)~posres, data=Publication)
CR # vrednost test statistike Likelihood ratio testa je ista kao i kod Log-rank testa.

#Log-rank

LK=survdiff(Surv(time,status)~posres, data=Publication)
LK
LK$chisq #ako test u osnovi pokaze vrednost samo 0.8, pa je potrebno uociti stvarno poklapanje

CR1=coxph(Surv(time,status)~.-mech, data=Publication) #obelezje mech je izbaceno jer ima veliki broj kategorija
CR1

#### 3. Call Center ######

#Ovde se simuliraju podaci, nema gotove baze

set.seed(4) #obezbedjuje da random generator svaki put generise iste vrednosti
N=2000
Operators=sample(c(5:15), N, replace = T) #funkcija 'sample' pravi uzorak duzine N od vektora unetih vrednosti
Center=sample(c("A", "B", "C"), N, replace = T)
Time=sample(c("Morn.", "After.", "Even."), N, replace = T)
X=model.matrix( ~ Operators + Center + Time) #kreira matricu od unetih kolona
X[1:15, ]
X=model.matrix( ~ Operators + Center + Time)[,-1]# da bismo neutralisali prvu kolonnu 'Intercept'

true.beta=c(0.04,-0.3,0,0.2,-0.2)
hazard.function=function(x) return(0.00001*x)

# U gornjem delu su generisane vrednosti prediktora u pocetnom trenutku. 
# Zatim je dat Beta vektor koji treba da generise vreme prezivljavanja. 
#Hazard funkcija je data u istu svrhu i predstavlja h_0.

install.packages("coxed")
library(coxed)
Prezivljavanje=sim.survdata(N = N, T = 1000, X = X,     # T=1000 znaci da posmatramo ispitanike u 1000 jedinica vremena
        beta = true.beta, hazard.fun = hazard.function)
names(Prezivljavanje)
Prezivljavanje$data[1:50,] #dodatno je generisano vreme y i status
mean(Prezivljavanje$data$failed) #True znaci da je doslo do kriticnog dogadjaja

# Kaplan-Majer krive razdvojene po obelezju 'Center'

par(mfrow = c(1, 2))

KM.Center=survfit(Surv(y, failed) ~ Center, data = Prezivljavanje$data)
plot(KM.Center , xlab = "Seconds",
       ylab = "Probability of Still Being on Hold", col = c(1,2,3))
legend("bottomleft", c("Call Center A", "Call Center B", "Call Center C"),
         col = c(1,2,3), lty = 1)

KM.Time <- survfit(Surv(y, failed) ~ Time, data = Prezivljavanje$data)
plot(KM.Time, xlab = "Seconds", ylab = "Probability of Still Being on Hold",
       col = c(4,5,6))
legend("topright", c("Morning", "Afternoon", "Evening"),
         col = c(4,5,6), lty = 1)

#Za povratak na punu sluku grafika, prvo moramo uneti
par(mfrow = c(1, 1))

#Log-rank test
LR=survdiff(Surv(y, failed) ~ Center , data = Prezivljavanje$data)
LR
LR1=survdiff(Surv(y, failed) ~ Time , data = Prezivljavanje$data)
LR1

# Cox-ova regresija (cisto da uporedis rezultate sa pocetnim vektorom 'true.beta')

CReg=coxph(Surv(y,failed)~.,data=Prezivljavanje$data)
CReg