health-economics

name	health-economics
description	Health economic analysis in R, including cost-effectiveness, QALYs, decision models, and budget impact.

Health Economics Evaluation in R

Overview

Health economic evaluation methods covering cost-effectiveness analysis (CEA), quality-adjusted life years (QALYs), incremental cost-effectiveness ratios (ICERs), budget impact analysis, Markov cohort models, partitioned survival analysis, probabilistic sensitivity analysis, and value of information analysis.

Cost-Effectiveness Fundamentals

Basic Calculations

# Treatment comparison data
# Intervention vs Comparator
costs_int <- 15000      # Mean cost of intervention
costs_comp <- 8000      # Mean cost of comparator
effects_int <- 5.2      # Mean QALYs intervention
effects_comp <- 4.5     # Mean QALYs comparator

# Incremental calculations
delta_cost <- costs_int - costs_comp      # Incremental cost
delta_effect <- effects_int - effects_comp # Incremental effect (QALYs)

# Incremental Cost-Effectiveness Ratio (ICER)
icer <- delta_cost / delta_effect
cat("ICER:", round(icer, 0), "per QALY gained\n")

# Net Monetary Benefit (NMB) at WTP threshold
wtp <- 50000  # Willingness-to-pay threshold
nmb <- delta_effect * wtp - delta_cost
cat("NMB at WTP $", wtp, ":", round(nmb, 0), "\n")

# Net Health Benefit (NHB)
nhb <- delta_effect - delta_cost / wtp
cat("NHB:", round(nhb, 3), "QALYs\n")

Cost-Effectiveness Plane

library(ggplot2)

# Simulated incremental costs and effects
set.seed(123)
n_sim <- 1000
delta_c <- rnorm(n_sim, delta_cost, 2000)
delta_e <- rnorm(n_sim, delta_effect, 0.3)

ce_data <- data.frame(
  delta_cost = delta_c,
  delta_effect = delta_e
)

# CE plane
ggplot(ce_data, aes(x = delta_effect, y = delta_cost)) +
  geom_point(alpha = 0.3, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_abline(slope = wtp, intercept = 0, color = "red", linetype = "dashed") +
  annotate("text", x = 1.5, y = 15000, label = paste0("WTP = $", wtp),
           color = "red") +
  labs(x = "Incremental Effect (QALYs)",
       y = "Incremental Cost ($)",
       title = "Cost-Effectiveness Plane") +
  theme_bw()

Cost-Effectiveness Analysis with BCEA

Using BCEA Package

library(BCEA)

# From PSA samples (effects and costs matrices)
# Each row = simulation, columns = interventions
n_sim <- 1000
n_int <- 2  # Number of interventions

# Effects matrix (QALYs)
effects <- matrix(
  c(rnorm(n_sim, 4.5, 0.5),   # Comparator
    rnorm(n_sim, 5.2, 0.6)),  # Intervention
  nrow = n_sim, ncol = n_int
)

# Costs matrix
costs <- matrix(
  c(rnorm(n_sim, 8000, 1500),   # Comparator
    rnorm(n_sim, 15000, 3000)), # Intervention
  nrow = n_sim, ncol = n_int
)

colnames(effects) <- colnames(costs) <- c("Comparator", "Intervention")

# Create BCEA object
bcea_result <- bcea(
  e = effects,
  c = costs,
  ref = 1,                            # Reference intervention
  interventions = c("Comparator", "Intervention"),
  Kmax = 100000                       # Max WTP for analysis
)

# Summary at specific WTP
summary(bcea_result, wtp = 50000)

# Key outputs
bcea_result$ICER                      # ICER
bcea_result$ceac                      # CEAC values

CE Plane and CEAC Plots

library(BCEA)

# Cost-effectiveness plane
ceplane.plot(bcea_result,
             wtp = 50000,
             graph = "ggplot2",
             title = "Cost-Effectiveness Plane")

# Cost-Effectiveness Acceptability Curve (CEAC)
ceac.plot(bcea_result,
          graph = "ggplot2",
          title = "Cost-Effectiveness Acceptability Curve")

# Cost-Effectiveness Acceptability Frontier (CEAF)
ceaf.plot(bcea_result, graph = "ggplot2")

# Expected Incremental Benefit (EIB) plot
eib.plot(bcea_result, graph = "ggplot2")

Multiple Interventions

library(BCEA)

# Three interventions
effects_3 <- matrix(
  c(rnorm(n_sim, 4.0, 0.4),   # Standard care
    rnorm(n_sim, 4.8, 0.5),   # Treatment A
    rnorm(n_sim, 5.5, 0.6)),  # Treatment B
  nrow = n_sim
)

costs_3 <- matrix(
  c(rnorm(n_sim, 5000, 1000),
    rnorm(n_sim, 12000, 2500),
    rnorm(n_sim, 20000, 4000)),
  nrow = n_sim
)

colnames(effects_3) <- colnames(costs_3) <- c("Standard", "Treatment_A", "Treatment_B")

bcea_multi <- bcea(
  e = effects_3,
  c = costs_3,
  ref = 1,
  interventions = colnames(effects_3)
)

# Multi-comparison CEAC
mce <- multi.ce(bcea_multi)
ceac.plot(mce, graph = "ggplot2")

# Contour plot
contour2(bcea_multi, wtp = 50000)

Markov Cohort Models

Using heemod Package

library(heemod)

# Define transition probabilities
mat_trans <- define_transition(
  state_names = c("Healthy", "Sick", "Dead"),
  # From Healthy
  C, 0.15, 0.01,
  # From Sick
  0.10, C, 0.05,
  # From Dead (absorbing)
  0, 0, 1
)

# Define states with costs and utilities
state_healthy <- define_state(
  cost = 0,
  utility = 1
)

state_sick <- define_state(
  cost = 5000,
  utility = 0.7
)

state_dead <- define_state(
  cost = 0,
  utility = 0
)

# Define strategy (no treatment)
strat_base <- define_strategy(
  transition = mat_trans,
  Healthy = state_healthy,
  Sick = state_sick,
  Dead = state_dead
)

# Run the model
result_base <- run_model(
  strat_base,
  cycles = 50,
  cost = cost,
  effect = utility,
  init = c(1000, 0, 0),     # Initial cohort distribution
  method = "beginning"       # Cycle correction
)

# Summary
summary(result_base)
plot(result_base)

Comparing Strategies

library(heemod)

# Define treatment strategy (reduced transition to sick)
mat_trans_trt <- define_transition(
  state_names = c("Healthy", "Sick", "Dead"),
  C, 0.10, 0.01,     # Lower transition to sick
  0.15, C, 0.04,     # Higher recovery
  0, 0, 1
)

state_healthy_trt <- define_state(
  cost = 500,        # Treatment cost
  utility = 1
)

strat_trt <- define_strategy(
  transition = mat_trans_trt,
  Healthy = state_healthy_trt,
  Sick = state_sick,
  Dead = state_dead
)

# Run both strategies
result_comp <- run_model(
  base = strat_base,
  treatment = strat_trt,
  cycles = 50,
  cost = cost,
  effect = utility,
  init = c(1000, 0, 0)
)

# Summary comparison
summary(result_comp)

# ICER
icer_result <- summary(result_comp)$res_comp
print(icer_result)

Time-Dependent Parameters

library(heemod)

# Parameters that change over time
param <- define_parameters(
  age_init = 50,
  age = age_init + model_time,
  mortality = 1 - exp(-0.0001 * age^2),  # Age-dependent mortality
  p_sick = 0.1 + 0.005 * model_time       # Increasing disease risk
)

# Transition matrix with parameters
mat_time <- define_transition(
  state_names = c("Healthy", "Sick", "Dead"),
  C, p_sick, mortality,
  0.05, C, mortality * 1.5,
  0, 0, 1
)

# Run with parameters
result_time <- run_model(
  define_strategy(
    transition = mat_time,
    Healthy = state_healthy,
    Sick = state_sick,
    Dead = state_dead
  ),
  cycles = 30,
  cost = cost,
  effect = utility,
  init = c(1000, 0, 0),
  parameters = param
)

Partitioned Survival Analysis

Using hesim Package

library(hesim)
library(flexsurv)

# Fit parametric survival models for each health state
# Overall Survival (OS)
fit_os <- flexsurvreg(
  Surv(time, status) ~ treatment,
  data = surv_data,
  dist = "weibull"
)

# Progression-Free Survival (PFS)
fit_pfs <- flexsurvreg(
  Surv(time_pfs, status_pfs) ~ treatment,
  data = surv_data,
  dist = "weibull"
)

# State probabilities from survival curves
# Pre-progression: S_PFS(t)
# Post-progression: S_OS(t) - S_PFS(t)
# Death: 1 - S_OS(t)

Building PSM with hesim

library(hesim)

# Treatment strategies
strategies <- data.table(
  strategy_id = 1:2,
  strategy_name = c("Standard", "New Treatment")
)

# Patients (can include heterogeneity)
patients <- data.table(
  patient_id = 1:100,
  age = rnorm(100, 60, 10)
)

# Health states
states <- data.table(
  state_id = 1:3,
  state_name = c("Stable", "Progressed", "Dead")
)

# Create hesim data object
hesim_data <- hesim_data(
  strategies = strategies,
  patients = patients,
  states = states
)

# Define input data for survival models
input_data <- expand(hesim_data, by = c("strategies", "patients"))

# Create survival curves (from fitted models)
survmods <- create_PsmCurves(
  input_data = input_data,
  params = params_surv_list  # From fitted flexsurv models
)

# State probabilities
stprobs <- survmods$sim_stateprobs(t = seq(0, 10, 0.1))

Costs and QALYs from PSM

library(hesim)

# State utilities
utility_tbl <- stateval_tbl(
  data.table(
    state_id = 1:3,
    est = c(0.8, 0.6, 0)  # Utilities by state
  ),
  dist = "fixed"
)

# State costs
cost_tbl <- stateval_tbl(
  data.table(
    state_id = 1:3,
    est = c(1000, 5000, 0)  # Annual costs
  ),
  dist = "fixed"
)

# Create PSM object
psm <- Psm$new(
  survival_models = survmods,
  utility_model = utility_tbl,
  cost_models = list(drug = cost_tbl)
)

# Simulate QALYs and costs
psm$sim_stateprobs(t = seq(0, 10, by = 0.1))
psm$sim_qalys(dr = 0.03)  # 3% discount rate
psm$sim_costs(dr = 0.03)

# Summarize
ce_results <- psm$summarize()

Probabilistic Sensitivity Analysis

Parameter Distributions

library(heemod)

# Define parameter distributions
param_psa <- define_parameters(
  # Beta distribution for probabilities
  p_sick = rbeta(1, shape1 = 20, shape2 = 80),
  p_death_healthy = rbeta(1, shape1 = 2, shape2 = 198),
  p_death_sick = rbeta(1, shape1 = 10, shape2 = 90),

  # Gamma distribution for costs
  cost_sick = rgamma(1, shape = 100, rate = 0.02),
  cost_treatment = rgamma(1, shape = 50, rate = 0.1),

  # Beta distribution for utilities
  utility_sick = rbeta(1, shape1 = 70, shape2 = 30)
)

# Run PSA
psa_result <- run_psa(
  model = result_comp,
  psa = param_psa,
  N = 1000  # Number of simulations
)

# Summary
summary(psa_result)

# CE plane from PSA
plot(psa_result, type = "ce")

# CEAC from PSA
plot(psa_result, type = "ac")

Custom PSA Implementation

# Manual PSA loop
n_sim <- 1000
psa_results <- data.frame(
  sim = 1:n_sim,
  cost_base = NA,
  cost_trt = NA,
  qaly_base = NA,
  qaly_trt = NA
)

for (i in 1:n_sim) {
  # Sample parameters
  p_sick <- rbeta(1, 20, 80)
  cost_sick <- rgamma(1, 100, 0.02)
  utility_sick <- rbeta(1, 70, 30)

  # Run model with sampled parameters
  # ... model calculations ...

  # Store results
  psa_results$cost_base[i] <- sum_cost_base