health-economics

What this skill does

# 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

```r
# 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

```r
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

```r
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

```r
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

```r
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

```r
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

```r
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

```r
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

```r
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

```r
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