In W-Mohammed/calibrater: Calibration of Decision-Analytic Models

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(devtools)
load_all()

\pagebreak

Introduction:

Model: Infectious disease model

The model (Nicolas et al.,, 2017) is adapted from approaches for modelling HIV in high-burden settings. The population is divided into five health states including non-susceptible (N), susceptible (S), early disease (E), late disease (L), and treatment (T). The number of individuals by state and year (t) is given by $N_t$, $S_t$, $E_t$, $L_t$, and $T_t$ respectively. Individuals enter the model distributed across the $N$ and $S$ states, and transition between states to allow for infection ($S$ to $E$), disease progression ($E$ to $L$), treatment initiation ($E$ and $L$ to $T$), and death ($N$, $S$, $E$, $L$ and $T$ to $D$) via background and disease-specific mortality. The diagram below represents the model.

HIV model

Calibration process specifications:

Calibration parameters:

mu_e Cause-specific mortality rate with early-stage disease
mu_l Cause-specific mortality rate with late-stage disease
mu_t Cause-specific mortality rate on treatment
p Transition rate from early to late-stage disease
r_l Rate of uptake onto treatment (r_l = late-stage disease)
rho Effective contact rate
b Fraction of population in at-risk group

Calibration targets:

Prev Prevalence at 10, 20, 30 years
Surv HIV survival without treatment
Trt_vol Treatment volume at 30 years

Parameter transformations:

Two calibration process was carried out twice, with and without parameter-transformation.

The calibration parameters were transformed as follows: - b was transformed to the entire real line $(-\infty, \infty)$ using logit. - All other parameters were transformed using log.

Search method:

Random search using:
Full factorial grid
Random grid
Latin-Hypercube Sampling
Directed search using:
Gradient-based (GB) (uses derivatives)
Nelder-Mead (NM) algorithm, aka simplex method (derivative free)
Global optimization techniques:
- Simulated Annealing (SANN)
- Genetic algorithms (GA)

Goodness-of-fit measure:

Sum of log-likelihoods
Sum of squared errors

Bayesian methods:

Sampling Importance Resampling (SIR)
Incremental Mixture Importance Sampling (IMIS)

The Calibration process:

The code blocks below show the code for the untransformed version of the model; while the results from both runs will be displayed below in the results sub-section.

Sampling from prior using different methods:

set.seed(1)
HID_results <- list()
# LHS:----
HID_results$Prior_samples[['LHS']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)
# FGS:----
HID_results$Prior_samples[['FGS']] <- sample_prior_FGS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)
# RGS:----
HID_results$Prior_samples[['RGS']] <- sample_prior_RGS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)

A quick look at a sample of the samples:

Values from both runs of the process are shown below.

# Original scale parameters:----
HID_results$Prior_samples[['LHS']] %>%
  head(5)
# Transformed scale parameters:----
HID2_results$Prior_samples[['LHS']] %>%
  head(5)

Using random-search with goodness-of-fit (without optimisation):

While any of the three random search methods can be used, the following code uses LHS.

# LHS with weighted sum of square errors:----
HID_results$Calib_results$Random[[1]] <- wSSE_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = HID_results$Prior_samples$LHS,
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")
# LHS with log likelihood:----
HID_results$Calib_results$Random[[2]] <- LLK_GOF(
  .func = HID_markov, .optim = FALSE,
  .args = NULL,
  .samples = HID_results$Prior_samples$LHS,
  .l_targets = HID_data$l_targets,
  .sample_method = "LHS")

A quick look at a sample of the calibration results:

Values from runs using original and scaled parameters are shown below.

# O_scale: LHS with weighted sum of square errors:----
HID_results$Calib_results$Random[[1]] %>%
  head(5)
# O_scale: LHS with log likelihood:----
HID_results$Calib_results$Random[[2]] %>%
  head(5)
# T_scale: LHS with log likelihood:----
HID2_results$Calib_results$Random[[1]] %>%
  head(5)
# T_scale: LHS with log likelihood:----
HID2_results$Calib_results$Random[[2]] %>%
  head(5)

Optimisation algorithms:

Both goodness-of-fit methods are used in the code below.

# Sample 5 starting parameter sets for the optimisation algorithms:----
HID_results$Prior_samples[['LHS_Directed']] <- sample_prior_LHS(
  .n_samples = 5,
  .l_params = HID_data$l_params)
# Nelder-Mead using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[1]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'log_likelihood',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'Nelder-Mead',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)
# Nelder-Mead using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[2]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'wSumSquareError',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'Nelder-Mead',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)
# Gradient-based using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[3]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'log_likelihood',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)
# Gradient-based using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[4]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'wSumSquareError',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'BFGS',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000)
# Simulated annealing using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[5]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'log_likelihood',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  fnscale = -1,
  temp = 10,
  tmax = 10,
  maxit = 1000)
# Simulated annealing using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[6]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'wSumSquareError',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'SANN',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)
# Genetic algorithm using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[7]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'log_likelihood',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)
# Genetic algorithm using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[8]] <- calibrateModel_directed(
  .l_params = HID_data$l_params,
  .func = HID_markov,
  .args = NULL,
  .gof = 'wSumSquareError',
  .samples = HID_results$Prior_samples$LHS_Directed,
  .s_method = 'GA',
  .maximise = TRUE,
  .l_targets = HID_data$l_targets,
  maxit = 1000,
  temp = 10,
  tmax = 10)

A quick look at a sample of the calibration results:

# Original scale:----
## O: Nelder-Mead using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[1]][[1]]
## O: Nelder-Mead using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[2]][[1]]
## O: Gradient-based using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[3]][[1]]
## O: Gradient-based using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[4]][[1]]
## O: Simulated annealing using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[5]][[1]]
## O: Simulated annealing using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[6]][[1]]
## O: Genetic algorithm using log_likelihood as goodness-of-fit:----
HID_results$Calib_results$Directed[[7]][[1]]
## O: Genetic algorithm using sum-squared-errors as goodness-of-fit:----
HID_results$Calib_results$Directed[[8]][[1]]

# Transformed scale:----
## T: Nelder-Mead using log_likelihood as goodness-of-fit:----
HID2_results$Calib_results$Directed[[1]][[1]]
## T: Nelder-Mead using sum-squared-errors as goodness-of-fit:----
HID2_results$Calib_results$Directed[[2]][[1]]
## T: Gradient-based using log_likelihood as goodness-of-fit:----
HID2_results$Calib_results$Directed[[3]][[1]]
## T: Gradient-based using sum-squared-errors as goodness-of-fit:----
HID2_results$Calib_results$Directed[[4]][[1]]
## T: Simulated annealing using log_likelihood as goodness-of-fit:----
HID2_results$Calib_results$Directed[[5]][[1]]
## T: Simulated annealing using sum-squared-errors as goodness-of-fit:----
HID2_results$Calib_results$Directed[[6]][[1]]
## T: Genetic algorithm using log_likelihood as goodness-of-fit:----
HID2_results$Calib_results$Directed[[7]][[1]]
## T: Genetic algorithm using sum-squared-errors as goodness-of-fit:----
HID2_results$Calib_results$Directed[[8]][[1]]

Bayesian methods:

# Sample values for the SIR methods:----
HID_results$Prior_samples[['LHS_Bayesian']] <- sample_prior_LHS(
  .n_samples = 10000,
  .l_params = HID_data$l_params)
# SIR:----
HID_results$Calib_results$Bayesian[[1]] = calibrateModel_beyesian(
  .b_method = 'SIR', .func = HID_markov,
  .args = NULL,
  .l_targets = HID_data$l_targets,
  .n_resample = 10000,
  .l_params = HID_data$l_params,
  .samples = HID_results$Prior_samples$LHS_Bayesian)
# IMIS:----
HID_results$Calib_results$Bayesian[[2]] = calibrateModel_beyesian(
  .b_method = 'IMIS', .func = HID_markov_2,
  .args = NULL,
  .l_targets = HID_data2$l_targets,
  .l_params = HID_data2$l_params,
  .transform = TRUE,
  .n_resample = 1000,
  .IMIS_iterations = 10,
  .IMIS_sample = 100)

A quick look at a sample of the calibration results:

# Original scale:----
## Effective sample size:----
### O: SIR:----
HID_results$Calib_results$Bayesian[[1]] %>%
  effective_sample_size(bayes_calib_output_list = .)
### O: IMIS:----
HID_results$Calib_results$Bayesian[[2]] %>%
  effective_sample_size(bayes_calib_output_list = .)
## Number of unique parameter sets:----
### O: SIR:----
HID_results$Calib_results$Bayesian[[1]]$Results %>%
  distinct() %>%
  nrow()
### O: IMIS:----
HID_results$Calib_results$Bayesian[[2]]$Results %>%
  distinct() %>%
  nrow()

# Transformed scale:----
## Effective sample size:----
### T: SIR:----
HID2_results$Calib_results$Bayesian[[1]] %>%
  effective_sample_size(bayes_calib_output_list = .)
### T: IMIS:----
HID2_results$Calib_results$Bayesian[[2]] %>%
  effective_sample_size(bayes_calib_output_list = .)
## Number of unique parameter sets:----
### T: SIR:----
HID2_results$Calib_results$Bayesian[[1]]$Results %>%
  distinct() %>%
  nrow()
### T: IMIS:----
HID2_results$Calib_results$Bayesian[[2]]$Results %>%
  distinct() %>%
  nrow()

PSA:

Sample PSA parameter draws using directed methods results where possible:

In addition to sampling PSA parameter draws, the following function ensures all results are in proper shape for the next step. Therefore, results from other methods are passed to the same function.

# Passing calibration results from directed search:
HID_results$PSA_samples[["Directed"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Directed,
  .search_method = 'Directed',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)
# Passing calibration results from random search:
HID_results$PSA_samples[["Random"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Random,
  .search_method = 'Random',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)
# Passing calibration results from Bayesian methods:
HID_results$PSA_samples[["Bayesian"]] <- PSA_calib_values(
  .l_calib_res_lists = HID_results$Calib_results$Bayesian,
  .search_method = 'Bayesian',
  .PSA_samples = 10000,
  .transform_ = FALSE,
  .l_params = HID_data$l_params)

Run PSA:

The following function runs PSA draws from each calibration process separately.

# Run all calibration results together:
HID_results$PSA_results <- run_PSA(
  .func_ = HID_markov,
  .PSA_calib_values_ = c(HID_results$PSA_samples$Directed,
                         HID_results$PSA_samples$Random,
                         HID_results$PSA_samples$Bayesian),
  .args_ = list(calibrate_ = FALSE),
  .PSA_unCalib_values_ = NULL)

Results:

In addition to the incremental net benefit (at £30,000), the values of incremental cost and life years gained are shown below.

# Original scale:----
HID_results$PSA_summary <- 
  map_df(
    .x = HID_results$PSA_results,
    .f = function(PSA) {
      data_ <- tibble(
        'mean_inc_Costs' = mean(PSA$inc_cost),
        'mean_inc_LY' = mean(PSA$inc_LY),
        'iNMB' = (mean_inc_LY * 30000) - mean_inc_Costs,
        'calibration_method' = if(nrow(PSA) == 1) paste(PSA$Label[[1]], "_*") else PSA$Label[[1]],
        'goodness_of_fit' = PSA$Overall_fit[[1]]
      )
    } 
  )

HID_results$PSA_summary

# Transformed scale:----
HID2_results$PSA_summary <- 
  map_df(
    .x = HID2_results$PSA_results,
    .f = function(PSA) {
      data_ <- tibble(
        'mean_inc_Costs' = mean(PSA$inc_cost),
        'mean_inc_LY' = mean(PSA$inc_LY),
        'iNMB' = (mean_inc_LY * 30000) - mean_inc_Costs,
        'calibration_method' = if(nrow(PSA) == 1) paste(PSA$Label[[1]], "_*") else PSA$Label[[1]],
        'goodness_of_fit' = PSA$Overall_fit[[1]]
      )
    } 
  )

HID2_results$PSA_summary

W-Mohammed/calibrater documentation built on Oct. 14, 2023, 1:57 a.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

W-Mohammed/calibrater
Calibration of Decision-Analytic Models

In W-Mohammed/calibrater: Calibration of Decision-Analytic Models

Introduction:

Model: Infectious disease model

Calibration process specifications:

Calibration parameters:

Calibration targets:

Parameter transformations:

Search method:

Goodness-of-fit measure:

Bayesian methods:

The Calibration process:

Sampling from prior using different methods:

A quick look at a sample of the samples:

Using random-search with goodness-of-fit (without optimisation):

A quick look at a sample of the calibration results:

Optimisation algorithms:

A quick look at a sample of the calibration results:

Bayesian methods:

A quick look at a sample of the calibration results:

PSA:

Sample PSA parameter draws using directed methods results where possible:

Run PSA:

Results:

R Package Documentation

Browse R Packages

We want your feedback!

W-Mohammed/calibrater Calibration of Decision-Analytic Models

In W-Mohammed/calibrater: Calibration of Decision-Analytic Models

Introduction:

Model: Infectious disease model

Calibration process specifications:

Calibration parameters:

Calibration targets:

Parameter transformations:

Search method:

Goodness-of-fit measure:

Bayesian methods:

The Calibration process:

Sampling from prior using different methods:

A quick look at a sample of the samples:

Using random-search with goodness-of-fit (without optimisation):

A quick look at a sample of the calibration results:

Optimisation algorithms:

A quick look at a sample of the calibration results:

Bayesian methods:

A quick look at a sample of the calibration results:

PSA:

Sample PSA parameter draws using directed methods results where possible:

Run PSA:

Results:

R Package Documentation

Browse R Packages

We want your feedback!

W-Mohammed/calibrater
Calibration of Decision-Analytic Models