Foundational knowledge for writing current PyMC models including syntax, distributions, sampling, and ArviZ diagnostics. Use when creating or reviewing PyMC models.

2026-05-277

stan-fundamentals

choxos/BiostatAgent

Foundational knowledge for writing modern Stan models including program structure, type system, distributions, and best practices. Use when creating or reviewing Stan models.

2026-05-277

group-sequential-methods

choxos/BiostatAgent

Group sequential design methods for interim analyses, alpha spending, and futility stopping. Use when designing trials with interim looks or implementing spending functions.

2026-05-277

mediana-fundamentals

choxos/BiostatAgent

Core Mediana package functions for Clinical Scenario Evaluation (CSE). Use when designing data models, analysis models, evaluation models, and running comprehensive trial simulations.

2026-05-277

simtrial-fundamentals

choxos/BiostatAgent

Core simtrial package functions for time-to-event clinical trial simulation. Use when generating survival data, performing weighted logrank tests, or running TTE simulations.

2026-05-277

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
  psa_results$cost_trt[i] <- sum_cost_trt
  psa_results$qaly_base[i] <- sum_qaly_base
  psa_results$qaly_trt[i] <- sum_qaly_trt
}

# Incremental results
psa_results <- psa_results |>
  mutate(
    delta_cost = cost_trt - cost_base,
    delta_qaly = qaly_trt - qaly_base,
    icer = delta_cost / delta_qaly
  )

# CEAC calculation
wtp_range <- seq(0, 100000, 1000)
ceac <- sapply(wtp_range, function(k) {
  mean(psa_results$delta_qaly * k - psa_results$delta_cost > 0)
})

# Plot CEAC
plot(wtp_range, ceac, type = "l",
     xlab = "Willingness-to-Pay ($/QALY)",
     ylab = "Probability Cost-Effective",
     main = "Cost-Effectiveness Acceptability Curve")

Value of Information Analysis

Expected Value of Perfect Information (EVPI)

library(BCEA)

# EVPI from BCEA object
evpi_result <- evpi(bcea_result)

# Plot EVPI
evi.plot(bcea_result, graph = "ggplot2")

# EVPI at specific WTP
evpi_50k <- evpi_result[evpi_result$k == 50000, "evpi"]
cat("EVPI at $50,000/QALY:", round(evpi_50k, 0), "\n")

# Population EVPI
pop_size <- 100000           # Affected population
time_horizon <- 10           # Years of technology relevance
discount_rate <- 0.03        # Annual discount rate

# Annuity factor
annuity <- (1 - (1 + discount_rate)^(-time_horizon)) / discount_rate

pop_evpi <- evpi_50k * pop_size * annuity
cat("Population EVPI:", round(pop_evpi / 1e6, 2), "million\n")

Expected Value of Partial Perfect Information (EVPPI)

library(BCEA)

# EVPPI for specific parameters
# Requires PSA parameter matrix
psa_params <- data.frame(
  p_sick = rbeta(n_sim, 20, 80),
  cost_sick = rgamma(n_sim, 100, 0.02),
  utility_sick = rbeta(n_sim, 70, 30)
)

# EVPPI for p_sick
evppi_result <- evppi(
  bcea_result,
  param_idx = "p_sick",
  input = psa_params,
  method = "gam"
)

# Plot
plot(evppi_result, graph = "ggplot2")

# EVPPI for multiple parameters
evppi_multi <- evppi(
  bcea_result,
  param_idx = c("p_sick", "cost_sick"),
  input = psa_params
)

Budget Impact Analysis

Basic Budget Impact Model

# Population and market share assumptions
pop_eligible <- 50000        # Eligible patient population
uptake_year1 <- 0.10
uptake_year2 <- 0.25
uptake_year3 <- 0.40
uptake_year4 <- 0.50
uptake_year5 <- 0.55

uptake <- c(uptake_year1, uptake_year2, uptake_year3, uptake_year4, uptake_year5)

# Annual treatment costs
cost_current <- 8000
cost_new <- 15000

# Budget impact calculation
budget_impact <- data.frame(
  Year = 1:5,
  Uptake = uptake,
  Patients_New = pop_eligible * uptake,
  Patients_Current = pop_eligible * (1 - uptake),
  Cost_New = pop_eligible * uptake * cost_new,
  Cost_Current = pop_eligible * (1 - uptake) * cost_current,
  Cost_Reference = pop_eligible * cost_current
)

budget_impact <- budget_impact |>
  mutate(
    Total_Cost = Cost_New + Cost_Current,
    Budget_Impact = Total_Cost - Cost_Reference
  )

# Summary
print(budget_impact)

# Total 5-year budget impact
total_impact <- sum(budget_impact$Budget_Impact)
cat("Total 5-year budget impact:", round(total_impact / 1e6, 2), "million\n")

Budget Impact Visualization

library(ggplot2)
library(tidyr)

# Reshape for plotting
bi_long <- budget_impact |>
  select(Year, Cost_New, Cost_Current, Cost_Reference) |>
  pivot_longer(
    cols = -Year,
    names_to = "Category",
    values_to = "Cost"
  )

# Stacked bar chart
ggplot(bi_long |> filter(Category != "Cost_Reference"),
       aes(x = factor(Year), y = Cost / 1e6, fill = Category)) +
  geom_bar(stat = "identity") +
  geom_line(data = bi_long |> filter(Category == "Cost_Reference"),
            aes(y = Cost / 1e6, group = 1), linetype = "dashed", size = 1) +
  labs(x = "Year", y = "Cost ($ millions)",
       title = "Budget Impact Analysis",
       fill = "Treatment") +
  scale_fill_brewer(palette = "Set2") +
  theme_bw()

Decision Trees

Using dampack Package

library(dampack)

# Define decision tree parameters
params <- list(
  p_disease = 0.10,          # Probability of disease
  p_cure_trt = 0.80,         # Cure probability with treatment
  p_cure_notrt = 0.50,       # Cure probability without treatment
  c_test = 100,              # Cost of diagnostic test
  c_treatment = 5000,        # Cost of treatment
  c_disease = 20000,         # Cost of disease (if not cured)
  u_healthy = 1.0,           # Utility healthy
  u_disease = 0.60           # Utility with disease
)

# Strategy 1: Treat all
cost_treat_all <- params$c_treatment +
  params$p_disease * (1 - params$p_cure_trt) * params$c_disease
qaly_treat_all <- params$p_disease * (
  params$p_cure_trt * params$u_healthy +
  (1 - params$p_cure_trt) * params$u_disease
) + (1 - params$p_disease) * params$u_healthy

# Strategy 2: Test then treat
cost_test_treat <- params$c_test +
  params$p_disease * (params$c_treatment +
  (1 - params$p_cure_trt) * params$c_disease)

# Calculate all strategies and compare

Discounting

# Discount costs and effects
discount <- function(values, rate, time_points) {
  values / (1 + rate)^time_points
}

# Example: 30-year costs
years <- 0:29
annual_costs <- rep(5000, 30)
discount_rate <- 0.03

# Present value
pv_costs <- sum(discount(annual_costs, discount_rate, years))
cat("Undiscounted total:", sum(annual_costs), "\n")
cat("Present value (3% discount):", round(pv_costs, 0), "\n")

# Differential discounting (costs vs QALYs)
# Some guidelines recommend different rates
discount_costs <- 0.03
discount_qalys <- 0.015

pv_costs <- sum(discount(annual_costs, discount_costs, years))
pv_qalys <- sum(discount(annual_qalys, discount_qalys, years))

Key Packages Summary

Package	Purpose
BCEA	Bayesian cost-effectiveness analysis
heemod	Markov cohort models for HE
hesim	Health economic simulation modeling
dampack	Decision-analytic modeling tools
survHE	Survival analysis for HE
flexsurv	Parametric survival for extrapolation
CEAutil	CEA utility functions
valueEQ5D	EQ-5D utility mapping

Best Practices

Model structure: Match model type to disease natural history
Time horizon: Sufficient to capture all relevant costs and effects
Discounting: Apply recommended rates (often 3% for both costs and effects)
Uncertainty: Always conduct PSA and report CIs/CrIs
Transparency: Document all assumptions and data sources
Validation: Internal consistency checks and external validation
Reporting: Follow CHEERS checklist for publications
Perspective: Clearly state healthcare payer vs societal perspective