This package implements basic methodology for evaluating one or more biomarkers for their ability to guide patient treatment recommendations.
The methodology assumes that the data come from a randomized and controlled trial (RCT) comparing two treatment options, which we refer to as "treatment" and "no treatment". These could be, for example, two different active prophylactic or therapeutic interventions, or an experimental intervention and a standard of care. Subjects are followed for the development of a binary clinical outcome or event within a specified time-frame following treatment/no treatment. The biomarker(s) are assumed to be measured at baseline; they could be anything from patient demographics or clinical characteristics, to traditional biomarker measurments or the results of imaging or genetic or proteomic analyses. The methodology accommodates settings where the biomarker(s) are measured on all RCT participants, and also settings in which participants are retrospectively sub-sampled based on clinical outcome and/or treatment assignment for marker measurement.
This tutorial uses the Treatment Selection package (v 2.1.0) to analyze the example data provided in the package.
The current version of the package is on github. Use the devtools
package to install the package:
if (!require("devtools")) install.packages("devtools") devtools::install_github("mdbrown/TreatmentSelection")
First, load the example data set for this package called tsdata
. Four markers are included in the data example, a ''weak'' and a ''strong'' marker ($Y1$ and $Y2$ respectively), along with a weak/strong discrete markers.
Load the example data set for this package called tsdata
. These are hypothetical data from an RCT with a binary outcome (outcome
) and treatment assignment (trt
). Four markers, measured on all RCT participants, are included-- a ''weak'' and a ''strong'' marker (Y1
and Y2
respectively), along with binary versions of these markers.
library(TreatmentSelection) set.seed(12321) data(tsdata) data(surv_tsdata) tsdata[1:5, ]
For a pre-defined marker-based treatment rule, we provide a simple function to estimate point estimates of performance measures. This is done using the trtsel_measures()
function. The user must specify a vector of clinical outcomes, a vector of treatment assigments, and a vector of marker-based treatment recommendations based on the pre-specified rule.
Here we let Y1_disc
represent a user-specified treatment rule and evaluate its performance.
trtsel_measures(event = tsdata$event, trt = tsdata$trt, trt.rule = tsdata$Y1_disc, default.trt = "trt all" ) trtsel_measures(event = surv_tsdata$di, time = surv_tsdata$xi, trt = surv_tsdata$trt, trt.rule = as.numeric(surv_tsdata$Y > 0), prediction.time = 1, default.trt = "trt none" )
We can also fit our own risk model using GLM, use this model to develop a marker-based treatment recommendation, and evaluate its performance. This allows us to obtain model-based estimates of performance:
mod <- glm(event~trt*Y1_disc, data = tsdata, family = binomial()) tsdata.0 <- tsdata; tsdata.0$trt = 0 tsdata.1 <- tsdata; tsdata.1$trt = 1 delta.hat <- predict(mod, newdata= tsdata.0, type = "response") - predict(mod, newdata= tsdata.1, type = "response") trtsel_measures(event = tsdata$event, trt = tsdata$trt, trt.rule = 1- tsdata$Y1_disc, trt.effect = delta.hat )
An alternative to providing a pre-specified treatment rule the user can use the data to fit a risk model and to develop a treatment rule based on estimated treatment effect.
The user can do this by first creating a treatment selection R object using the function trtsel
. Creating this object entails estimating risk of the clinical outcome given treatment and marker, i.e. fitting a "risk model". Here, logistic regression is used (family = binomial(link="logit")
) and the formula is specified in the first argument. The default treatment strategy-- that which would be recommended in the absence of marker measurements-- is specified as default.trt="trt all"
. The methodology considers marker-specific treatment rules or policies that recommend treatment if the marker-specific treatment effect-- on the risk difference scale-- is above a specified thresh
; the default is thresh = 0
. (Alternatively, see the tab on "additional features"" to see how an arbitrary user-specified treatment rule can be evaluated.) The study design="RCT"
indicates that the markers are measured on all RCT participants.
trtsel.Y1 <- trtsel(event ~ Y1*trt, treatment.name = "trt", data = tsdata, family = binomial(), study.design = "RCT", default.trt = "trt all") trtsel.Y1
As we see above, the trtsel.Y1
object contains information about the study design, the fitted risk model, the fitted risk given each treatment, and the estimated marker-specific treatment effect (on the risk difference scale) for each individual. This object also contains the treatment recommendation for each subject, based on a treatment rule that recommends treatment if the estimated treatment effect is above thresh
; here thresh = 0
.
This is what a trtsel
object looks like for a binary marker.
# Y2_disc = as.numeric(Y2>0) trtsel.Y2_disc <- trtsel(event ~ trt*(Y2_disc), treatment.name = "trt", data = tsdata, study.design = "RCT", family = binomial("logit"))
See ?trtsel
for more information.
Now that we have created trtsel objects, we can assess the calibration of the risk models; and plot, evaluate, and compare markers.
It is important to assess the fit of the risk model to the data. The package implements methods for examining calibration.
Here we compare (average) predicted and observed risk values by decile ('group = 10') of fitted risk, for each treatment group. The risks are shown on the log scale. The Hosmer-Lemeshow goodness of fit test is used to detect evidence of mis-calibration by treatment group.
calibrate(trtsel.Y1, groups = 10, plot.type = "treatment effect", trt.names = c("chemo.", "no chemo."), point.color = "lightcoral", line.color = "black") calibrate(trtsel.Y1, groups = 10, plot.type = "risk.t0", trt.names = c("chemo.", "no chemo."), point.color = "lightcoral", line.color = "black")
See ?calibrate.trtsel
for more plot options.
The plot()
function allows us to visually assess the ability of the marker or markers to predict the treatment effect.
Risk curves apply to individual markers. They are plots of estimated risk of the clinical outcome as a function of the marker, for each treatment group. The two risk curves are aligned with respect to marker percentile, but the option show.marker.axis=TRUE
adds a second x-axis showing the marker values corresponding to each percentile. Pointwise confidence intervals, obtained using the bootstrap, can be added (ci="vertical"
).
plot(trtsel.Y1, main = "Y1: Oncotype-DX-like marker", plot.type = "risk", ci = "vertical", conf.bands = TRUE, bootstraps = 50, #more bootstraps should be run than this in practice! trt.names = c("chemo.", "no chemo."), show.marker.axis = FALSE)
This is what the plot looks like for a binary marker:
tmp <- plot(trtsel.Y2_disc, main = "Discrete version of Y2", plot.type = "risk", ci = "vertical", conf.bands = TRUE, bootstraps = 50, trt.names = c("chemo.", "no chemo."))
Note that tmp
is a list with elements plot
that holds the ggplot output, and ci.bounds
which holds the information regarding the confidence bounds.
tmp$ci.bounds
Treatment effect curves apply both to univariate and multivariate markers. They show the distribution of marker-specific treatment effects. Using plot.type = "treatment effect"
produces a reverse-CDF-type plot:
plot(trtsel.Y1, plot.type = "treatment effect", ci = "horizontal", conf.bands = TRUE, bootstraps = 50)
Alternatively, using plot.type = "cdf"
produces a traditional CDF of the treatment effects:
plot(trtsel.Y1, plot.type = "cdf", ci = "vertical", conf.bands = TRUE, bootstraps = 50)
Here is an example of a treatment effect plot for a discrete marker:
plot(trtsel.Y2_disc, plot.type = "treatment effect", conf.bands = TRUE, bootstraps = 50)
The selection impact plot also applies to both univariate and multivariate markers. They show the population rate of the clinical outcome as a function of the proportion of the population recommended treatment according a marker-specific treatment rule. Rules that recommend treatment based on the marker-specific treatment effect (on the risk difference scale) are considered, and are indexed by the proportion of the population recommended treatment according to the rule. The outcome rates under "treat all" and "treat none" policies are shown for comparison.
plot(trtsel.Y1, plot.type = "selection impact", ci = "vertical", conf.bands = TRUE, bootstraps = 50)
The function evaluate()
provides estimates of measures of marker performance. A wide variety of performance measures are estimated-- some which depend on a specified treatment rule and some which do not. For the former, the rule is inherited from the trtsel
object. Confidence intervals are calculated using the quantile bootstrap; in every bootstrap sample the risk model is re-fit and the model performance measures are re-estimated.
tmp <- evaluate(trtsel.Y1, bias.correct = TRUE, bootstraps = 50) tmp
The evaluate
function used this way both fit and evaluates a risk model using the same data, which is known to bias performance measure estimates to be overly optimistic. To correct for this bias, the evaluate function by default sets bias.correct=TRUE
. This implements bootstrap bias correction via the ``refined bootstrap" method described in chapter 17 of An introduction to the Bootstrap by Efron and Tibshirani 1994.
In short, we sample $B$ bootstrap datasets. For each, obtain a new treatment selection rule based on the re-fit model and calculate the difference in estimated performance of this rule using the bootstrap vs. original data. The average of these differences estimates the bias. We shift naive performance estimates and confidence intervals down by the estimated bias.
Here is what performance evaluation looks like for a discrete marker:
# discrete marker evaluate(trtsel.Y2_disc, bootstraps = 50)
Two markers-- or marker-specific treatment rules-- can be compared using the compare()
function. Note that the compare()
function requires that that the treatment and outcome labels are identical for the two markers/treatment rules.
Here, we compare the continuous markers Y1
and Y2
. The first step is to create a trtsel
object for Y2
.
# trtsel object for markers 1 and 2 trtsel.Y2 <- trtsel(event~trt*(Y2 + Y1), treatment.name = "trt", data = tsdata, default.trt = "trt all")
A suite of performance measures is estimated for each marker, and marker-comparisons are based on Wald hypothesis tests for differences in measures with variances estimated using the bootstrap.
# Compare the markers based on summary measures mycompare <- compare(x = trtsel.Y1, x2 = trtsel.Y2, bias.correct =TRUE, model.names = c("Y1", "Y1 + Y2"), bootstraps = 50, ci = "vertical") mycompare
## Compare two discrete markers Y1_disc = as.numeric(Y1>mean(Y1)) trtsel.Y1_disc <- trtsel(event~Y1_disc*trt, treatment.name = "trt", data = tsdata) compare(x = trtsel.Y1_disc, x2 = trtsel.Y2_disc, bias.correct = TRUE, ci = "vertical", offset = 0.2, bootstraps = 50, conf.bands = TRUE, annotate.plot = TRUE)
See ?compare.trtsel
for more options.
Version 2.1.0 of the package allows for continous and time-to-event outcomes.
When using a continuous outcome, it is assumed that higher values are associated with worse outcomes.
tsdata$event.c <- tsdata$event + rnorm( nrow(tsdata), mean = 0, sd = .1 ) trtsel.Yc <- trtsel(event.c ~ Y1*trt, treatment.name = "trt", data = tsdata, study.design = "RCT", family = gaussian("identity"), #specify model error dist. default.trt = "trt all") #plot, calibrate, evaluate and compare all work as usual
When evaluating a time-to-event outcome, a landmark prediction time must be specified. Cox proportional hazards models are used for all model-based estimates. Absolute risk estimates are calculated by pairing the Cox model with the Nelson-Aalen baseline hazard estimate.
library(survival) data("surv_tsdata") #simulated sample data set with time-to-event outcome trtsel.Ysurv <- trtsel(Surv(xi, di) ~ Y*trt, treatment.name = "trt", data = surv_tsdata, prediction.time = 4, #landmark prediction time must be set study.design = "RCT", default.trt = "trt all") #notice 'censoring weights' can now be found in trtsel.Ysurv$'derived.data' #plot, calibrate, evaluate and compare all work as usual
An alternative to specifing one or more markers and fitting a risk model is for the user to specify fitted risks given treatment and no treatment for each subject. Importantly, in this case the inference provided based on bootstrap confidence intervals is conditional on these user-specified fitted risks.
Here, we illustrate this by fitting a risk model for Y2
using the glm()
function, and then feed the fitted risks into the trtsel()
function:
#calculate model fit mymod <- glm(event~trt*(Y2), data= tsdata, family = binomial("logit")) tsdata$fitted.t0 <- predict(mymod, newdata=data.frame(trt = 0,Y1 = tsdata$Y1, Y2 = tsdata$Y2), type = "response") tsdata$fitted.t1 <- predict(mymod, newdata=data.frame(trt = 1,Y1 = tsdata$Y1, Y2 = tsdata$Y2), type = "response") myfitted.trtsel <- trtsel( event~trt, treatment.name = "trt", data = tsdata, fittedrisk.t0 = "fitted.t0", fittedrisk.t1 = "fitted.t1", study.design = "RCT", default.trt = "trt all")
We can now use this trtsel
object to plot and evaluate performance just as before, but confidence intervals are narrower because they are conditional on the fitted risk model.
plot(myfitted.trtsel, bootstraps = 50, plot.type = "risk", ci = "vertical", show.marker.axis = FALSE)
Efron B, Tibshirani R. An Introduction to the Bootstrap. New York: Chapman& Hall, 1994.
Song X, Pepe MS. Evaluating markers for selecting a patient's treatment. Biometrics. 2004;60(4):874-883.
Huang Y, Sullivan Pepe M, Feng Z. Evaluating the predictiveness of a continuous marker. Biometrics. 2007;63:1181-8.
Janes H, Pepe MS, Bossuyt PM, et al. Measuring the performance of markers for guiding treatment decisions. Ann Intern Med. 2011;154(4):253-259.
Janes H, Brown MD, Huang Y, et al. An approach to evaluating and comparing biomarkers for patient treatment selection. Int J Biostat. 2014;10(1):99- 121.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.