Survival model

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

plot(cars)

##  The hmohiv data-set is drawn from a study of HIV positive patients. 
#The study examined whether there was a difference in survival times of HIV 
#positive patients between those who had used intravenous drugs and those 
#who had not. This data-set has been taken from 

## Quesion: Is there a difference in survival times of HIV positive patients between 
#those who had used intravenous drugs and those who had not?


## Read data
hmohiv<-read.table("https://stats.idre.ucla.edu/stat/r/examples/asa/hmohiv.csv", sep=",", header = TRUE)
attach(hmohiv)
hmohiv

###  Read and view data
##  To check to see what object (format) that you have created
class(hmohiv)

##  To confirm the data that is contained within hmohiv
hmohiv 

##  In order to access variables 
hmohiv$ID 
hmohiv$time 

##  To View patient ID and age    
hmohiv[,c(1,3)]

###  Calculate survival time  
##  1) Create a new column
hmohiv["TIME_CHECK"]<-NA 

##  2) Calculate the difference, in days between end date and start date 
hmohiv$TIME_CHECK<-as.numeric(as.Date(hmohiv[,7], "%m/%d/%Y")-as.Date(hmohiv[,6], "%m/%d/%Y")) 

##  3) Check on output
hmohiv$TIME_CHECK 

##  4) Aggregate from days to months  
hmohiv$TIME_CHECK<-round(hmohiv$TIME_CHECK/30.5) 

##  5) Double-check your answer against the original time in the table
hmohiv[,c(2,8)] 

###  Survival anaysis
##  Load the survival package: 
library(survival) 

###  Create a Survival Object 
##  The initial step is to create a survival object
s_obj<-Surv(time, censor) 

##  Look at it in the context of time/censoring. Create a new column and copy the results of survival object into it 
hmohiv["s_obj"]<-NA 
hmohiv$s_obj<-Surv(time, censor) 

##  Compare it to the time and censor columns 
hmohiv[,c(2,5,9)] 

##  In the survival object, patients with censored times are represented by 4+, 1+ etc. 
##  Patients who experienced an event have unaltered times, eg 7, 9 etc. 
with(hmohiv, Surv(time, censor)) 

## Descriptive Statistics
survfit(Surv(time, censor)~ 1) 
summary(survfit(Surv(time, censor)~1))

# 100 people in the HIV dataset
# 80 people are uncensored (followed for the entire time, until occurrence of event)
# Among those 80 people, there was a median survival time of 7 months
# The 95% confidence interval for the median survival time for the 80 uncensored individuals is (5, 10)

# n.risk: the number of people at risk
# n.event: the number of events at that time point
# survival: the proportion of individuals who survived up until that point

## Nonparametric plots
surv.aml <- survfit(Surv(time,censor)~1)
plot(surv.aml, main = "Plot of Survival Curve for HIV Patients",
     xlab = "Length of Survival", ylab = "Proportion of Individuals who have Survived")

# Survival curve is also known a Kaplan-Meier plot
# The proportion of individuals who have survived up until that particular time as a solid black line
# 95% confidence interval as the dashed lines

## Going back to original questions
print(survfit(Surv(time,censor)~drug))

# Is there a difference in survival times of HIV positive patients between 
# those who had used intravenous drugs and those who had not?
# The description of data is broken into drug (=0) and drug (=1)

# 49 people on drug 
# 38 people died

plot(survfit(Surv(time,censor)~drug), main = "Plot of Survival Curves by Drug", xlab = "Length of Survival",ylab="Proportion of Individuals who have Survived",col=c("blue","red"))

legend("topright", legend=c("Drug=0", "Drug=1"),fill=c("blue","red"),bty="n")

# Is there a difference in survival times of HIV positive patients between those 
# who had used intravenous drugs and those who had not?
# Log rank statistics
# where ∑▒𝑂_𝑗𝑡  represents the sum of the observed number of events in the jth group over t
# ime and ∑▒𝐸_𝑗𝑡  represents the sum of the expected number of events in the jth group over time

## Statistic 
survdiff(Surv(time,censor)~drug)

# We see that when testing the null hypothesis that there is a significant difference in the survival 
# function for those who were on drug versus those who were not on drug. We reject the null hypothesis 
# Chisq = 11.9 with a p-value <.05.

## Cox proportional hazard model
# Univariate
coxph(Surv(time,censor)~factor(drug))

# This is therefore a linear model for the log-hazard or a multiplicative model for the hazard itself


# We see that when using the Cox regression to perform the test, the results are very similar to the 
# log rank test (Chisq = 11.9 with a p-value = 0.000575)
# The Cox regression estimates the hazard ratio of dying when comparing “Drug” with “No drugs”
# People have significantly better survival without drug (p-value = 0.000659)
# “No drug” is a reference level. Using drug has 2.295 times the hazard of dying in comparison to “Drug”


# Multivariate
coxph(Surv(time,censor)~factor(drug) + age)

# We consider additional risk factors. Here we consider “age”.


# After adjusting for age, peoples still have significantly better survival without drug (p-value < .05)
# The p-value for age <.05, with a hazard ratio HR = exp(coef)=1.102, indicating a strong relationship 
# between age and increased risk of death. Holding the other covariates constant, a higher value of age 
# is associated with a poor survival

detach(hmohiv)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.