# Practical Multivariate Analysis, 5th Edition
# Chapter 13


ct1 <- crosstabs( ~ staget+death)
print(ct1, marginal.totals=F)
ct2 <- crosstabs( ~ perfbl[perfbl != "NA"]+death[perfbl != "NA"])
print(ct2, marginal.totals=F)
ct3 <- crosstabs( ~ treat+death)
print(ct3, marginal.totals=F)
ct4 <- crosstabs( ~ poinf[poinf != "NA"]+death[poinf != "NA"])
print(ct4, marginal.totals=F)

#table 13.3, p. 318.
surv.weibull <- censorReg( censor(days, death) ~ staget+perfbl+poinf+treat,
                  surv, distritution="weibull", na.action=na.omit)
summary(surv.weibull)

#table 13.4, p. 320.
surv.cox <- coxph( Surv(days, death) ~ staget+perfbl+poinf+treat,
                  surv, na.action=na.omit)
summary(surv.cox)

#Obtaining the survival functions.
fit <- survfit(surv.cox, conf.type="none")
fit$surv2 <- fit$surv^exp(.54)

#calculating the log(-log(S(t))
fit$log.surv <- log( -log(fit$surv))
fit$log.surv2 <- log( -log(fit$surv2))

#Plotting
plot( fit$time, fit$log.surv2 , type="s", xlab="time",
      ylab="log(-log(S(t)))" )
points(fit$time, fit$log.surv,  type="s")