source("http://192.38.117.59/~linearpredictors/datafiles/readPbc3.R")
library(survival)

## percent of each gender in the treatment groups
prop.table(xtabs(~ tment + sex, data = pbc3), margin = 1)

## mean and SD for variables within treatment groups
with(pbc3, 
     apply(cbind(alb, log(bili), log(alkph), log(asptr), age),
           2,
           function(x)tapply(x, tment,
                             function(i)c(mean(i, na.rm = T),
                                          sqrt(var(i, na.rm = T))))
           )
     )

## Cox proportional hazards models
## for each covariate we fit a model with this alone (gives "Covariate Effect"),
## as well as one with the covariate and treatment (gives "Treatment Effect").

coxph(Surv(followup, status != 0) ~ alb, data = pbc3)
coxph(Surv(followup, status != 0) ~ alb + tment, data = pbc3)

coxph(Surv(followup, status != 0) ~ log(bili), data = pbc3)
coxph(Surv(followup, status != 0) ~ log(bili) + tment, data = pbc3)

coxph(Surv(followup, status != 0) ~ log(alkph), data = pbc3)
coxph(Surv(followup, status != 0) ~ log(alkph) + tment, data = pbc3)

coxph(Surv(followup, status != 0) ~ asptr, data = pbc3)
coxph(Surv(followup, status != 0) ~ log(asptr) + tment, data = pbc3)

coxph(Surv(followup, status != 0) ~ age, data = pbc3)
coxph(Surv(followup, status != 0) ~ age + tment, data = pbc3)

coxph(Surv(followup, status != 0) ~ sex, data = pbc3)
coxph(Surv(followup, status != 0) ~ sex + tment, data = pbc3)