### survival function for the lifetimes of guinea pigs - Example A p.381
rm(list=ls())
library(survival)

# c1 lifetime c2 status c3 group
### SUESS's WAY ###
#gpig = matrix(scan("C:\\whale\\gpig1.dat"),ncol=3,byrow=T)
#dimnames(gpig) = list(NULL, c("time","status","group"))		# dimnames fills names row by cols
#dimnames(gpig) 

### MY WAY ###
gpig = read.table("http://www.sci.csueastbay.edu/~esuess/Statistics_65016502/Handouts/2013/6501/survivalfn/GPIG1.DAT",
header=FALSE, quote="")
colnames(gpig)=c("time","status","group")

gpig.data = data.frame(gpig)

#gpig.data

#gpig.data$time
#gpig.data$status
#gpig.data$group

fit = survfit(Surv(time, status) ~ 1, data=gpig.data)
X11()
plot(fit)

# plot the survival functions for the groups on the same plot, Kaplan-Meier

fit = survfit(Surv(time, status) ~ group, data=gpig.data)
X11()
plot(fit, lty=1:6, col=2:7, lwd=2 , 
    main="Survival Functions for Guinea Pig Lifetimes 
          with Doses of Tubercle Bacilli",
    xlab="Days elapsed",ylab="Proportion of live animals")
legend(500,0.9, c("Control","I","II","III","IV","V"),lty=1:6, col=2:7, lwd=2)
summary(fit)

# comparing Kaplan-Meier Survival Curves

survdiff(Surv(time, status) ~ group, data = gpig.data)

# Fit an exponential model

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