health-economics

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