In mdbrown/TreatmentSelection: Evaluate Treatment Selection Biomarkers

{.tabset .tabset-fade .tabset-pills}

Overview of the package

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.

Install the package

This tutorial uses the Treatment Selection package (v 2.1.1) to analyze the example data provided in the package.

The current version of the package is on CRAN.

#download from CRAN
install.packages("TreatmentSelection")

Alternatively, use the devtools package to install the developmental version of the package directly from github:

#download from github (using the package 'devtools')
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.

First steps

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, ]

Evaluate performance of user-specified marker-based treatment rule

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 )

Create TrtSel objects

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.

Assess model calibration

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.

Plot

The plot() function allows us to visually assess the ability of the marker or markers to predict the treatment effect.

Risk curves

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

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)

Selection impact plot

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)

Evaluate marker performance

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)

Compare markers

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.

Additional features

Other outcome types

Version 2.1.0 of the package allows for continous and time-to-event outcomes.

Continuous

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 

Time-to-event

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 

Specifying fitted risks

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)

References

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.

mdbrown/TreatmentSelection documentation built on May 22, 2019, 3:23 p.m.