Basic examples

Introduction

This vignette provides a walkthrough of the basic functionality of the {blendR} package. We’ll show how to combine survival models from different sources, such as observed trial data and external evidence, using various fitting methods.

library(blendR)

library(survHE)
library(INLA)

Example 1: INLA piece-wise exponential and external knowledge

This first example demonstrates a common use case: blending a semi-parametric model fitted to observed data with a parametric model representing external knowledge.

Prepare the models

First, we load the package’s example dataset and inspect its structure. This dataset contains individual patient data from a clinical trial.

data("TA174_FCR", package = "blendR")
head(dat_FCR)
#> # A tibble: 6 × 5
#>   patid treat death death_t death_ty
#>   <int> <int> <int>   <dbl>    <dbl>
#> 1     1     1     0  32       2.67  
#> 2     2     1     0  30.6     2.55  
#> 3     3     1     0  28       2.33  
#> 4     8     1     0  30       2.5   
#> 5    10     1     1   0.458   0.0382
#> 6    11     1     1   1.57    0.131

Next, we fit a piece-wise exponential model to the observed trial data. This is done using Integrated Nested Laplace Approximation (INLA) via the fit_inla_pw() helper function. We define the cut-points for the piece-wise hazard function.

# observed estimate
obs_Surv <- fit_inla_pw(data = dat_FCR,
                        cutpoints = seq(0, 180, by = 5))

For the external model, we first generate a synthetic dataset that is consistent with some external information or expert opinion. Here, we simulate data where the survival probability is known to be 5% at 144 months, up to a maximum follow-up time of 180 months.

# external estimate
data_sim <- ext_surv_sim(t_info = 144,
                         S_info = 0.05,
                         T_max = 180)

Now, we fit a parametric Gompertz model to this synthetic external data. This is done using Hamiltonian Monte Carlo (HMC) via the {survHE} package. Note that this step requires the {survHEhmc} package to be installed.

ext_Surv <- fit.models(formula = Surv(time, event) ~ 1,
                       data = data_sim,
                       distr = "gompertz",
                       method = "hmc",
                       priors = list(gom = list(a_alpha = 0.1,
                                                b_alpha = 0.1)))
#> The survHEhmc version loaded is: 0.0.11
#> 
#> SAMPLING FOR MODEL 'Gompertz' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 3.8e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.38 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1.925 seconds (Warm-up)
#> Chain 1:                0.232 seconds (Sampling)
#> Chain 1:                2.157 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'Gompertz' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 2.5e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.362 seconds (Warm-up)
#> Chain 2:                0.239 seconds (Sampling)
#> Chain 2:                0.601 seconds (Total)
#> Chain 2:

Blend the curves

With both the observed and external survival models fitted, we can now use the core blendsurv() function. We must specify the blending interval (blend_interv) and the parameters of the Beta distribution (beta_params) that will control the blending weight.

blend_interv <- list(min = 48, max = 150)
beta_params <- list(alpha = 3, beta = 3)

ble_Surv <- blendsurv(obs_Surv, ext_Surv, blend_interv, beta_params)

Finally, we can easily visualize the observed curve, the external curve, and the final blended curve using the default plot method.

plot(ble_Surv)
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the blendR package.
#>   Please report the issue at
#>   <https://github.com/StatisticsHealthEconomics/blendR/issues/>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Example 2: Blending two HMC survHE models

This example shows how to blend two survival curves that were both fitted using HMC with the {survHE} package.

We’ll fit an exponential model to both the observed trial data and the synthetic external data created in the previous example.

obs_Surv2 <- fit.models(formula = Surv(death_t, death) ~ 1,
                        data = dat_FCR,
                        distr = "exponential",
                        method = "hmc")
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 9.2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.92 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.208 seconds (Warm-up)
#> Chain 1:                0.142 seconds (Sampling)
#> Chain 1:                0.35 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 3.2e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.32 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.225 seconds (Warm-up)
#> Chain 2:                0.142 seconds (Sampling)
#> Chain 2:                0.367 seconds (Total)
#> Chain 2:

ext_Surv2 <- fit.models(formula = Surv(time, event) ~ 1,
                       data = data_sim,
                       distr = "exponential",
                       method = "hmc")
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 2e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.087 seconds (Warm-up)
#> Chain 1:                0.048 seconds (Sampling)
#> Chain 1:                0.135 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 1e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.079 seconds (Warm-up)
#> Chain 2:                0.048 seconds (Sampling)
#> Chain 2:                0.127 seconds (Total)
#> Chain 2:

The blending step is identical to before, demonstrating the consistent interface. We can then plot the result and add a non-parametric Kaplan-Meier estimate for comparison.

ble_Surv2 <- blendsurv(obs_Surv2, ext_Surv2, blend_interv, beta_params)

# kaplan-meier
km <- survfit(Surv(death_t, death) ~ 1, data = dat_FCR)

plot(ble_Surv2) +
  geom_line(aes(km$time, km$surv, colour = "Kaplan-Meier"),
            size = 1.25, linetype = "dashed")

Example 3: Blending HMC and frequentist models

This final example highlights the flexibility of {blendR} by blending a Bayesian model (from HMC) with a frequentist one fitted using the {flexsurv} package.

First, we fit the observed data model using HMC via {survHE}. Then, we fit the external data model using flexsurv::flexsurvreg().

obs_Surv3 <- fit.models(formula = Surv(death_t, death) ~ 1,
                        data = dat_FCR,
                        distr = "exponential",
                        method = "hmc")
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 4e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.209 seconds (Warm-up)
#> Chain 1:                0.162 seconds (Sampling)
#> Chain 1:                0.371 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'Exponential' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 3.1e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.31 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.202 seconds (Warm-up)
#> Chain 2:                0.125 seconds (Sampling)
#> Chain 2:                0.327 seconds (Total)
#> Chain 2:

ext_Surv3 <- flexsurv::flexsurvreg(formula = Surv(time, event) ~ 1,
                       data = data_sim,
                       dist = "gompertz")

ble_Surv3 <- blendsurv(obs_Surv3, ext_Surv3, blend_interv, beta_params)

plot(ble_Surv3)