# PBC07.S
# xref: PBC.DOC
# input:   PBC (S object)
# output:
# does: - exlporing the PBC data set
#       - illustrate survival plotting in S
#       - illustrate fitting a P.H. model
#******************************************************
cut4 <- function (x)    # cut x into 4 categories of equal size
{
  cut (x, breaks=c(-1, 0, 0, 0, 1) + summary (x) [c(1:3,5:6)])
}
#******************************************************
# An overall survival curve
all  <- surv.fit(pbc$time, pbc$delta)
plot.surv.fit (all)
title (main="Estimated S(t) - overall")
scan()
#******************************************************
# consider the covariate spider naevi: 0=no, 1=yes.

# frequencies on z6 (spider-naevi)
table(pbc$z6)

# estimate survival for each z6 group
by.z6 <- surv.fit (pbc$time, pbc$delta, pbc$z6)

# plot survival by group of z6
plot.surv.fit ( by.z6, lty=1:2 )
title (main="Estimated S(t) by  spider-naevi")
scan()

# plot survival by group of z6, with C.I.
plot.surv.fit ( by.z6, conf.int=T, lty=1:2 )
title (main="Estimated S(t) by  spider-naevi")
scan()

#plot on the log scale i.e H(.)
plot.surv.fit ( by.z6, log=T, lty=1:2 )
title (main="Estimated S(t) by  spider-naevi")
scan()

#******************************************************
# consider the covariate prothrombin time (z16)
# plot survival by z16 prothrombin time
by.z16 <- surv.fit (pbc$time, pbc$delta, cut4(pbc$z16))
plot.surv.fit ( by.z16, log=T, lty=1:4 )
title (main="Estimated S(t) by prothrombin time")
scan()
#******************************************************
# fit a p.h. model
coxreg(pbc$time, pbc$delta, pbc$z6)

# fit a p.h. model, save results
pbc.ph  <- coxreg(pbc$time, pbc$delta, cbind(pbc$z6, pbc$z16))
print(pbc.ph)