## "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))