## "lkr" dataset (Garrett, 1997) includes 5 variables for a leukemia remission study:
## "weeks", "relapse", "treatment1", "treatment2" and "wbc3cat".
## "weeks": the time that a patient was remained in remission (in weeks).
## "relapse": whether the patient was relapsed - 1 = yes, 0 = no.
## "treatment1": 1=drug A, 0 = standard drug for treatment1.
## "treatment2": 1=drug B, 0 = standard drug for treatment2.
## "wbc3cat": white blood cell count, recorded in three categories -
## 1 = normal, 2= moderate, 3 = high.
## load the "threg" package
library("threg")
## load the data "lkr"
data("lkr", package = "threg")
## Transform the "treatment2" variable into factor variable "f.treatment2" .
lkr$f.treatment2 <- factor(lkr$treatment2)
## generate the Kaplan-Meier survival curves for the drug B group and
## the standard group.
fit <- survfit(Surv(weeks, relapse) ~ f.treatment2, data = lkr)
plot(fit, lty = 1:2, xlab="time",ylab="Estimated S(t)")
legend(20, 1, c("Standard","Drug B"), lty = 1:2)
## fit the threshold regression model on the factor variable "f.treatment2",
fit <- threg(Surv(weeks, relapse) ~ f.treatment2 | f.treatment2, data = lkr)
fit
## generate the threshold regression predicted survival curves for the drug B group and
## the standard group.
plot(fit, var = f.treatment2, graph = "sv", nolegend = 1, nocolor = 1)
legend(20, 1, c("Standard", "Drug B"), lty = 1:2)
## generate the threshold regression predicted hazard curves for the drug B group and
## the standard group.
plot(fit, var = f.treatment2, graph = "hz", nolegend = 1, nocolor = 1)
legend(1, .18, c("Standard","Drug B"), lty = 1:2)
## calculate the hazard ratio of the drug B group v.s. the standard group at
## week 5 (this hazard ratio is calculated as 2.08)
hr(fit, var = f.treatment2, timevalue = 5)
## calculate the hazard ratio of the drug B group v.s. the standard group at
## week 20 (this hazard ratio is calculated as 0.12)
hr(fit, var = f.treatment2, timevalue = 20)
## As a comparison, fit the Cox proportion hazards model on "f.treatment2",
## and the Cox model gives a constant hazard ratio, 0.73.
library("survival")
summary(coxph(Surv(weeks, relapse) ~ f.treatment2, data = lkr))
## generate the Cox model predicted survival curves for the drug B group and
## the standard group.
fit <- coxph(Surv(weeks, relapse) ~ f.treatment2, data = lkr)
plot(survfit(fit, newdata=lkr[lkr$f.treatment2==0,]
),lty = 1,conf.int=F,xlab="time",ylab="Estimated S(t)")
lines( survfit(fit, newdata=lkr[lkr$f.treatment2==1,]
), lty = 2)
legend(15, 1, c("Standard","Drug B"), lty = 1:2)
## the bone marrow transplantation dataset "bmt"
## (Klein and Moeschberger, 2003) includes 12 variables, but we will only use
## the following 5 variables in our example:
## "time", "indicator", "recipient_age", "group" and "fab".
## "time": time to relapse, death or end of study (in days).
## "indicator": censoring indicator variable - 1 = dead or replapse, 0 = otherwise.
## "recipient_age": patient age.
## "group": risk categories based on their status at the time of transplantation.
## "fab": French-American-British (FAB) classification based on standard
## morphological criteria.
## load the data "bmt"
data("bmt", package = "threg")
## Transform the "group" and "fab" variables into factor variables
## "f.group" and "f.fab".
bmt$f.group <- factor(bmt$group)
bmt$f.fab <- factor(bmt$fab)
## fit a threshold regression model on the "bmt" dataset, by using "recipient_age" and
## "f.fab" as the predictors for ln(y0), and "f.group" and "f.fab" as predictors for mu.
fit <- threg(Surv(time, indicator) ~ recipient_age + f.fab | f.group + f.fab,
data = bmt)
fit
## Calculate the hazard ratio for
## "f.group" for the specified scenario that "the patient age is 18 years old and
## the FAB classification is 0", at the time ``500 days''.
hr(fit, var = f.group, timevalue = 500, scenario = recipient_age(18) + f.fab1(0))
## fit the same model as above, but additionally overlay curves of survival functions
## corresponding to different levels of "f.group'.
plot(fit, var = f.group, scenario = recipient_age(18) + f.fab1(0), graph = "sv", nocolor = 1)
## fit the same model as above, but additionally overlay curves of hazard functions
## corresponding to different levels of "f.group'.
plot(fit, var = f.group, scenario = recipient_age(18) + f.fab1(0), graph = "hz", nocolor = 1)
## fit the same model as above, but additionally overlay curves of probability density
## functions corresponding to different levels of "f.group'.
plot(fit, var = f.group, scenario = recipient_age(18) + f.fab1(0), graph = "ds", nocolor = 1)
## predict lny0, y0, mu, f, S and h for the specified scenario that "the patient age is
## 18 years old, the FAB classification is 0 and the risk category is 3", at the
## time ``2000 days''
predict(fit, timevalue = 2000, scenario = recipient_age(18) + f.fab1(0) + f.group2(0) + f.group3(1))