Calibrating heemod models

The parameters for health economic models can be difficult to measure, either because they cannot be observed directly, or because appropriate data are not systematically gathered in the area of interest. When expected model results are know, model calibration is the search for the appropriate value of initially unknown parameters that allow to obtain these results.

For example the shape and scale parameters of a Weibull survival model can be unknown parameter values. But from the litterature we can know the expected probability of being alive at time t. If this probability is a result from the model, we can find the value of the shape and scale parameters that allow the model results to match, as closely as possible, the observed probability of being alive.

In order to perform calibration, the user must provide:

1. A heemod object from run_model() of update()¹.
2. The names of the parameters of the model to calibrate (the parameters for which we seek appropriate values).
3. A function that when applied to the model returns the result we want to match with reference values.
4. The target values we would like the model results to match.

For this example we will use the result from the assessment of a new total hip replacement previously described in vignette("d-non-homogeneous", "heemod").

We will calibrate the parameters gamma (a Weibull survival parameter) and rrNP1 (the relative risk associated with the new treatment), which originally have values of 1.45 and 0.26 respectively.

The original number of patients with a THR revision after 20 cycles are found in this way:

library(dplyr)
get_counts(res_mod) |> 
  dplyr::filter(model_time == 20 & state_names == "RevisionTHR")

## # A tibble: 2 × 4
##   .strategy_names model_time state_names count
##   <chr>                <int> <chr>       <dbl>
## 1 standard                20 RevisionTHR 2.69 
## 2 np1                     20 RevisionTHR 0.714

We want to calibrate gamma and rrNP1 to obtain 3 patients for the standard strategy and 1 patient for the np1 strategy at time 20. We need to define a function to extract the values we want to change from the model and return them as a numeric vector:

extract_values <- function(x) {
  dplyr::filter(
    get_counts(x),
    model_time == 20 & state_names == "RevisionTHR"
  )$count
}
extract_values(res_mod)

## [1] 2.687124 0.714282

Any arbitrary function of any model output would work, as long as it returns numeric values.

A convenience function define_calibration_fn() exists to help easily define calibration functions.

calib_fn <- define_calibration_fn(
  type = "count",
  strategy_names = c("standard", "np1"),
  element_names = c("RevisionTHR", "RevisionTHR"),
  cycles = c(20, 20)
)
calib_fn(res_mod)

## [1] 2.687124 0.714282

We can now call calibrate_model(), and give the values we want to reach as target_values.

res_cal <- calibrate_model(
  res_mod,
  parameter_names = c("gamma", "rrNP1"),
  fn_values = extract_values,
  target_values = c(2.5, 0.8)
)

## Loading required namespace: optimx

res_cal

##      gamma     rrNP1        value convcode
## 1 1.431919 0.3146302 3.346417e-10        0

The new parameter values are 1.43 for gamma and 0.31 for rrNP1. The convcode code at 0 indicates the calibration was successful.

It is possible to specify several possible starting values for the calibration procedure in order to explore the parameter space:

start <- data.frame(
  gamma = c(1.0, 1.5, 2.0),
  rrNP1 = c(0.2, 0.3, 0.4)
)

res_cal_2 <- calibrate_model(
  res_mod,
  parameter_names = c("gamma", "rrNP1"),
  fn_values = extract_values,
  target_values = c(3, 1),
  initial_values = start,
  lower = c(0, 0), upper = c(2, 1)
)

Additional options to control the optimization process can be passed to calibrate_model(). These options are parameters of the optimx function, such as upper and lower to specify upper and lower values, method to change the optimization method, etc.

Calibration uses optimization to minimize the sum of squared errors between calculated and desired values, and so is subject to all the many difficulties of optimization. Different optimization methods, for example Nelder-Mead (which does not require gradients) and BFGS or conjugate gradient methods, which do require gradients but can approximate them numerically, may work better for different problems. Some attempted optimizations may not converge.

It may be impossible to evaluate the function at badly-specified initial parameter values (for example, if a negative initial value is given for a parameter that must be positive); using box limits on some parameters may help with this.

Even if the calibration converges from different initial values, it may not converge to the same parameter values every time; in general, an underconstrained model can have different parameter sets that fit equally well. For these and other reasons, the user is advised to carefully check the results of calibrations.