### Example: Stratified Medicine Using Partitioned Cox-Models
### A combination of and
### based on infrastructure in the mlt R add-on package described in
### https://cran.r-project.org/web/packages/mlt.docreg/vignettes/mlt.pdf
library("trtf")
library("survival")
### German Breast Cancer Study Group 2 data set
data("GBSG2", package = "TH.data")
GBSG2$y <- with(GBSG2, Surv(time, cens))
### set-up Cox model with overall treatment effect in hormonal therapy
cmod <- Coxph(y ~ horTh, data = GBSG2, support = c(100, 2000), order = 5)
### overall log-hazard ratio
coef(cmod)
### roughly the same as
coef(coxph(y ~ horTh, data = GBSG2))
### partition the model, ie both the baseline hazard function AND the
### treatment effect
(part_cmod <- trafotree(cmod, formula = y ~ horTh | age + menostat + tsize +
tgrade + pnodes + progrec + estrec, data = GBSG2))
### compare the log-likelihoods
logLik(cmod)
logLik(part_cmod)
### stronger effects in nodes 2 and 4 and no effect in node 5
coef(part_cmod)[, "horThyes"]
### plot the conditional survivor functions; blue is untreated
### and green is hormonal therapy
nd <- data.frame(horTh = sort(unique(GBSG2$horTh)))
plot(part_cmod, newdata = nd,
tp_args = list(type = "survivor", col = c("cadetblue3", "chartreuse4")))
### same model, but with explicit intercept term and max-type statistic
### for _variable_ selection
(part_cmod_max <- trafotree(cmod, formula = y ~ horTh | age + menostat + tsize +
tgrade + pnodes + progrec + estrec, data = GBSG2, intercept = "shift",
control = ctree_control(teststat = "max")))
logLik(part_cmod_max)
coef(part_cmod_max)[, "horThyes"]
### the trees (and log-likelihoods are the same) but the
### p-values are sometimes much smaller in the latter tree
cbind(format.pval(info_node(node_party(part_cmod))$criterion["p.value",]),
format.pval(info_node(node_party(part_cmod_max))$criterion["p.value",]))
Run the code above in your browser using DataLab