rm(list=ls())

library(survival)

### survival function for the survival times of hepetitis patients 


#### MY WAY #####
hepa = read.table("http://www.sci.csueastbay.edu/~esuess/Statistics_65016502/Handouts/2011/6502/survivalfn/hepatitis.dat",
header=FALSE, quote="")

colnames(hepa) = c("time","status","group")
#lty=0 lty="blank"
#lty=1 lty="solid" 
#lty=2 lty="dashed"
#lty=3 lty="dotted"
#lty=4 lty="dotdash"
#lty=5 lty="longdash"
#lty=6 lty="twodash"

#### SUESS WAY ####
#hepa = matrix(scan("C:\\whale\\hepatitis.dat"),ncol=3,byrow=T)
#dimnames fills names row("NULL") by cols
#dimnames(hepa) = list(NULL, c("time","status","group"))		
#dimnames(hepa) 

hepa.data = data.frame(hepa)

hepa.data

hepa.data$time
hepa.data$status
hepa.data$group

#?survival::survfit 
#Compute a Survival Curve for Censored Data
#Computes an estimate of a survival curve for censored data using either the Kaplan-Meier 
#or the Fleming-Harrington method or computes the predicted survivor function. 
#For competing risks data it computes the cumulative incidence curve. 
#See survfit.coxph for survival curves from a fitted Cox model. 

fit = survfit(Surv(time, status) ~ 1, data=hepa.data)

X11()
plot(fit)

# plot the survival functions for the hepatitis data, Kaplan-Meier

fit = survfit(Surv(time, status) ~ group, data=hepa.data)

X11()
plot(fit, lty = 1:6, xlab="Weeks elapsed",ylab="Proportion of live patients")
summary(fit)

# Fit an exponential model

fit.exp = survreg(Surv(time, status) ~ group, data = hepa.data, dist="exponential")
summary(fit.exp)