Decision Curve Analysis

Univariate Decision Curve Analysis

We’ll be working with the example dataset df_binary included with the package. The dataset includes information on 750 patients who have recently discovered they have a gene mutation that puts them at a higher risk for harboring cancer. Each patient has been biopsied and we know their cancer status. It is known that older patients with a family history of cancer have a higher probability of harboring cancer. A clinical chemist has recently discovered a marker that she believes can distinguish between patients with and without cancer. We wish to assess whether or not the new marker does indeed identify patients with and without cancer. If the marker does indeed predict well, many patients will not need to undergo a painful biopsy.

First, we want to confirm family history of cancer is indeed associated with the biopsy result.

mod <- glm(cancer ~ famhistory, df_binary, family = binomial)
tbl <- tbl_regression(mod, exponentiate = TRUE)
tbl

Via logistic regression with cancer as the outcome, we can see that family history is related to biopsy outcome OR 1.80 (95% CI 1.07, 2.96; p=0.022). The DCA can help us address the clinical utility of using family history to predict biopsy outcome.

dca(cancer ~ famhistory, df_binary) %>%
  plot(smooth = TRUE)

First, note that there are many threshold probabilities shown here that are not of interest. For example, it is unlikely that a patient would demand that they had at least a 50% risk of cancer before they would accept a biopsy. Let’s do the DCA again, this time restricting the output to threshold probabilities a more clinically reasonable range, between 0% and 35%.

dca(cancer ~ famhistory, df_binary, thresholds = seq(0, 0.35, by = 0.01)) %>%
  plot(smooth = TRUE)

Now that the graph is showing a more reasonable range of threshold probabilities, let’s assess the clinical utility of family history alone. We can see here that although family history is significantly associated with biopsy outcome, it only adds value to a small range of threshold probabilities near 13% - 20%. If your personal threshold probability is 15% (i.e. you would undergo a biopsy if your probability of cancer was greater than 15%), then family history alone can be beneficial in making the decision to undergo biopsy. However, if your threshold probability is less than 13% or higher than 20%, then family history adds no more benefit than a biopsy all, or biopsy none scheme.

Multivariable Decision Curve Analysis

glm(cancer ~ marker + age + famhistory, df_binary, family = binomial) %>%
  broom::augment(newdata = df_binary, type.predict = "response")
#> # A tibble: 750 × 9
#>    patientid cancer  dead risk_group     age famhistory marker cancerpredmarker
#>        <dbl> <lgl>  <dbl> <chr>        <dbl>      <dbl>  <dbl>            <dbl>
#>  1         1 FALSE      0 low           64.0          0 0.776           0.0372 
#>  2         2 FALSE      0 high          78.5          0 0.267           0.579  
#>  3         3 FALSE      0 low           64.1          0 0.170           0.0216 
#>  4         4 FALSE      0 low           58.5          0 0.0240          0.00391
#>  5         5 FALSE      0 low           64.0          0 0.0709          0.0188 
#>  6         6 FALSE      0 intermediate  65.7          0 0.428           0.0426 
#>  7         7 FALSE      1 intermediate  71.9          0 0.942           0.296  
#>  8         8 TRUE       0 intermediate  66.6          0 2.91            0.387  
#>  9         9 FALSE      0 low           64.3          0 0.715           0.0373 
#> 10        10 FALSE      0 intermediate  65.7          0 0.129           0.0320 
#> # ℹ 740 more rows
#> # ℹ 1 more variable: .fitted <dbl>

dca(cancer ~ famhistory + cancerpredmarker, 
    data = df_binary, 
    thresholds = seq(0, 0.35, by = 0.01)) %>%
  plot(smooth = TRUE)

df_binary_updated <-
  df_binary %>%
  mutate(
    # Use the coefficients from the Brown model
    logodds_Brown = 0.75 * (famhistory) + 0.26 * (age) - 17.5,
    # Convert to predicted probability
    phat_Brown = exp(logodds_Brown) / (1 + exp(logodds_Brown))
  )

# Run the decision curve
dca(cancer ~ phat_Brown, 
    data = df_binary_updated, 
    thresholds = seq(0, 0.35, by = 0.01),
    label = list(phat_Brown = "Brown Model")) %>%
  plot(smooth = TRUE)

df_binary_updated <-
  df_binary_updated %>%
  mutate(
    # Create a variable for the strategy of treating only high risk patients
    # This will be 1 for treat and 0 for don't treat
    high_risk = ifelse(risk_group == "high", 1, 0),
    # Treat based on Joint Approach
    joint = ifelse(risk_group == "high" | cancerpredmarker > 0.15, 1, 0),
    # Treat based on Conditional Approach
    conditional = ifelse(risk_group == "high" | 
                           (risk_group == "intermediate" & cancerpredmarker > 0.15), 1, 0)
  )

dca(cancer ~ high_risk + joint + conditional, 
    data = df_binary_updated, 
    thresholds = seq(0, 0.35, by = 0.01)) %>%
  plot(smooth = TRUE)

Incorporating Harms into Model Assessment

To incorporate the harm of testing and measuring the marker, we ask a clinician, who tells us that, even if the marker were perfectly accurate, few clinicians would conduct more than 30 tests to predict one cancer diagnosis. This might be because the test is expensive, or requires some invasive procedure to obtain. The “harm” of measuring the marker is the reciprocal of 30, or 0.0333. To construct the decision curves for each strategy we now incorporate harm. We have to calculate harm specially for the conditional test, because only patients at intermediate risk are measured for the marker. Then incorporate it into our decision curve. The strategy for incorporating harm for the conditional test is by multiplying the proportion scanned by the harm of the scan.

# the harm of measuring the marker is stored in a scalar
harm_marker <- 0.0333

# in the conditional test, only patients at intermediate risk have their marker measured
intermediate_risk <- df_binary_updated$risk_group == "intermediate"

# harm of the conditional approach is proportion of patients who have the marker measured multiplied by the harm
harm_conditional <- mean(intermediate_risk) * harm_marker

# Run the decision curve 
dca_with_harms <-
  dca(cancer ~ high_risk + joint + conditional, 
    data = df_binary_updated, 
    harm = list(joint = harm_marker, conditional = harm_conditional),
    thresholds = seq(0, 0.35, by = 0.01))
plot(dca_with_harms, smooth = TRUE)

Here the conditional test is clearly the best option (above the 8% treatment threshold), in fact, once you take into account the harm of measuring the marker, it is clearly not worth measuring the marker for everyone: the net benefit of just treating high risk patients is often higher than that of the joint test.

Saving out Net Benefit Values

For any model assessed through decision curve analysis, if we also wanted to show the net benefits in a table, we could save them out. For a particular range of value, we would only need to specify what threshold to start, stop, and the increment we’d like to use.

Let us imagine that we want to view the net benefit of using only the marker to predict whether a patient has cancer, compared with the net benefits of biopsying all patients at thresholds of 5%, 10%, 15%, … 35%. For the model itself, we would actually need to first specify that the marker variable—unlike those of any of the models before—is not a probability. Based on our thresholds, we’d want to begin at 0.05, and by increments of 0.05, stop at 0.35.

dca_with_harms %>%
  as_tibble() %>%
  filter(threshold %in% seq(0.05, 0.35, by = 0.05)) %>%
  select(variable, threshold, net_benefit) %>%
  pivot_wider(id_cols = threshold, 
              names_from = variable,
              values_from = net_benefit)
#> # A tibble: 6 × 6
#>   threshold     all  none high_risk    joint conditional
#>       <dbl>   <dbl> <dbl>     <dbl>    <dbl>       <dbl>
#> 1      0.05  0.0947     0    0.0168  0.0671      0.0740 
#> 2      0.1   0.0444     0    0.0161  0.0571      0.0662 
#> 3      0.2  -0.075      0    0.0147  0.0334      0.0475 
#> 4      0.25 -0.147      0    0.0138  0.0191      0.0362 
#> 5      0.3  -0.229      0    0.0128  0.00289     0.0234 
#> 6      0.35 -0.323      0    0.0116 -0.0159      0.00861

Net benefit has a ready clinical interpretation. The value of 0.03 at a threshold probability of 20% can be interpreted as: “Comparing to conducting no biopsies, biopsying on the basis of the marker is the equivalent of a strategy that found 3 cancers per hundred patients without conducting any unnecessary biopsies.”

Interventions Avoided

As part of assessing the usefulness of the marker, we would be interested in whether using this marker to identify patients with and without cancer would help reduce unnecessary biopsies. This value was the “intervention avoided” column in the table that was saved out. To view it graphically, we would only need to specify it in our command.

dca(cancer ~ marker, 
    data = df_binary, 
    as_probability = "marker") %>%
  net_intervention_avoided() %>%
  plot(smooth = TRUE)

At a probability threshold of 15%, the net reduction in interventions is about 0.33. In other words, at this probability threshold, biopsying patients on the basis of the marker is the equivalent of a strategy that reduced the biopsy rate by 33%, without missing any cancers.

Motivating Example

Patients went to have a marker measured and patients were followed to determine whether they were eventually diagnosed with cancer. We want to build a model using age, family history, and the marker, and assess how well the model predicts cancer within 1.5 years. The variables that will be used from the example dataset for this section of the tutorial are the following:

library(survival)
df_surv %>% head()
#> # A tibble: 6 × 9
#>   patientid cancer ttcancer risk_group    age famhistory marker cancerpredmarker
#>       <dbl> <lgl>     <dbl> <chr>       <dbl>      <dbl>  <dbl>            <dbl>
#> 1         1 FALSE    3.01   low          64.0          0 0.776           0.0372 
#> 2         2 FALSE    0.249  high         78.5          0 0.267           0.579  
#> 3         3 FALSE    1.59   low          64.1          0 0.170           0.0216 
#> 4         4 FALSE    3.46   low          58.5          0 0.0240          0.00391
#> 5         5 FALSE    3.33   low          64.0          0 0.0709          0.0188 
#> 6         6 FALSE    0.0488 intermedia…  65.7          0 0.428           0.0426 
#> # ℹ 1 more variable: cancer_cr <fct>

Basic Data Set-up

The survival probability to any time-point can be derived from any type of survival model; here we use a Cox as this is the most common model in statistical practice. The formula for a survival probability from a Cox model is given by

S(t) = S₀(t)^eXβ

Where X is matrix of covariates in the Cox model, β is a vector containing the parameter estimates from the Cox model, and S₀(t) is the baseline survival probability to time t. To get such values within our code, we will run a Cox model with age, family history, and the marker, as predictors, save out the baseline survival function in a new variable, and obtaining the linear prediction from the model for each subject.

We then obtain the baseline survival probability to our time point of interest. If no patient was observed at the exact time of interest, we can use the baseline survival probability to the observed time closest to, but not after, the time point. We can then calculate the probability of failure at the specified time point. For our example, we will use a time point of 1.5 years.

# build Cox model
cox_model <- coxph(Surv(ttcancer, cancer) ~ age + famhistory + marker, data = df_surv)

# show summary of model results
tbl_regression(cox_model, exponentiate = TRUE)

# add model prediction to our data frame
df_surv_updated <-
  broom::augment(
    cox_model, 
    newdata = df_surv %>% mutate(ttcancer =  1.5),
    type.predict = "expected"
  ) %>%
    mutate(
      pr_failure18 = 1 - exp(-.fitted)
    ) %>%
    select(-.fitted, -.se.fit)

The code for running the decision curve analysis is straightforward after the probability of failure is calculated. All we have to do is specify the time point we are interested in. For our example, let us not only set the threshold from 0% to 50%, but also add smoothing. Results of a smoothed curve should always be compared with the unsmoothed curve to ensure accuracy.

dca(Surv(ttcancer, cancer) ~ pr_failure18, 
    data = df_surv_updated,
    time = 1.5,
    thresholds = 1:50 / 100) %>%
  plot(smooth = TRUE)

This shows that using the model to inform clinical decisions will lead to superior outcomes for any decision associated with a threshold probability of above 2% or so.

Competing Risks

Competing risks endpoints are handled similarly to survival endpoints. The outcome must be defined as a factor, with the lowest level called "censor", and the other levels defining the events of interest. The dca() function will treat the first outcome listed as the outcome of interest.

dca(Surv(ttcancer, cancer_cr) ~ pr_failure18, 
    data = df_surv_updated,
    time = 1.5,
    thresholds = 1:50 / 100) %>%
  plot(smooth = TRUE)