library(survival) download_data = TRUE ## Read in the data: ## if( download_data ){ DF = read.table('http://lib.stat.cmu.edu/datasets/pbc', header=FALSE, skip=60, nrows=418, na.strings=c('NA', '.')) colnames(DF) = c('ID', 'FU_Days', 'Status', 'Drug', 'Age', 'Sex', 'Ascites', 'Hepatomegaly', 'Spiders', 'Edema', 'Bile', 'Chol', 'Albumim', 'Copper', 'Alk_Phos', 'SGOT', 'TRIG', 'Platelet', 'Protime', 'Stage') DF$ID = NULL } ## Focus on the trial patients: ## DF = DF[1:312, ] liver_tx = (DF$Status==0) | (DF$Status==1) ## samples censored due to liver tx death = (DF$Status==0) | (DF$Status==2) ## samples censored by death ## Study the samples censored due to liver tx: ## if( FALSE ){ df_liver = DF[liver_tx, ] res = coxph(Surv(FU_Days, Status) ~ Drug + Protime + Alk_Phos, data=df_liver) print(res) } ## Study the samples censored due to death: ## df_death = DF[death, ] df_death$Status[df_death$Status==2] = 1 ## Split this group into treated and placebo groups: ## df_treated = df_death[df_death$Drug==1, ] ## df_placebo = df_death[df_death$Drug==2, ] ## placebo res_treated = coxph(Surv(FU_Days, Status) ~ Protime + Alk_Phos, data=df_treated) print(res_treated) res_placebo = coxph(Surv(FU_Days, Status) ~ Protime + Alk_Phos, data=df_placebo) print(res_placebo) ## Log-rank survival difference: ## res = survdiff(Surv(FU_Days, Status) ~ Drug, data=df_death) print(res)