# ============================================================
# GLP-1 Bayesian Generalized Bass Model
# Three-country scenario simulation: U.S., Canada, Mexico
#
# New scenario structure:
#
# Market potential m_t:
#   1. Current Clinical Coverage
#   2. Expanded Weight-Management Access
#
# Affordability and awareness:
#   1. Low Adoption
#   2. Mid Adoption
#   3. High Adoption
#
# Each country:
#   2 coverage scenarios x 3 adoption scenarios = 6 paths
#
# U.S.:
#   Estimate p, q, b_affordability, b_awareness using historical
#   cumulative ever-use data, then project 2026-2030.
#
# Canada and Mexico:
#   Prior-predictive simulations using country-specific p and q priors.
#   b_affordability and b_awareness are transferred from the U.S.
#   posterior median estimates.
# ============================================================

remove(list = ls())

# install.packages("R2jags")
# install.packages("coda")
# install.packages("dplyr")
# install.packages("tidyr")
# install.packages("ggplot2")

library(R2jags)
library(coda)
library(dplyr)
library(tidyr)
library(ggplot2)
library(purrr)

# ============================================================
# 0. Helper functions
# ============================================================

draw_truncated_normal <- function(n, mean, sd, lower, upper) {
  x <- rnorm(n, mean, sd)
  
  while (any(x < lower | x > upper)) {
    bad <- x < lower | x > upper
    x[bad] <- rnorm(sum(bad), mean, sd)
  }
  
  x
}

fitted_bayes_gbm_mt <- function(pars, Z, m_path) {
  
  p <- pars[1]
  q <- pars[2]
  b1 <- pars[3]
  b2 <- pars[4]
  
  Time <- nrow(Z) - 1
  unit <- 1:Time
  
  y_hat <- numeric(Time)
  
  for (i in 1:Time) {
    
    X_now <- unit[i] +
      b1 * Z[i + 1, 1] +
      b2 * Z[i + 1, 2]
    
    X_lag <- unit[i] - 1 +
      b1 * Z[i, 1] +
      b2 * Z[i, 2]
    
    F_now <- (1 - exp(-X_now * (p + q))) /
      (1 + (q / p) * exp(-X_now * (p + q)))
    
    F_lag <- (1 - exp(-X_lag * (p + q))) /
      (1 + (q / p) * exp(-X_lag * (p + q)))
    
    y_hat[i] <- m_path[i + 1] * F_now - m_path[i] * F_lag
  }
  
  y_hat[y_hat < 0] <- 0
  return(y_hat)
}

simulate_gbm_mt <- function(p, q, b1, b2, Z, m_path, years) {
  
  Time <- nrow(Z)
  flow <- numeric(Time)
  F_prev <- 0
  
  for (i in 1:Time) {
    
    X_now <- i +
      b1 * Z[i, "z_affordability"] +
      b2 * Z[i, "z_awareness"]
    
    F_now <- (1 - exp(-X_now * (p + q))) /
      (1 + (q / p) * exp(-X_now * (p + q)))
    
    if (i == 1) {
      flow[i] <- m_path[i] * F_now
    } else {
      flow[i] <- m_path[i] * F_now - m_path[i - 1] * F_prev
    }
    
    F_prev <- F_now
  }
  
  flow[flow < 0] <- 0
  cum <- cumsum(flow)
  
  data.frame(
    year = years,
    flow_adult_pct = flow * 100,
    cum_adult_pct = cum * 100
  )
}

# Manuscript equations:
#
# D_t = d_t A_{t-1}
# R_t = rho_t B_{t-1}
# A_t = A_{t-1}(1 - d_t) + S_t + rho_t B_{t-1}
# B_t = (1 - rho_t)B_{t-1} + lambda D_t
#
# where:
# A_t = current active adoption
# S_t = new first-time adoption from GBM
# D_t = annual discontinuers
# R_t = annual reinitiators
# B_t = returnable inactive prior-user stock
# lambda = share of discontinuers likely to reinitiate eventually
# rho_t = annual reinitiation rate among returnable inactive prior users

calculate_active_adoption <- function(df, persistence_table, lambda_default = 0.50) {
  
  df <- df %>%
    arrange(year) %>%
    left_join(persistence_table, by = "year") %>%
    mutate(
      new_adopter_rate_pct = flow_adult_pct
    )
  
  if (!"returnable_share" %in% names(df)) {
    df$returnable_share <- lambda_default
  }
  
  if (any(is.na(df$discontinuation_rate))) {
    stop("Missing discontinuation_rate after joining persistence_table.")
  }
  
  if (any(is.na(df$reinitiation_rate))) {
    stop("Missing reinitiation_rate after joining persistence_table.")
  }
  
  if (any(is.na(df$returnable_share))) {
    stop("Missing returnable_share after joining persistence_table.")
  }
  
  n <- nrow(df)
  
  active <- numeric(n)
  retained_previous_active <- numeric(n)
  gross_discontinued <- numeric(n)
  reinitiated_from_returnable_stock <- numeric(n)
  returnable_inactive_stock <- numeric(n)
  
  for (i in seq_len(n)) {
    
    if (i == 1) {
      
      # No prior active or inactive stock before first model year.
      gross_discontinued[i] <- 0
      retained_previous_active[i] <- 0
      reinitiated_from_returnable_stock[i] <- 0
      active[i] <- df$new_adopter_rate_pct[i]
      returnable_inactive_stock[i] <- 0
      
    } else {
      
      # D_t = d_t A_{t-1}
      gross_discontinued[i] <-
        df$discontinuation_rate[i] * active[i - 1]
      
      # A_{t-1}(1 - d_t)
      retained_previous_active[i] <-
        active[i - 1] * (1 - df$discontinuation_rate[i])
      
      # R_t = rho_t B_{t-1}
      reinitiated_from_returnable_stock[i] <-
        df$reinitiation_rate[i] * returnable_inactive_stock[i - 1]
      
      # A_t = A_{t-1}(1 - d_t) + S_t + R_t
      active[i] <-
        retained_previous_active[i] +
        df$new_adopter_rate_pct[i] +
        reinitiated_from_returnable_stock[i]
      
      # B_t = (1 - rho_t)B_{t-1} + lambda D_t
      returnable_inactive_stock[i] <-
        (1 - df$reinitiation_rate[i]) * returnable_inactive_stock[i - 1] +
        df$returnable_share[i] * gross_discontinued[i]
    }
  }
  
  df %>%
    mutate(
      retained_previous_active_pct = retained_previous_active,
      gross_discontinued_pct = gross_discontinued,
      reinitiated_from_returnable_stock_pct = reinitiated_from_returnable_stock,
      returnable_inactive_stock_pct = returnable_inactive_stock,
      current_active_adoption_pct = active,
      current_active_total_population_pct =
        current_active_adoption_pct * adult_share
    )
}

# ============================================================
# 1. U.S. historical estimation
# ============================================================

years_US <- 2014:2030

mt_table_US <- data.frame(
  year = years_US,
  
  adult_share = c(
    0.772, 0.773, 0.774, 0.775, 0.776, 0.777,
    0.779, 0.778, 0.780, 0.783, 0.785, 0.789,
    0.793, 0.795, 0.797, 0.798, 0.800
  ),
  
  adult_obesity_rate = c(
    0.379, 0.384, 0.390, 0.395, 0.400, 0.405,
    0.410, 0.415, 0.420, 0.425, 0.430, 0.435,
    0.440, 0.445, 0.450, 0.455, 0.460
  ),
  
  # Current Clinical Coverage:
  # This replaces the previous middle m_t path.
  m_current_clinical = c(
    0.045, 0.046, 0.047, 0.051, 0.056, 0.061,
    0.066, 0.080, 0.100, 0.125, 0.155, 0.190,
    0.230, 0.275, 0.330, 0.395, 0.460
  ),
  
  # Expanded Weight-Management Access:
  # This replaces the previous high m_t path.
  m_expanded_weight_management = c(
    0.045, 0.046, 0.047, 0.051, 0.056, 0.061,
    0.066, 0.080, 0.100, 0.125, 0.155, 0.190,
    0.240, 0.310, 0.405, 0.530, 0.680
  )
)

zt_raw_table_US <- data.frame(
  year = years_US,
  
  # Low Adoption: use current/original path
  low_affordability = c(
    0.03126, 0.03792, 0.04539, 0.05205, 0.06442,
    0.07598, 0.09244, 0.12046, 0.17460, 0.23540,
    0.26560, 0.30014, 0.32000, 0.36000, 0.40000,
    0.43000, 0.45000
  ),
  
  low_awareness = c(
    0.01, 0.01, 0.02, 0.03, 0.04,
    0.05, 0.06, 0.09, 0.13, 0.19,
    0.32, 0.46, 0.56, 0.64, 0.71,
    0.77, 0.82
  ),
  
  # Mid Adoption
  mid_affordability = c(
    0.03126, 0.03792, 0.04539, 0.05205, 0.06442,
    0.07598, 0.09244, 0.12046, 0.17460, 0.23540,
    0.26560, 0.30014, 0.35, 0.40, 0.45,
    0.500, 0.55
  ),
  
  mid_awareness = c(
    0.01, 0.01, 0.02, 0.03, 0.04,
    0.05, 0.06, 0.09, 0.13, 0.19,
    0.32, 0.46, 0.59, 0.77, 0.80,
    0.85, 0.87
  ),
  
  # High Adoption
  high_affordability = c(
    0.03126, 0.03792, 0.04539, 0.05205, 0.06442,
    0.07598, 0.09244, 0.12046, 0.17460, 0.23540,
    0.26560, 0.30014, 0.38, 0.46, 0.54,
    0.62, 0.70
  ),
  
  high_awareness = c(
    0.01, 0.01, 0.02, 0.03, 0.04,
    0.05, 0.06, 0.09, 0.13, 0.19,
    0.32, 0.46, 0.62, 0.79, 0.82,
    0.86, 0.92
  )
)

datos_US <- data.frame(
  year = 2014:2025,
  C_obs_pct = c(
    0.04, 0.12, 0.20, 0.30, 0.60, 1.00,
    1.50, 2.00, 4.00, 7.50, 12.00, 18.00
  )
)

datos_US <- datos_US %>%
  left_join(
    mt_table_US %>%
      select(
        year,
        adult_share,
        adult_obesity_rate,
        m_t = m_current_clinical
      ),
    by = "year"
  ) %>%
  left_join(
    zt_raw_table_US %>%
      select(
        year,
        affordability = low_affordability,
        awareness = low_awareness
      ),
    by = "year"
  ) %>%
  mutate(
    C_obs = C_obs_pct / 100,
    y_flow = C_obs - lag(C_obs, default = 0),
    y_flow_pct = y_flow * 100,
    F_cum = C_obs,
    F_cum_pct = C_obs_pct,
    unit = 1:n()
  )

Time_US <- nrow(datos_US)

affordability_base_US <- datos_US$affordability[1]
awareness_base_US <- datos_US$awareness[1]

Z_hist_year_US <- datos_US %>%
  mutate(
    z_affordability = log(affordability / affordability_base_US),
    z_awareness = log(awareness / awareness_base_US)
  ) %>%
  select(z_affordability, z_awareness)

Z_US <- rbind(
  c(0, 0),
  as.matrix(Z_hist_year_US)
)

colnames(Z_US) <- c("z_affordability", "z_awareness")

m_hist_US <- datos_US$m_t
m_jags_US <- c(m_hist_US[1], m_hist_US)

model_string <- "
model {

  for (i in 1:Time) {

    y[i] ~ dnorm(mu[i], tau)

    X_now[i] <- unit[i] +
      b[1] * z1[i + 1] +
      b[2] * z2[i + 1]

    X_lag[i] <- unit[i] - 1 +
      b[1] * z1[i] +
      b[2] * z2[i]

    F_now[i] <- (1 - exp(-X_now[i] * (p + q))) /
      (1 + (q / p) * exp(-X_now[i] * (p + q)))

    F_lag[i] <- (1 - exp(-X_lag[i] * (p + q))) /
      (1 + (q / p) * exp(-X_lag[i] * (p + q)))

    mu[i] <- m_t[i + 1] * F_now[i] - m_t[i] * F_lag[i]
  }


tempvar ~ dbeta(1,1)
sigma2 <- kappa*tempvar/(1-tempvar)
tau <- 1/sigma2
sigma <- sqrt(sigma2)


  p ~ dunif(0.003, 0.015)
  q ~ dunif(0.20, 0.55)

  b[1] ~ dnorm(0.25, 1 / pow(0.25, 2)) T(0, )
  b[2] ~ dnorm(0.25, 1 / pow(0.25, 2)) T(0, )
}
"

writeLines(model_string, "bass_gbm_model_positive_betas.txt")

kappa_US <- 5 * var(datos_US$y_flow)
n.chains <- 3
inits_gbm_US <- list(
  list(p=0.015,q=0.35,b=c(0.15,0.15),tempvar=0.30),
  list(p=0.025,q=0.50,b=c(0.25,0.25),tempvar=0.50),
  list(p=0.035,q=0.70,b=c(0.40,0.40),tempvar=0.70)
)

monitor <- c(
  "p",
  "q",
  "b",
  "tempvar",
  "tau",
  "sigma",
  "sigma2"
)

datos_list_US <- list(
  Time=Time_US,
  y=datos_US$y_flow,
  z1=Z_US[,1],
  z2=Z_US[,2],
  unit=1:Time_US,
  m_t=m_jags_US,
  kappa=kappa_US
)
n.iter <- 100000
n.burnin <- n.iter/2
n.thin <- 20
set.seed(123)
jagsfit_US <- jags(
  data=datos_list_US,
  inits=inits_gbm_US,
  parameters.to.save=monitor,
  model.file="bass_gbm_model_positive_betas.txt",
  n.chains=n.chains,
  n.iter=n.iter,
  n.burnin=n.burnin,
  n.thin=n.thin,
  progress.bar="gui"
)
print(jagsfit_US)


PARAMETROS_POSTERIORES_US <- data.frame(
  p = jagsfit_US$BUGSoutput$sims.list$p,
  q = jagsfit_US$BUGSoutput$sims.list$q,
  b_affordability = jagsfit_US$BUGSoutput$sims.list$b[, 1],
  b_awareness = jagsfit_US$BUGSoutput$sims.list$b[, 2]
)

print(summary(mcmc(PARAMETROS_POSTERIORES_US)))

coef_med_US <- c(
  p  = median(PARAMETROS_POSTERIORES_US$p),
  q  = median(PARAMETROS_POSTERIORES_US$q),
  b1 = median(PARAMETROS_POSTERIORES_US$b_affordability),
  b2 = median(PARAMETROS_POSTERIORES_US$b_awareness)
)

print(coef_med_US)


# ============================================================
#Adding 2 figures and 2 tables
# ============================================================
# ============================================================
# 1B. U.S. Bayesian estimation deliverables
#
# Deliverables:
#   Table 1: Posterior estimates and convergence diagnostics
#   Table 2: Historical U.S. adoption series and provenance
#   Figure 1: Observed versus fitted U.S. adoption
#   Figure 2: Prior versus posterior distributions
# ============================================================


# ------------------------------------------------------------
# Create output folder
# ------------------------------------------------------------

deliverable_dir_US <- "GLP1_GBM_estimation_deliverables_US"

dir.create(
  deliverable_dir_US,
  recursive = TRUE,
  showWarnings = FALSE
)


# ============================================================
# TABLE 1
# Posterior estimates and convergence diagnostics
# ============================================================

# R2jags already calculates:
# posterior mean,
# quantiles,
# R-hat,
# effective sample size
#
# These are stored in:
# jagsfit_US$BUGSoutput$summary


jags_summary_US <- jagsfit_US$BUGSoutput$summary


# Parameters to report.
#
# p                  = innovation parameter
# q                  = imitation parameter
# b[1]               = beta_AF, affordability effect
# b[2]               = beta_AW, awareness effect
# sigma              = residual standard deviation
# tempvar            = underlying sampled residual-variance parameter

parameters_to_report_US <- c(
  "p",
  "q",
  "b[1]",
  "b[2]",
  "sigma",
  "tempvar"
)


# Make sure all requested parameters are available

missing_parameters_US <- setdiff(
  parameters_to_report_US,
  rownames(jags_summary_US)
)

if (length(missing_parameters_US) > 0) {
  
  stop(
    paste(
      "These parameters are missing from the JAGS output:",
      paste(missing_parameters_US, collapse = ", "),
      "\nMake sure they are included in 'monitor' and rerun JAGS."
    )
  )
}


# Extract Table 1

Table1_US <- data.frame(
  
  parameter_raw = parameters_to_report_US,
  
  posterior_mean =
    jags_summary_US[parameters_to_report_US, "mean"],
  
  posterior_median =
    jags_summary_US[parameters_to_report_US, "50%"],
  
  percentile_2.5 =
    jags_summary_US[parameters_to_report_US, "2.5%"],
  
  percentile_97.5 =
    jags_summary_US[parameters_to_report_US, "97.5%"],
  
  R_hat =
    jags_summary_US[parameters_to_report_US, "Rhat"],
  
  effective_sample_size =
    jags_summary_US[parameters_to_report_US, "n.eff"],
  
  row.names = NULL
)


# Give parameters publication-friendly names

Table1_US$parameter <- c(
  "p",
  "q",
  "beta_AF (affordability)",
  "beta_AW (awareness)",
  "sigma (residual SD)",
  "tempvar"
)


# Reorder columns

Table1_US <- Table1_US %>%
  select(
    parameter,
    posterior_mean,
    posterior_median,
    percentile_2.5,
    percentile_97.5,
    R_hat,
    effective_sample_size
  )


# Print Table 1

print(Table1_US)


# Export Table 1

write.csv(
  Table1_US,
  file.path(
    deliverable_dir_US,
    "Table1_posterior_estimates_convergence_diagnostics.csv"
  ),
  row.names = FALSE
)


# ============================================================
# TABLE 2
# Historical U.S. adoption series and data provenance
# ============================================================
#
# IMPORTANT:
#
# The current R script contains the numerical adoption series,
# but it DOES NOT contain the actual historical source citations
# or explanations showing how each number was constructed.
#
# Therefore:
#
# 1. The model values are populated automatically below.
# 2. You MUST manually fill in:
#       original_source
#       source_reported_statistic
#       construction_flag
#       how_model_value_obtained
#
# Do NOT invent these.
#
# Suggested construction_flag values:
#
#   "Direct"
#   "Interpolated"
#   "Extrapolated"
#   "Calibrated"
#   "Combined sources"
#   "Converted"
#   "Other constructed"
#
# If more than one source contributes to one year,
# list all sources separated by semicolons.
# ============================================================


Table2_US <- data.frame(
  
  year = datos_US$year,
  
  # Historical cumulative percentage as shown in the paper/table
  cumulative_adoption_pct = datos_US$C_obs_pct,
  
  # Exact decimal quantity stored in C_obs
  model_C_obs_proportion = datos_US$C_obs,
  
  # The JAGS likelihood actually uses the annual increment
  annual_increment_pct = datos_US$y_flow_pct,
  
  annual_increment_proportion = datos_US$y_flow,
  
  # ---------------------------------------------
  # MANUALLY COMPLETE THE FOLLOWING FOUR COLUMNS
  # ---------------------------------------------
  
  original_source = rep(NA_character_, nrow(datos_US)),
  
  source_reported_statistic = rep(NA_character_, nrow(datos_US)),
  
  construction_flag = rep(NA_character_, nrow(datos_US)),
  
  how_model_value_obtained = rep(NA_character_, nrow(datos_US)),
  
  stringsAsFactors = FALSE
)


# ------------------------------------------------------------
# IMPORTANT:
# Replace the NA values above with the actual provenance.
#
# For example, once you know the sources, you could write:
#
# Table2_US$original_source[Table2_US$year == 2024] <-
#   "Author/organization, year, report/table/page"
#
# Table2_US$source_reported_statistic[Table2_US$year == 2024] <-
#   "Source reported X% of U.S. adults had ever used GLP-1 drugs"
#
# Table2_US$construction_flag[Table2_US$year == 2024] <-
#   "Direct"
#
# Table2_US$how_model_value_obtained[Table2_US$year == 2024] <-
#   "Used the reported cumulative ever-use percentage directly."
# ------------------------------------------------------------


# Check that Table 2 exactly reproduces the adoption series
# used in datos_US

stopifnot(
  isTRUE(
    all.equal(
      Table2_US$cumulative_adoption_pct,
      datos_US$C_obs_pct
    )
  )
)

stopifnot(
  isTRUE(
    all.equal(
      Table2_US$model_C_obs_proportion,
      datos_US$C_obs
    )
  )
)


# Check whether provenance information is still missing

missing_provenance_US <- Table2_US %>%
  filter(
    is.na(original_source) |
      is.na(source_reported_statistic) |
      is.na(construction_flag) |
      is.na(how_model_value_obtained)
  )


if (nrow(missing_provenance_US) > 0) {
  
  warning(
    paste0(
      "Table 2 still has missing provenance information for ",
      nrow(missing_provenance_US),
      " year(s). Complete these fields before using the table ",
      "in the manuscript."
    )
  )
}


# Print Table 2

print(Table2_US)


# Export Table 2

write.csv(
  Table2_US,
  file.path(
    deliverable_dir_US,
    "Table2_US_historical_adoption_and_provenance.csv"
  ),
  row.names = FALSE
)


# ============================================================
# FIGURE 1
# Observed versus fitted U.S. cumulative adoption
# ============================================================
#
# fitted_bayes_gbm_mt() returns annual adoption FLOWS.
#
# Because the observed series C_obs_pct is CUMULATIVE adoption,
# we need to:
#
# 1. calculate fitted annual flows for every posterior draw;
# 2. cumulatively sum the fitted flows;
# 3. calculate the posterior median and 95% credible interval.
# ============================================================


fitted_draws_US <- lapply(
  1:nrow(PARAMETROS_POSTERIORES_US),
  function(j) {
    
    pars_j <- c(
      PARAMETROS_POSTERIORES_US$p[j],
      PARAMETROS_POSTERIORES_US$q[j],
      PARAMETROS_POSTERIORES_US$b_affordability[j],
      PARAMETROS_POSTERIORES_US$b_awareness[j]
    )
    
    
    fitted_flow_j <- fitted_bayes_gbm_mt(
      pars = pars_j,
      Z = Z_US,
      m_path = m_jags_US
    )
    
    
    fitted_cumulative_j <- cumsum(fitted_flow_j) * 100
    
    
    data.frame(
      draw = j,
      year = datos_US$year,
      fitted_cumulative_pct = fitted_cumulative_j
    )
  }
) %>%
  bind_rows()


# Summarize fitted cumulative adoption by year

fitted_summary_US <- fitted_draws_US %>%
  group_by(year) %>%
  summarize(
    
    posterior_mean_fitted_pct =
      mean(fitted_cumulative_pct),
    
    posterior_median_fitted_pct =
      median(fitted_cumulative_pct),
    
    q2.5_fitted_pct =
      quantile(
        fitted_cumulative_pct,
        probs = 0.025,
        names = FALSE
      ),
    
    q97.5_fitted_pct =
      quantile(
        fitted_cumulative_pct,
        probs = 0.975,
        names = FALSE
      ),
    
    .groups = "drop"
  )


# Add observed values

Figure1_data_US <- fitted_summary_US %>%
  left_join(
    datos_US %>%
      select(
        year,
        observed_cumulative_pct = C_obs_pct
      ),
    by = "year"
  )


# Print underlying figure data

print(Figure1_data_US)


# Export figure data

write.csv(
  Figure1_data_US,
  file.path(
    deliverable_dir_US,
    "Figure1_observed_vs_fitted_data.csv"
  ),
  row.names = FALSE
)


# Create Figure 1

Figure1_US <- ggplot(
  Figure1_data_US,
  aes(x = year)
) +
  
  geom_ribbon(
    aes(
      ymin = q2.5_fitted_pct,
      ymax = q97.5_fitted_pct
    ),
    alpha = 0.20
  ) +
  
  geom_line(
    aes(
      y = posterior_median_fitted_pct,
      linetype = "Posterior median fitted"
    ),
    linewidth = 1.2
  ) +
  
  geom_point(
    aes(
      y = observed_cumulative_pct,
      shape = "Observed"
    ),
    size = 3
  ) +
  
  geom_line(
    aes(
      y = observed_cumulative_pct,
      linetype = "Observed"
    ),
    linewidth = 0.7
  ) +
  
  scale_x_continuous(
    breaks = datos_US$year
  ) +
  
  labs(
    title = "Observed and Bayesian GBM-Fitted U.S. GLP-1 Adoption",
    subtitle = "Posterior median fitted trajectory with 95% credible interval",
    x = "Year",
    y = "Cumulative Adult Adoption (%)",
    linetype = NULL,
    shape = NULL
  ) +
  
  theme_minimal(base_size = 15) +
  
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1
    )
  )


print(Figure1_US)


# Save Figure 1

ggsave(
  filename = file.path(
    deliverable_dir_US,
    "Figure1_observed_vs_fitted_US_adoption.png"
  ),
  plot = Figure1_US,
  width = 10,
  height = 6,
  dpi = 300
)


# ============================================================
# FIGURE 2
# Prior versus posterior distributions
#
# Parameters:
#   p
#   q
#   beta_AF = affordability
#   beta_AW = awareness
# ============================================================


set.seed(123)

n_prior_draws_US <- 50000


# ------------------------------------------------------------
# Generate draws from the SAME priors used in JAGS
# ------------------------------------------------------------

prior_draws_US <- data.frame(
  
  p = runif(
    n_prior_draws_US,
    min = 0.01,
    max = 0.04
  ),
  
  q = runif(
    n_prior_draws_US,
    min = 0.30,
    max = 0.80
  ),
  
  beta_AF = draw_truncated_normal(
    n = n_prior_draws_US,
    mean = 0.25,
    sd = 0.25,
    lower = 0,
    upper = Inf
  ),
  
  beta_AW = draw_truncated_normal(
    n = n_prior_draws_US,
    mean = 0.25,
    sd = 0.25,
    lower = 0,
    upper = Inf
  )
)


# Convert prior draws to long format

prior_long_US <- prior_draws_US %>%
  pivot_longer(
    cols = everything(),
    names_to = "parameter",
    values_to = "value"
  ) %>%
  mutate(
    distribution = "Prior"
  )


# ------------------------------------------------------------
# Convert posterior draws to the same format
# ------------------------------------------------------------

posterior_long_US <- PARAMETROS_POSTERIORES_US %>%
  transmute(
    p = p,
    q = q,
    beta_AF = b_affordability,
    beta_AW = b_awareness
  ) %>%
  pivot_longer(
    cols = everything(),
    names_to = "parameter",
    values_to = "value"
  ) %>%
  mutate(
    distribution = "Posterior"
  )


# Combine

prior_posterior_US <- bind_rows(
  prior_long_US,
  posterior_long_US
)


# Give parameter labels a fixed order

prior_posterior_US$parameter <- factor(
  prior_posterior_US$parameter,
  levels = c(
    "p",
    "q",
    "beta_AF",
    "beta_AW"
  ),
  labels = c(
    "p",
    "q",
    "\u03B2_AF: Affordability",
    "\u03B2_AW: Awareness"
  )
)


prior_posterior_US$distribution <- factor(
  prior_posterior_US$distribution,
  levels = c(
    "Prior",
    "Posterior"
  )
)


# ------------------------------------------------------------
# Create Figure 2
# ------------------------------------------------------------

Figure2_US <- ggplot(
  prior_posterior_US,
  aes(
    x = value,
    color = distribution,
    fill = distribution
  )
) +
  
  geom_density(
    alpha = 0.20,
    linewidth = 1.0
  ) +
  
  facet_wrap(
    ~ parameter,
    scales = "free",
    ncol = 2
  ) +
  
  labs(
    title = "Prior and Posterior Distributions of Bayesian GBM Parameters",
    x = "Parameter Value",
    y = "Density",
    color = NULL,
    fill = NULL
  ) +
  
  theme_minimal(base_size = 15) +
  
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold"),
    strip.text = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )


print(Figure2_US)


# Save Figure 2

ggsave(
  filename = file.path(
    deliverable_dir_US,
    "Figure2_prior_vs_posterior_US.png"
  ),
  plot = Figure2_US,
  width = 10,
  height = 7,
  dpi = 300
)


# ============================================================
# Final checks
# ============================================================

cat("\n")
cat("==============================================\n")
cat("U.S. Bayesian estimation deliverables created\n")
cat("==============================================\n")
cat("\n")

cat(
  "Output folder:",
  deliverable_dir_US,
  "\n\n"
)

cat(
  "Table 1 rows:",
  nrow(Table1_US),
  "\n"
)

cat(
  "Table 2 rows:",
  nrow(Table2_US),
  "\n"
)

cat(
  "Historical years:",
  min(Table2_US$year),
  "-",
  max(Table2_US$year),
  "\n"
)

cat(
  "Posterior draws used for Figure 1:",
  nrow(PARAMETROS_POSTERIORES_US),
  "\n"
)

cat("\n")

print(Table1_US)

cat("\n")


























# ============================================================
# 2. Generic country simulation function
# ============================================================

run_country_scenarios <- function(country_name,
                                  country_code,
                                  years_country,
                                  mt_table,
                                  zt_raw_table,
                                  param_draws,
                                  persistence_table,
                                  output_dir) {
  #Create output folder before saving plots.
  dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
  
  scenario_grid <- expand.grid(
    coverage_scenario = c(
      "Current Clinical Coverage",
      "Expanded Weight-Management Access"
    ),
    adoption_scenario = c(
      "Low Adoption",
      "Mid Adoption",
      "High Adoption"
    ),
    stringsAsFactors = FALSE
  ) %>%
    mutate(
      m_col = ifelse(
        coverage_scenario == "Current Clinical Coverage",
        "m_current_clinical",
        "m_expanded_weight_management"
      ),
      affordability_col = case_when(
        adoption_scenario == "Low Adoption" ~ "low_affordability",
        adoption_scenario == "Mid Adoption" ~ "mid_affordability",
        adoption_scenario == "High Adoption" ~ "high_affordability"
      ),
      awareness_col = case_when(
        adoption_scenario == "Low Adoption" ~ "low_awareness",
        adoption_scenario == "Mid Adoption" ~ "mid_awareness",
        adoption_scenario == "High Adoption" ~ "high_awareness"
      )
    )
  
  simulate_one_path <- function(coverage_scenario,
                                adoption_scenario,
                                m_col,
                                affordability_col,
                                awareness_col) {
    
    affordability_base <- zt_raw_table[[affordability_col]][1]
    awareness_base <- zt_raw_table[[awareness_col]][1]
    
    Z <- zt_raw_table %>%
      mutate(
        z_affordability = log(.data[[affordability_col]] / affordability_base),
        z_awareness = log(.data[[awareness_col]] / awareness_base)
      ) %>%
      select(z_affordability, z_awareness) %>%
      as.matrix()
    
    m_path <- mt_table[[m_col]]
    adult_share_path <- mt_table$adult_share
    
    lapply(1:nrow(param_draws), function(j) {
      
      tmp <- simulate_gbm_mt(
        p = param_draws$p[j],
        q = param_draws$q[j],
        b1 = param_draws$b_affordability[j],
        b2 = param_draws$b_awareness[j],
        Z = Z,
        m_path = m_path,
        years = years_country
      )
      
      tmp %>%
        mutate(
          draw = j,
          country = country_name,
          coverage_scenario = coverage_scenario,
          adoption_scenario = adoption_scenario,
          scenario = paste(coverage_scenario, adoption_scenario, sep = " | "),
          m_t_adult_pct = m_path * 100,
          adult_share = adult_share_path
        )
    }) %>%
      bind_rows()
  }
  
  posterior_forecast_all <- pmap_dfr(
    list(
      scenario_grid$coverage_scenario,
      scenario_grid$adoption_scenario,
      scenario_grid$m_col,
      scenario_grid$affordability_col,
      scenario_grid$awareness_col
    ),
    simulate_one_path
  )
  
  active_posterior_all <- posterior_forecast_all %>%
    group_by(country, coverage_scenario, adoption_scenario, scenario, draw) %>%
    group_modify(~ calculate_active_adoption(.x, persistence_table)) %>%
    ungroup()
  
  cumulative_adult_summary_2026_2030 <- posterior_forecast_all %>%
    filter(year >= 2026, year <= 2030) %>%
    group_by(country, coverage_scenario, adoption_scenario, scenario, year) %>%
    summarize(
      mean_cumulative_adult_pct = mean(cum_adult_pct),
      median_cumulative_adult_pct = median(cum_adult_pct),
      q2.5_cumulative_adult_pct = quantile(cum_adult_pct, 0.025),
      q97.5_cumulative_adult_pct = quantile(cum_adult_pct, 0.975),
      .groups = "drop"
    )
  
  new_adult_summary_2026_2030 <- active_posterior_all %>%
    filter(year >= 2026, year <= 2030) %>%
    group_by(country, coverage_scenario, adoption_scenario, scenario, year) %>%
    summarize(
      mean_new_adult_pct = mean(new_adopter_rate_pct),
      median_new_adult_pct = median(new_adopter_rate_pct),
      q2.5_new_adult_pct = quantile(new_adopter_rate_pct, 0.025),
      q97.5_new_adult_pct = quantile(new_adopter_rate_pct, 0.975),
      .groups = "drop"
    )
  
  current_active_adult_summary_2026_2030 <- active_posterior_all %>%
    filter(year >= 2026, year <= 2030) %>%
    group_by(country, coverage_scenario, adoption_scenario, scenario, year) %>%
    summarize(
      mean_current_active_adult_pct = mean(current_active_adoption_pct),
      median_current_active_adult_pct = median(current_active_adoption_pct),
      q2.5_current_active_adult_pct =
        quantile(current_active_adoption_pct, 0.025),
      q97.5_current_active_adult_pct =
        quantile(current_active_adoption_pct, 0.975),
      .groups = "drop"
    )
  
  current_active_total_population_summary_2026_2030 <- active_posterior_all %>%
    filter(year >= 2026, year <= 2030) %>%
    group_by(country, coverage_scenario, adoption_scenario, scenario, year) %>%
    summarize(
      mean_current_active_total_population_pct =
        mean(current_active_total_population_pct),
      median_current_active_total_population_pct =
        median(current_active_total_population_pct),
      q2.5_current_active_total_population_pct =
        quantile(current_active_total_population_pct, 0.025),
      q97.5_current_active_total_population_pct =
        quantile(current_active_total_population_pct, 0.975),
      .groups = "drop"
    )
  
  make_adoption_plot <- function(data, y_var, y_label, title_prefix) {
    
    data %>%
      mutate(
        coverage_scenario = factor(
          coverage_scenario,
          levels = c(
            "Current Clinical Coverage",
            "Expanded Weight-Management Access"
          )
        ),
        adoption_scenario = factor(
          adoption_scenario,
          levels = c("Low Adoption", "Mid Adoption", "High Adoption")
        )
      ) %>%
      ggplot(
        aes(
          x = year,
          y = .data[[y_var]],
          color = coverage_scenario
        )
      ) +
      geom_point(
        alpha = 0.04,
        size = 0.8,
        position = position_jitter(width = 0.08, height = 0)
      ) +
      stat_summary(
        fun = median,
        geom = "line",
        aes(group = coverage_scenario),
        linewidth = 1.1
      ) +
      stat_summary(
        fun = median,
        geom = "point",
        size = 2.2
      ) +
      facet_wrap(~ adoption_scenario) +
      labs(
        title = paste0(country_name, ": ", title_prefix),
        x = "Year",
        y = y_label,
        color = "Coverage Scenario"
      ) +
      theme_minimal(base_size = 15) +
      theme(
        legend.position = "bottom",
        plot.title = element_text(face = "bold"),
        panel.grid.minor = element_blank(),
        
        axis.title = element_text(size = 20),   # "Year"
        axis.text = element_text(size = 20),    # numbers
        strip.text = element_text(size = 20),   # Low/Mid/High Adoption
        legend.title = element_text(size = 20), # Coverage Scenario
        legend.text = element_text(size = 20)   # Current / Expanded
      )
  }
  
  plot_new_by_adoption <- make_adoption_plot(
    data = active_posterior_all,
    y_var = "new_adopter_rate_pct",
    y_label = "New Adult Adoption Rate, S(t) (%)",
    title_prefix = "Simulated GLP-1 New Adoption Among Adults"
  )
  
  plot_cum_by_adoption <- make_adoption_plot(
    data = posterior_forecast_all,
    y_var = "cum_adult_pct",
    y_label = "Cumulative Adult Adoption Rate, C(t) (%)",
    title_prefix = "Simulated GLP-1 Cumulative Adoption Among Adults"
  )
  
  plot_active_by_adoption <- make_adoption_plot(
    data = active_posterior_all,
    y_var = "current_active_adoption_pct",
    y_label = "Current Active Adult Adoption Rate, A(t) (%)",
    title_prefix = "Simulated GLP-1 Current Active Adoption Among Adults"
  )
  
  print(plot_new_by_adoption)
  print(plot_cum_by_adoption)
  print(plot_active_by_adoption)
  
  ggsave(
    filename = file.path(output_dir, paste0(country_code, "_new_adult_adoption_by_adoption_scenario.png")),
    plot = plot_new_by_adoption,
    width = 12,
    height = 6,
    dpi = 300
  )
  
  ggsave(
    filename = file.path(output_dir, paste0(country_code, "_cumulative_adult_adoption_by_adoption_scenario.png")),
    plot = plot_cum_by_adoption,
    width = 12,
    height = 6,
    dpi = 300
  )
  
  ggsave(
    filename = file.path(output_dir, paste0(country_code, "_current_active_adult_adoption_by_adoption_scenario.png")),
    plot = plot_active_by_adoption,
    width = 12,
    height = 6,
    dpi = 300
  )
  
  list(
    posterior_forecast_all = posterior_forecast_all,
    active_posterior_all = active_posterior_all,
    cumulative_adult_summary_2026_2030 = cumulative_adult_summary_2026_2030,
    new_adult_summary_2026_2030 = new_adult_summary_2026_2030,
    current_active_adult_summary_2026_2030 = current_active_adult_summary_2026_2030,
    current_active_total_population_summary_2026_2030 =
      current_active_total_population_summary_2026_2030,
    plot_new_by_adoption = plot_new_by_adoption,
    plot_cum_by_adoption = plot_cum_by_adoption,
    plot_active_by_adoption = plot_active_by_adoption
  )
}

# ============================================================
# 3. U.S. simulation inputs and run
# ============================================================


set.seed(123)

draw_id_US <- sample(
  1:nrow(PARAMETROS_POSTERIORES_US),
  min(5000, nrow(PARAMETROS_POSTERIORES_US)),
  replace = FALSE
)

param_draws_US <- PARAMETROS_POSTERIORES_US[draw_id_US, ]

persistence_table_US <- data.frame(
  year = years_US,
  
  discontinuation_rate = rep(0.648, length(years_US)),
  
  # rho_t:
  # Annual reinitiation rate among returnable inactive prior users.
  # This is calibrated from 36.3% one-year reinitiation among all discontinuers
  # and lambda = 0.50:
  # rho = 0.363 / 0.50 = 0.726, rounded to 0.73.
  reinitiation_rate = rep(0.73, length(years_US)),
  
  # lambda:
  # Share of discontinuers who enter the returnable inactive-user stock.
  returnable_share = rep(0.50, length(years_US))
)

results_US <- run_country_scenarios(
  country_name = "United States",
  country_code = "US",
  years_country = years_US,
  mt_table = mt_table_US,
  zt_raw_table = zt_raw_table_US,
  param_draws = param_draws_US,
  persistence_table = persistence_table_US,
  output_dir = "GLP1_GBM_outputs_US_2"
)





# ============================================================
# 4. Canada inputs and run
# ============================================================

years_CA <- 2015:2030

mt_table_CA <- data.frame(
  year = years_CA,
  
  adult_share = c(
    0.815, 0.816, 0.817, 0.818,
    0.819, 0.820, 0.821, 0.822,
    0.823, 0.824, 0.825, 0.826,
    0.827, 0.828, 0.829, 0.830
  ),
  
  # Current Clinical Coverage:
  # This replaces the previous middle m_t path.
  m_current_clinical = c(
    0.020, 0.020, 0.022, 0.023,
    0.025, 0.026, 0.029, 0.032,
    0.036, 0.041, 0.048, 0.058,
    0.071, 0.087, 0.107, 0.130
  ),
  
  # Expanded Weight-Management Access:
  # This replaces the previous high m_t path.
  m_expanded_weight_management = c(
    0.020, 0.020, 0.022, 0.023,
    0.025, 0.026, 0.029, 0.032,
    0.036, 0.041, 0.048, 0.065,
    0.085, 0.110, 0.140, 0.175
  )
)

zt_raw_table_CA <- data.frame(
  year = years_CA,
  
  # Low Adoption: use current/original path.
  low_affordability = c(
    0.020, 0.021, 0.022, 0.023,
    0.024, 0.025, 0.026, 0.027,
    0.029, 0.031, 0.033,
    0.040, 0.060, 0.085, 0.120, 0.150
  ),
  
  low_awareness = c(
    0.010, 0.011, 0.012, 0.013,
    0.014, 0.015, 0.016, 0.017,
    0.020, 0.050, 0.100,
    0.150, 0.200, 0.250, 0.300, 0.350
  ),
  
  # Mid Adoption.
  mid_affordability = c(
    0.020, 0.021, 0.022, 0.023,
    0.024, 0.025, 0.026, 0.027,
    0.029, 0.031, 0.033,
    0.050, 0.1, 0.150, 0.2, 0.250
  ),
  
  mid_awareness = c(
    0.010, 0.011, 0.012, 0.013,
    0.014, 0.015, 0.016, 0.017,
    0.020, 0.050, 0.100,
    0.200, 0.3, 0.4, 0.5, 0.6
  ),
  
  # High Adoption.
  high_affordability = c(
    0.020, 0.021, 0.022, 0.023,
    0.024, 0.025, 0.026, 0.027,
    0.029, 0.031, 0.033,
    0.080, 0.150, 0.220, 0.29, 0.36
  ),
  
  high_awareness = c(
    0.010, 0.011, 0.012, 0.013,
    0.014, 0.015, 0.016, 0.017,
    0.020, 0.050, 0.100,
    0.250, 0.4, 0.55, 0.7, 0.85
  )
)

set.seed(123)

n_draws_CA <- 5000

param_draws_CA <- data.frame(
  draw = 1:n_draws_CA,
  
  p = draw_truncated_normal(
    n = n_draws_CA,
    mean = 0.0045,
    sd = 0.004,
    lower = 0.005,
    upper = 0.020
  ),
  
  q = draw_truncated_normal(
    n = n_draws_CA,
    mean = 0.300,
    sd = 0.100,
    lower = 0.250,
    upper = 0.650
  ),
  
  b_affordability = rep(as.numeric(coef_med_US["b1"]), n_draws_CA),
  b_awareness = rep(as.numeric(coef_med_US["b2"]), n_draws_CA)
)

persistence_table_CA <- data.frame(
  year = years_CA,
  
  discontinuation_rate = rep(0.648, length(years_CA)),
  
  reinitiation_rate = rep(0.73, length(years_CA)),
  
  returnable_share = rep(0.50, length(years_CA))
)

results_CA <- run_country_scenarios(
  country_name = "Canada",
  country_code = "CA",
  years_country = years_CA,
  mt_table = mt_table_CA,
  zt_raw_table = zt_raw_table_CA,
  param_draws = param_draws_CA,
  persistence_table = persistence_table_CA,
  output_dir = "GLP1_GBM_outputs_CA_2"
)


# ============================================================
# 5. Mexico inputs and run
# ============================================================

years_MX <- 2015:2030

mt_table_MX <- data.frame(
  year = years_MX,
  
  adult_share = c(
    0.700, 0.703, 0.706, 0.709,
    0.712, 0.715, 0.718, 0.721,
    0.724, 0.727, 0.730,
    0.733, 0.736, 0.739, 0.742, 0.745
  ),
  
  # Current Clinical Coverage:
  # This replaces the previous middle m_t path.
  m_current_clinical = c(
    0.013, 0.013, 0.015, 0.017,
    0.020, 0.022, 0.025, 0.029,
    0.034, 0.042, 0.055, 0.065,
    0.078, 0.094, 0.113, 0.135
  ),
  
  # Expanded Weight-Management Access:
  # This replaces the previous high m_t path.
  m_expanded_weight_management = c(
    0.013, 0.013, 0.015, 0.017,
    0.020, 0.022, 0.025, 0.029,
    0.034, 0.042, 0.055, 0.072,
    0.092, 0.118, 0.150, 0.190
  )
)

zt_raw_table_MX <- data.frame(
  year = years_MX,
  
  # Low Adoption: use current/original path.
  low_affordability = c(
    0.012, 0.0125, 0.013, 0.0135,
    0.014, 0.0145, 0.015, 0.0155,
    0.0165, 0.018, 0.020,
    0.024, 0.031, 0.040, 0.052, 0.066
  ),
  
  low_awareness = c(
    0.006, 0.0065, 0.007, 0.0075,
    0.008, 0.009, 0.010, 0.011,
    0.014, 0.025, 0.050,
    0.075, 0.105, 0.140, 0.180, 0.225
  ),
  
  # Mid Adoption.
  mid_affordability = c(
    0.012, 0.0125, 0.013, 0.0135,
    0.014, 0.0145, 0.015, 0.0155,
    0.0165, 0.018, 0.020,
    0.030, 0.046, 0.065, 0.088, 0.115
  ),
  
  mid_awareness = c(
    0.006, 0.0065, 0.007, 0.0075,
    0.008, 0.009, 0.010, 0.011,
    0.014, 0.025, 0.050,
    0.095, 0.150, 0.215, 0.290, 0.375
  ),
  
  # High Adoption.
  high_affordability = c(
    0.012, 0.0125, 0.013, 0.0135,
    0.014, 0.0145, 0.015, 0.0155,
    0.0165, 0.018, 0.020,
    0.040, 0.065, 0.095, 0.130, 0.170
  ),
  
  high_awareness = c(
    0.006, 0.0065, 0.007, 0.0075,
    0.008, 0.009, 0.010, 0.011,
    0.014, 0.025, 0.050,
    0.120, 0.200, 0.290, 0.390, 0.500
  )
)

set.seed(123)

n_draws_MX <- 5000

param_draws_MX <- data.frame(
  draw = 1:n_draws_MX,
  
  p = draw_truncated_normal(
    n = n_draws_MX,
    mean = 0.003,
    sd = 0.0012,
    lower = 0.0005,
    upper = 0.008
  ),
  
  q = draw_truncated_normal(
    n = n_draws_MX,
    mean = 0.200,
    sd = 0.050,
    lower = 0.080,
    upper = 0.350
  ),
  
  b_affordability = rep(as.numeric(coef_med_US["b1"]), n_draws_MX),
  b_awareness = rep(as.numeric(coef_med_US["b2"]), n_draws_MX)
)

persistence_table_MX <- data.frame(
  year = years_MX,
  
  discontinuation_rate = rep(0.648, length(years_MX)),
  
  reinitiation_rate = rep(0.73, length(years_MX)),
  
  returnable_share = rep(0.50, length(years_MX))
)

results_MX <- run_country_scenarios(
  country_name = "Mexico",
  country_code = "MX",
  years_country = years_MX,
  mt_table = mt_table_MX,
  zt_raw_table = zt_raw_table_MX,
  param_draws = param_draws_MX,
  persistence_table = persistence_table_MX,
  output_dir = "GLP1_GBM_outputs_MX_3"
)



# ============================================================
# 6. Plot current active adult adoption rate for US, CA, and MX
#    Current vs Expanded scenarios, mid pathway only
# ============================================================

active_US_CA_MX_mid <- bind_rows(
  results_US$active_posterior_all,
  results_CA$active_posterior_all,
  results_MX$active_posterior_all
) %>%
  filter(adoption_scenario == "Mid Adoption") %>%
  mutate(
    country = case_when(
      country == "United States" ~ "U.S.",
      country == "Canada" ~ "Canada",
      country == "Mexico" ~ "Mexico",
      TRUE ~ country
    ),
    country = factor(country, levels = c("U.S.", "Canada", "Mexico")),
    coverage_scenario = factor(
      coverage_scenario,
      levels = c(
        "Current Clinical Coverage",
        "Expanded Weight-Management Access"
      )
    )
  )

active_US_CA_MX_mid_summary <- active_US_CA_MX_mid %>%
  group_by(country, coverage_scenario, year) %>%
  summarize(
    median_active_adult_pct = median(current_active_adoption_pct),
    q2.5_active_adult_pct = quantile(current_active_adoption_pct, 0.025),
    q97.5_active_adult_pct = quantile(current_active_adoption_pct, 0.975),
    .groups = "drop"
  )

plot_active_US_CA_MX_mid <- ggplot(
  active_US_CA_MX_mid_summary,
  aes(
    x = year,
    y = median_active_adult_pct,
    color = coverage_scenario,
    fill = coverage_scenario
  )
) +
  geom_ribbon(
    aes(
      ymin = q2.5_active_adult_pct,
      ymax = q97.5_active_adult_pct
    ),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  geom_point(size = 2.4) +
  facet_wrap(~ country, scales = "free_y") +
  labs(
    title = "Current Active GLP-1 Adoption Among Adults",
    x = "Year",
    y = "Current Active Adult Adoption Rate, A(t) (%)",
    color = "Coverage Scenario",
    fill = "Coverage Scenario"
  ) +
  theme_minimal(base_size = 15) +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", size = 22),
    axis.title = element_text(size = 22),
    axis.text = element_text(size = 22),
    strip.text = element_text(size = 22),
    legend.title = element_text(size = 22),
    legend.text = element_text(size = 22),
    panel.grid.minor = element_blank()
  )

print(plot_active_US_CA_MX_mid)

ggsave(
  filename = "GLP1_current_active_adoption_US_CA_MX.png",
  plot = plot_active_US_CA_MX_mid,
  width = 13,
  height = 6.5,
  dpi = 300
)



# ============================================================
# Export selected adoption summaries for US, CA, and MX
# ============================================================

output_dir_combined <- "GLP1_GBM_outputs_combined"

if (!dir.exists(output_dir_combined)) {
  dir.create(output_dir_combined, recursive = TRUE)
}

# Combine cumulative adult adoption summaries
combined_cumulative_adult_summary_2026_2030 <- bind_rows(
  results_US$cumulative_adult_summary_2026_2030,
  results_CA$cumulative_adult_summary_2026_2030,
  results_MX$cumulative_adult_summary_2026_2030
) %>%
  arrange(country, coverage_scenario, adoption_scenario, year)

# Combine new adult adoption summaries
combined_new_adult_summary_2026_2030 <- bind_rows(
  results_US$new_adult_summary_2026_2030,
  results_CA$new_adult_summary_2026_2030,
  results_MX$new_adult_summary_2026_2030
) %>%
  arrange(country, coverage_scenario, adoption_scenario, year)

# Combine current active adult adoption summaries
combined_current_active_adult_summary_2026_2030 <- bind_rows(
  results_US$current_active_adult_summary_2026_2030,
  results_CA$current_active_adult_summary_2026_2030,
  results_MX$current_active_adult_summary_2026_2030
) %>%
  arrange(country, coverage_scenario, adoption_scenario, year)

# Export only the three needed CSV files
write.csv(
  combined_cumulative_adult_summary_2026_2030,
  file.path(output_dir_combined, "combined_cumulative_adult_adoption_rate_2026_2030.csv"),
  row.names = FALSE
)

write.csv(
  combined_new_adult_summary_2026_2030,
  file.path(output_dir_combined, "combined_new_adult_adoption_rate_2026_2030.csv"),
  row.names = FALSE
)

write.csv(
  combined_current_active_adult_summary_2026_2030,
  file.path(output_dir_combined, "combined_current_active_adult_adoption_rate_2026_2030.csv"),
  row.names = FALSE
)