Survival Analysis | R

#Set working directory
setwd("C:/Users/kumar/Desktop/ABA/Master")

#Read in data
#install.packages("readr")
library(readr)
survdat<-read.csv("Survival.csv")
names(survdat)

#create censoring variable
survdat$censored[survdat$Turnover==1]<-1
survdat$censored[survdat$Turnover==0|survdat$Turnover==2]<-0

#Inspect Distribution of LOS Variable
hist(survdat$LOS)
hist(subset(survdat,Turnover=0)$LOS)
hist(subset(survdat,Turnover=1)$LOS)
hist(subset(survdat,Turnover=2)$LOS)

#Kaplan-Meier Analysis(KM Analysis)- LifeTable
library(survival)
km_fit1<-survfit(Surv(LOS,censored)~1,
                data=survdat,
                type="kaplan-meier")
print(km_fit1)
summary(km_fit1)

#summarize results by pre-specified time intervals
summary(km_fit1,times=c(30,60,90*(1:30)))

#install.packages(("survminer"))
library("survminer")

ggsurvplot(km_fit1,data=survdat,risk.table=TRUE,conf.int=TRUE,ggtheme=theme_minimal())

#Categorical KM Analysis:Race
km_fit2<-survfit(Surv(LOS,censored)~Race,
                 data=survdat,
                 type="kaplan-meier")
print(km_fit2)
summary(km_fit2)

#Categorical KM Analysis: Gender
km_fit3<-survfit(Surv(LOS,censored)~Gender,
                 data=survdat,
                 type="kaplan-meier")
print(km_fit3)
summary(km_fit3)

cox_reg<-coxph(Surv(LOS,censored)~1,data=survdat)#null-model
summary(cox_reg)

cox_reg1<-coxph(Surv(LOS,censored)~Race,data=dropna(survdat,LOS,censored,Race,Pay_hourly,Pay_sat))
summary(cox_reg1)
#Note1: Alphabetically, Black is the reference group.Survival Rate of survival of Hispanic Latino and  white wrt to Black
#Note2: concordance is how well did the model predict how long people are going to survive before the event.
#Note3: logrrank test compares model with covariate with model with no covariate/null model
#Obsv1: Since the p-values>0.05, the variables are insignificant. Therefore, we do not discuss the hazard wrt to them.

cox_reg2<-coxph(Surv(LOS,censored)~Race+Pay_hourly+Pay_sat,data=dropna(survdat,LOS,censored,Race,Pay_hourly,Pay_sat))
summary(cox_reg2)
#Note1: Sample size has changed to 663. It is because Pay_sat column has missing values.
#Obsv1: Hourly pay is significant. It is negative.Since it is significant,we look at the hazard ratio.
#1- hazard ratio=1-0.8961=0.1039. Interpretably,The risk of an individual voluntarily turning over/experiencing the focal
#event is going to decrease by 10.4% for every additional $ in hourly pay that the individual actually earns

#Application of the regression model
RaceHL=1#ishispaniclatino
RaceW=0
PH=16.00
PS=4.00

Log_overallrisk=.121*RaceHL-0.044*RaceW-0.109*PH+0.14135*PS
print(exp(Log_overallrisk))
#Risk of this individual turning over is 0.347.This individual is 65.27% less likely to quit when compared to individual
#with scores of 0 on each of the covariance in the model. But in reality Pay_Hourly and Pay_sat cannot be 0.

#Therefore, we mean-center the covariates
survdat$c_Pay_hourly<-scale(survdat$Pay_hourly,center=TRUE,scale=FALSE)
survdat$c_Pay_sat<-scale(survdat$Pay_sat,center=TRUE,scale=FALSE)

#Also,change the reference group from black to HispanicLatino
survdat$HL_Race<-factor(survdat$Race,levels=c("Hispaniclatino","Black","White"))

cox_reg3<-coxph(Surv(LOS,censored)~HL_Race+c_Pay_hourly+c_Pay_sat,data=dropna(survdat,LOS,censored,Race,Pay_hourly,Pay_sat))
summary(cox_reg3)

RaceB=0
RaceW=0
PH=16.00-mean(survdat$Pay_hourly,na.rm=TRUE )#na.rm will drop the missing data
PS=4.00-mean(survdat$Pay_sat,na.rm=TRUE)

Log_overallrisk=.121*RaceB-0.044*RaceW-0.109*PH+0.14135*PS
print(exp(Log_overallrisk))
##Risk of this individual turning over is 0.851.This individual is 14.8% less likely to quit when compared to individual
#with scores of 0 on dummy variables and mean Pay_hourly and Pay_Sat

library(tidyr)
anova(cox_reg,cox_reg1)#number of observations should be same