This analysis presents a Markov model for evaluating malaria intervention strategies in children under 5 years across urban and rural settings in sub-Saharan Africa. The model compares four intervention strategies across urban and rural settings in children under 5 years old.

Interventions Being Compared The analysis evaluates four intervention strategies (all layered on top of ITNs as baseline prevention):

Baseline (ITNs only): Universal coverage with insecticide-treated nets (current standard of care)

Urban: 70% coverage

Rural: 45% coverage

SMC + ITNs: Seasonal malaria chemoprevention (4 cycles/year) plus ITNs

Target: Children 2-5 years (primary SMC candidates)

Efficacy: 85% against clinical malaria

Cost: $6.00 per child per year

Vaccine + ITNs: RTS,S/AS01 malaria vaccine plus ITNs

Target: Children 0-24 months (age-appropriate for vaccination)

Efficacy: 30% against clinical infection, 70% against severe disease

Cost: $26.00 for full 4-dose series

Combined Strategy: SMC + Vaccine + ITNs (maximal intervention)

Layered protection across all age groups

Accounts for geographic delivery costs (rural multiplier: 1.5×)

Model Overview Time horizon: 10 years

Cohort size: 10,000 children (scaled from birth to age 5+)

Discount rate: 3% per annum for both costs and health outcomes (DALYs)

Perspective: Healthcare system (direct medical costs + intervention delivery)

Outcomes: Clinical cases, severe cases, deaths, DALYs (discounted), costs (discounted)

# ===================================================================
# STEP 1: ENVIRONMENT SETUP AND PACKAGE LOADING
# ===================================================================

# Clear environment to prevent conflicts from previous runs
rm(list = ls())

# Force garbage collection to free memory
gc()
##           used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells  618606 33.1    1417320 75.7         NA   700251 37.4
## Vcells 1177539  9.0    8388608 64.0      49152  1963400 15.0
# Load required packages for analysis and visualization
library(ggplot2)   # For creating publication-quality plots
library(dplyr)     # For data manipulation and aggregation
library(tidyr)     # For reshaping data between wide/long formats

# Set random seed for reproducibility of any stochastic elements
set.seed(2026)

# Helper function for formatted numbers with commas
#  Use vectorized operations with ifelse() 
format_number <- function(x) {
  # Handle vectors properly
  result <- ifelse(is.na(x) | is.infinite(x), "NA", 
                   format(round(x, 0), big.mark = ",", scientific = FALSE))
  return(result)
}

Model Structure

# ===================================================================
# MARKOV MODEL STRUCTURE DIAGRAM 
# ===================================================================
# This code chunk visualizes the Markov model structure showing 
# health states and transitions used in the simulation

# Load required packages
if(!require(diagram)) {
  install.packages("diagram")
  library(diagram)
}

# Set up plotting parameters with larger figure
par(mar = c(2, 2, 3, 2), bg = "white")

# Define the model structure with states and transitions
states <- c(
  "Susceptible\n(Uninfected)",
  "Clinical Malaria\n(Non-severe)",
  "Severe Malaria",
  "Death"
)

# Define colors for different state types
state_colors <- c(
  "Susceptible" = "#4DAF4A",      # Green for healthy state
  "Clinical" = "#FFD92F",          # Yellow for mild disease
  "Severe" = "#F46D43",            # Orange-red for severe disease
  "Death" = "#377EB8"              # Blue for absorbing state
)

# Define transition colors
transition_colors <- c(
  "progression" = "#E41A1C",       # Red for disease progression
  "mortality" = "#000000",         # Black for death
  "recovery" = "#4DAF4A"           # Green for recovery
)

# Create the plot with larger dimensions
openplotmat(xlim = c(-0.2, 1.2), ylim = c(-0.1, 1.1))

# Position states with increased spacing between Severe and Death
pos <- matrix(nrow = 4, ncol = 2)
pos[1, ] <- c(0.5, 0.85)   # Susceptible (top)
pos[2, ] <- c(0.5, 0.55)   # Clinical (upper middle)
pos[3, ] <- c(0.25, 0.25)  # Severe (bottom left - moved further left)
pos[4, ] <- c(0.75, 0.25)  # Death (bottom right - moved further right)

# ============================================================
# DRAW TRANSITION ARROWS USING straightarrow()
# ============================================================

# 1. Susceptible → Clinical (disease progression)
straightarrow(from = pos[1, ], to = pos[2, ], lwd = 2.5, 
              arr.pos = 0.6, arr.type = "triangle",
              lcol = transition_colors["progression"], 
              arr.col = transition_colors["progression"])

# Add label for this transition
textrect(mid = c(0.55, 0.70), radx = 0.10, rady = 0.025, 
         lab = expression(bold(paste("Infection ×", ~delta[clinical]))), 
         cex = 0.9, shadow.size = 0, box.col = "white", 
         col = transition_colors["progression"])

# 2. Clinical → Severe (progression to severe disease)
straightarrow(from = pos[2, ], to = pos[3, ], lwd = 2.5,
              arr.pos = 0.6, arr.type = "triangle",
              lcol = transition_colors["progression"],
              arr.col = transition_colors["progression"])

# Adding label for this transition 
textrect(mid = c(0.35, 0.42), radx = 0.08, rady = 0.025, 
         lab = expression(bold(paste(delta[severe]))), 
         cex = 0.9, shadow.size = 0, box.col = "white",
         col = transition_colors["progression"])

# 3. Severe → Death (mortality)
straightarrow(from = pos[3, ], to = pos[4, ], lwd = 2.5,
              arr.pos = 0.6, arr.type = "triangle",
              lcol = transition_colors["mortality"],
              arr.col = transition_colors["mortality"])

# Add label for this transition - centered between boxes
textrect(mid = c(0.5, 0.19), radx = 0.10, rady = 0.025, 
         lab = expression(bold(paste("Case fatality =", ~CFR))), 
         cex = 0.9, shadow.size = 0, box.col = "white",
         col = transition_colors["mortality"])

# 4. Clinical → Susceptible (recovery with immunity) - LEFT ARC
curvedarrow(from = pos[2, ], to = pos[1, ], lwd = 2.5,
            arr.pos = 0.6, curve = -0.4, 
            lcol = transition_colors["recovery"], 
            arr.col = transition_colors["recovery"],
            arr.type = "triangle")

# Add label for recovery with immunity - placed further left
textrect(mid = c(0.30, 0.70), radx = 0.12, rady = 0.025, 
         lab = "Recovery +\nImmunity", cex = 0.85, 
         shadow.size = 0, box.col = "#A8E6A3", 
         col = transition_colors["recovery"],
         font = 2)

# 5. Severe → Susceptible (recovery from severe) - LEFT ARC
curvedarrow(from = pos[3, ], to = pos[1, ], lwd = 2.5,
            arr.pos = 0.6, curve = -0.5,
            lcol = transition_colors["recovery"],
            arr.col = transition_colors["recovery"],
            arr.type = "triangle")

# Add label for recovery from severe - placed further left
textrect(mid = c(0.10, 0.50), radx = 0.10, rady = 0.025, 
         lab = "Recovery", cex = 0.85, 
         shadow.size = 0, box.col = "#A8E6A3",
         col = transition_colors["recovery"],
         font = 2)

# ============================================================
# DRAW STATE BOXES
# ============================================================

# 1. Susceptible state (top)
textrect(mid = pos[1, ], radx = 0.15, rady = 0.07, 
         lab = states[1], cex = 1.0, shadow.size = 0.01,
         box.col = state_colors["Susceptible"], col = "white", 
         lwd = 3, font = 2)

# 2. Clinical state (middle)
textrect(mid = pos[2, ], radx = 0.15, rady = 0.07, 
         lab = states[2], cex = 1.0, shadow.size = 0.01,
         box.col = state_colors["Clinical"], col = "#333333", 
         lwd = 3, font = 2)

# 3. Severe state (bottom left)
textrect(mid = pos[3, ], radx = 0.13, rady = 0.07, 
         lab = states[3], cex = 1.0, shadow.size = 0.01,
         box.col = state_colors["Severe"], col = "white", 
         lwd = 3, font = 2)

# 4. Death state (bottom right)
textrect(mid = pos[4, ], radx = 0.12, rady = 0.07, 
         lab = states[4], cex = 1.0, shadow.size = 0.01,
         box.col = state_colors["Death"], col = "white", 
         lwd = 3, font = 2)

# ============================================================
# ADD TITLE AND ANNOTATIONS
# ============================================================

# Main title
title(main = "Markov Model Structure for Malaria Progression\n", 
      cex.main = 1.4, line = 0.8, font.main = 2)

# Subtitle
text(x = 0.5, y = 0.96, 
     labels = "Monthly transitions with time-dependent discounting (3% annually)",
     cex = 0.85, font = 3, col = "#555555")

# ============================================================
# ADD LEGEND AT TOP RIGHT
# ============================================================

# Legend for transition types positioned at top right
legend(x = 1.2, y = 0.75, 
       legend = c("Disease progression", "Mortality", "Recovery"),
       col = c(transition_colors["progression"], 
               transition_colors["mortality"], 
               transition_colors["recovery"]), 
       lwd = 3, cex = 0.85, bty = "o", bg = "white",
       box.col = "black", title = "Transition Types",
       horiz = FALSE, xjust = 1, yjust = 1)

# ============================================================
# ADD DISCOUNTING NOTE AT BOTTOM LEFT
# ============================================================

text(x = 0.02, y = 0.08, 
     labels = "NOTE: Only costs and DALYs are discounted (3% annually)",
     cex = 0.75, col = "#E41A1C", adj = 0, font = 2)

text(x = 0.02, y = 0.05, 
     labels = "Clinical cases, severe cases, and deaths are NOT discounted",
     cex = 0.70, col = "#E41A1C", adj = 0, font = 2)

# ============================================================
# ADD PARAMETER DEFINITIONS AT BOTTOM
# ============================================================

text(x = 0.02, y = 0.02, 
     labels = expression(paste(delta[clinical], " = Probability clinical malaria | infection")),
     cex = 0.65, col = "#555555", adj = 0)

text(x = 0.35, y = 0.02, 
     labels = expression(paste(delta[severe], " = Probability severe | clinical")),
     cex = 0.65, col = "#555555", adj = 0)

text(x = 0.68, y = 0.02, 
     labels = "CFR = Case Fatality Rate",
     cex = 0.65, col = "#555555", adj = 0)

# ============================================================
# ADD AGE AND GEOGRAPHY NOTE AT TOP RIGHT BELOW LEGEND
# ============================================================

text(x = 0.98, y = 0.88, 
     labels = "Age-specific risks\nGeography-specific modifiers",
     cex = 0.65, col = "#377EB8", adj = 1, font = 2)

Part 1: Model Parameters

Geographic and Demographic Settings- To be updated

# ===================================================================
# STEP 2: DEFINE GEOGRAPHIC AND DEMOGRAPHIC PARAMETERS
# ===================================================================

# Geographic strata reflecting transmission heterogeneity
geo_settings <- c("Urban", "Rural")
geo_pop_weights <- c(0.30, 0.70)  # 30% urban, 70% rural population distribution

# Geographic-specific parameters that influence transmission and outcomes
geo_params <- list(
  Urban = list(
    eir_multiplier = 0.4,        # Lower transmission intensity in urban settings
    itn_coverage = 0.70,          # Higher ITN access in urban areas
    vaccine_coverage = 0.65,      # Moderate vaccine uptake
    smc_coverage = 0.55,          # Moderate SMC coverage
    cfr_multiplier = 1.0,         # Baseline case fatality (reference)
    access_delay_days = 0.5       # Rapid healthcare access (12 hours)
  ),
  Rural = list(
    eir_multiplier = 2.0,         # Higher transmission (5× urban baseline)
    itn_coverage = 0.45,          # Lower ITN coverage due to access barriers
    vaccine_coverage = 0.40,      # Lower vaccine coverage
    smc_coverage = 0.45,          # Lower SMC coverage
    cfr_multiplier = 1.6,         # 60% higher case fatality from delayed care
    access_delay_days = 4.5       # Delayed healthcare access (4.5 days)
  )
)

# ===================================================================
# STEP 3: AGE GROUP DEFINITIONS
# ===================================================================

# Age groups reflect key epidemiological periods in malaria immunity development
age_groups <- c("0-6m", "6-24m", "2-5y", "5+y")
# Population weights based on typical demographic profiles in high-burden countries
age_weights <- c(0.05, 0.11, 0.16, 0.68)

# ===================================================================
# STEP 4: EPIDEMIOLOGICAL AND INTERVENTION PARAMETERS
# ===================================================================

base_eir <- 25                    # Baseline EIR: 25 infectious bites/person/year
smc_efficacy <- 0.85              # SMC reduces clinical malaria by 85%
vaccine_ve_infection <- 0.30      # RTS,S efficacy against clinical infection
vaccine_ve_severe <- 0.70         # RTS,S efficacy against severe disease
itn_efficacy <- 0.45              # ITNs reduce transmission by 45%

# ===================================================================
# STEP 5: COST PARAMETERS (2024 USD)
# ===================================================================

# Intervention costs
cost_vaccine_full <- 26.00        # Full 4-dose series per child
cost_smc_per_child_per_year <- 6.00  # 4 cycles at $1.50 per cycle

# Treatment costs (direct medical)
cost_act <- 1.50                  # Artemisinin combination therapy course
cost_outpatient <- 8.00           # Outpatient visit (clinical case)
cost_hospital <- 150.00           # Hospitalization (severe case)

# ===================================================================
# STEP 6: DISABILITY-ADJUSTED LIFE YEAR (DALY) PARAMETERS
# ===================================================================

# Based on Global Burden of Disease (GBD) 2021 estimates
dw_clinical <- 0.006              # Disability weight for uncomplicated malaria
dw_severe <- 0.133                # Disability weight for severe malaria
duration_clinical_years <- 14/365  # 14 days duration (converted to years)
duration_severe_years <- 21/365    # 21 days duration (converted to years)

# Life expectancy for YLL calculations (standard for sub-Saharan Africa)
life_expectancy <- 65

# ===================================================================
# STEP 7: DISCOUNTING PARAMETERS
# ===================================================================

# Standard 3% annual discount rate for costs and DALYs only
discount_rate <- 0.03  # 3% per annum

# Display confirmation of loaded parameters
cat("✓ Parameters loaded successfully\n")
## ✓ Parameters loaded successfully
cat("✓ Discount rate: 3% per annum (applied to COSTS and DALYs ONLY)\n")
## ✓ Discount rate: 3% per annum (applied to COSTS and DALYs ONLY)
cat("✓ Clinical cases, severe cases, and deaths are NOT discounted\n")
## ✓ Clinical cases, severe cases, and deaths are NOT discounted
cat("✓ Population distribution:\n")
## ✓ Population distribution:
print(data.frame(Geography = geo_settings, Population_Weight = geo_pop_weights))
##   Geography Population_Weight
## 1     Urban               0.3
## 2     Rural               0.7

Part 2: Markov Model Implementation

# ===================================================================
# STEP 8: MAIN MODEL FUNCTION - RUNS MARKOV SIMULATION
# ===================================================================
# NOTE: Only costs and DALYs are discounted at 3% per annum
# Clinical cases, severe cases, and deaths are reported as undiscounted counts

run_stable_model <- function(intervention = "baseline", 
                             time_horizon_years = 10,
                             cohort_size = 10000) {
  
  # Initialize accumulators for population-level outcomes (UNDISCOUNTED)
  total_clinical <- 0       
  total_severe <- 0         
  total_deaths <- 0         
  
  # Discounted cost and DALY accumulators
  total_treatment_cost_discounted <- 0 
  discounted_dalys <- 0
  
  # Data frame to store detailed results
  results <- data.frame()
  
  # Track intervention-specific costs
  intervention_cost_undiscounted <- 0
  intervention_cost_discounted <- 0
  
  # Loop over geographic settings
  for(geo in c("Urban", "Rural")) {
    geo_pop <- cohort_size * geo_pop_weights[ifelse(geo == "Urban", 1, 2)]
    gp <- geo_params[[geo]]
    
    # Calculate effective transmission intensity
    effective_eir <- base_eir * gp$eir_multiplier
    monthly_infection_prob <- 1 - exp(-effective_eir / 12 * 0.5)
    itn_effect <- 1 - (itn_efficacy * gp$itn_coverage)
    monthly_infection_prob <- monthly_infection_prob * itn_effect
    
    # Apply intervention-specific effects
    delivery_multiplier <- 1
    if(intervention == "smc") {
      smc_effect <- 1 - (smc_efficacy * gp$smc_coverage)
      monthly_infection_prob <- monthly_infection_prob * smc_effect
      intervention_cost_undiscounted <- intervention_cost_undiscounted + 
        geo_pop * cost_smc_per_child_per_year * time_horizon_years
      
    } else if(intervention == "vaccine") {
      vaccine_effect <- 1 - (vaccine_ve_infection * gp$vaccine_coverage)
      monthly_infection_prob <- monthly_infection_prob * vaccine_effect
      vaccine_cost <- geo_pop * cost_vaccine_full * gp$vaccine_coverage
      intervention_cost_undiscounted <- intervention_cost_undiscounted + vaccine_cost
      # Vaccine at year 0, no discounting
      intervention_cost_discounted <- intervention_cost_discounted + vaccine_cost
      
    } else if(intervention == "combined") {
      smc_effect <- 1 - (smc_efficacy * gp$smc_coverage)
      vaccine_effect <- 1 - (vaccine_ve_infection * gp$vaccine_coverage)
      monthly_infection_prob <- monthly_infection_prob * smc_effect * vaccine_effect
      
      smc_cost <- geo_pop * cost_smc_per_child_per_year * time_horizon_years
      vaccine_cost <- geo_pop * cost_vaccine_full * gp$vaccine_coverage
      delivery_multiplier <- ifelse(geo == "Rural", 1.5, 1)
      intervention_cost_undiscounted <- intervention_cost_undiscounted + 
        (smc_cost + vaccine_cost) * delivery_multiplier
      # Vaccine at year 0 (no discount), SMC will be discounted monthly
      intervention_cost_discounted <- intervention_cost_discounted + 
        vaccine_cost * delivery_multiplier
    }
    
    # Loop over age groups
    for(a in 1:length(age_groups)) {
      age <- age_groups[a]
      age_pop <- geo_pop * age_weights[a]
      
      # Age-specific clinical progression parameters
      if(age == "0-6m") {
        clinical_risk <- 0.25
        severe_risk <- 0.08
        cfr <- 0.18 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
      } else if(age == "6-24m") {
        clinical_risk <- 0.20
        severe_risk <- 0.06
        cfr <- 0.12 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
      } else if(age == "2-5y") {
        clinical_risk <- 0.10
        severe_risk <- 0.02
        cfr <- 0.06 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
      } else {
        clinical_risk <- 0.04
        severe_risk <- 0.005
        cfr <- 0.03 * gp$cfr_multiplier * (1 + gp$access_delay_days/10)
      }
      
      # Apply vaccine effect on severe disease
      if(intervention %in% c("vaccine", "combined")) {
        severe_risk <- severe_risk * (1 - vaccine_ve_severe * gp$vaccine_coverage)
      }
      
      # Initialize susceptible population
      susceptible <- age_pop
      clinical_cases <- 0      # UNDISCOUNTED count
      severe_cases <- 0        # UNDISCOUNTED count
      deaths <- 0              # UNDISCOUNTED count
      
      # Discounted accumulators for costs and DALYs
      discounted_treatment_cost_group <- 0
      discounted_smc_cost_group <- 0
      discounted_yld <- 0
      discounted_yll <- 0
      
      # Monthly simulation
      for(month in 1:(time_horizon_years * 12)) {
        # Calculate discount factor for costs and DALYs (3% annual)
        discount_factor <- 1 / (1 + discount_rate)^(month/12)
        
        # Disease transmission (UNDISCOUNTED - counts are NOT discounted)
        new_infections <- susceptible * monthly_infection_prob
        new_clinical <- new_infections * clinical_risk
        new_severe <- new_clinical * severe_risk
        new_deaths <- new_severe * cfr
        
        # Accumulate UNDISCOUNTED counts
        clinical_cases <- clinical_cases + new_clinical
        severe_cases <- severe_cases + new_severe
        deaths <- deaths + new_deaths
        
        # ============================================================
        # DISCOUNTED COST CALCULATIONS (only costs are discounted)
        # ============================================================
        transport_cost <- ifelse(geo == "Urban", 2, 8)
        monthly_treatment_cost <- (new_clinical * (cost_outpatient + cost_act + transport_cost)) +
                                  (new_severe * (cost_hospital + cost_act + transport_cost))
        discounted_treatment_cost_group <- discounted_treatment_cost_group + 
                                          (monthly_treatment_cost * discount_factor)
        
        # ============================================================
        # DISCOUNTED SMC COSTS (for combined and SMC-only strategies)
        # ============================================================
        if(intervention %in% c("smc", "combined")) {
          if(intervention == "smc") {
            monthly_smc_cost <- geo_pop * (cost_smc_per_child_per_year / 12)
          } else {
            monthly_smc_cost <- geo_pop * (cost_smc_per_child_per_year / 12) * delivery_multiplier
          }
          discounted_smc_cost_group <- discounted_smc_cost_group + 
                                       (monthly_smc_cost * discount_factor)
        }
        
        # ============================================================
        # DISCOUNTED DALY CALCULATIONS (DALYs are discounted)
        # ============================================================
        # YLD (Years Lived with Disability) - discounted
        discounted_yld <- discounted_yld + 
          (new_clinical * dw_clinical * duration_clinical_years * discount_factor) +
          (new_severe * dw_severe * duration_severe_years * discount_factor)
        
        # YLL (Years of Life Lost) - discounted
        if(age == "0-6m") {
          yll_per_death <- life_expectancy - 0.25
        } else if(age == "6-24m") {
          yll_per_death <- life_expectancy - 1.25
        } else if(age == "2-5y") {
          yll_per_death <- life_expectancy - 3.5
        } else {
          yll_per_death <- life_expectancy - 10
        }
        discounted_yll <- discounted_yll + (new_deaths * yll_per_death * discount_factor)
        
        # Update susceptible population (UNDISCOUNTED)
        susceptible <- susceptible - new_infections
        
        # Numerical stability - maintain minimum susceptible population
        if(susceptible < age_pop * 0.1) {
          susceptible <- age_pop * 0.1
        }
      }
      
      # Total discounted DALYs for this group
      total_dalys_discounted <- discounted_yll + discounted_yld
      
      # Aggregate discounted costs
      if(intervention %in% c("smc", "combined")) {
        intervention_cost_discounted <- intervention_cost_discounted + discounted_smc_cost_group
      }
      total_treatment_cost_discounted <- total_treatment_cost_discounted + 
                                        discounted_treatment_cost_group
      discounted_dalys <- discounted_dalys + total_dalys_discounted
      
      # Store results (cases/deaths are UNDISCOUNTED, costs/DALYs are DISCOUNTED)
      results <- rbind(results, data.frame(
        geography = geo,
        age_group = age,
        clinical_cases = clinical_cases,           # UNDISCOUNTED
        severe_cases = severe_cases,               # UNDISCOUNTED
        deaths = deaths,                           # UNDISCOUNTED
        dalys_discounted = total_dalys_discounted, # DISCOUNTED
        treatment_cost_discounted = discounted_treatment_cost_group  # DISCOUNTED
      ))
      
      # Update undiscounted totals
      total_clinical <- total_clinical + clinical_cases
      total_severe <- total_severe + severe_cases
      total_deaths <- total_deaths + deaths
    }
  }
  
  total_costs_discounted <- total_treatment_cost_discounted + intervention_cost_discounted
  
  return(list(
    results = results,
    total_clinical = total_clinical,                      # UNDISCOUNTED
    total_severe = total_severe,                          # UNDISCOUNTED
    total_deaths = total_deaths,                          # UNDISCOUNTED
    total_dalys_discounted = discounted_dalys,            # DISCOUNTED
    total_costs_discounted = total_costs_discounted,      # DISCOUNTED
    intervention_cost_discounted = intervention_cost_discounted,
    treatment_cost_discounted = total_treatment_cost_discounted,
    intervention_cost_undiscounted = intervention_cost_undiscounted
  ))
}

Part 3: Run All Strategies

# ===================================================================
# STEP 9: EXECUTE MODEL FOR ALL FOUR INTERVENTION STRATEGIES
# ===================================================================

strategies <- c("baseline", "smc", "vaccine", "combined")
strategy_names <- c("Baseline (ITNs)", "SMC + ITNs", "Vaccine + ITNs", "Combined")

results_list <- list()

cat("\n========================================\n")
## 
## ========================================
cat("RUNNING MALARIA MODEL SIMULATIONS\n")
## RUNNING MALARIA MODEL SIMULATIONS
cat("========================================\n")
## ========================================
cat("NOTE: Cases and deaths are UNDISCOUNTED (raw counts)\n")
## NOTE: Cases and deaths are UNDISCOUNTED (raw counts)
cat("      Costs and DALYs are DISCOUNTED at 3% per annum\n")
##       Costs and DALYs are DISCOUNTED at 3% per annum
cat("========================================\n\n")
## ========================================
for(i in 1:length(strategies)) {
  cat(sprintf("▶ Executing %s...\n", strategy_names[i]))
  
  result <- run_stable_model(intervention = strategies[i], 
                             time_horizon_years = 10,
                             cohort_size = 10000)
  
  results_list[[i]] <- result
  
  cat(sprintf("  ✓ Clinical cases (undiscounted): %s\n", format_number(result$total_clinical)))
  cat(sprintf("  ✓ Severe cases (undiscounted): %s\n", format_number(result$total_severe)))
  cat(sprintf("  ✓ Deaths (undiscounted): %s\n", format_number(result$total_deaths)))
  cat(sprintf("  ✓ DALYs (discounted at 3%%): %s\n", format_number(result$total_dalys_discounted)))
  cat(sprintf("  ✓ Total costs (discounted at 3%%): $%s\n", format_number(result$total_costs_discounted)))
  cat("\n")
}
## ▶ Executing Baseline (ITNs)...
##   ✓ Clinical cases (undiscounted): 5,791
##   ✓ Severe cases (undiscounted): 207
##   ✓ Deaths (undiscounted): 58
##   ✓ DALYs (discounted at 3%): 3,243
##   ✓ Total costs (discounted at 3%): $113,424
## 
## ▶ Executing SMC + ITNs...
##   ✓ Clinical cases (undiscounted): 3,715
##   ✓ Severe cases (undiscounted): 133
##   ✓ Deaths (undiscounted): 37
##   ✓ DALYs (discounted at 3%): 2,099
##   ✓ Total costs (discounted at 3%): $2,148,584
## 
## ▶ Executing Vaccine + ITNs...
##   ✓ Clinical cases (undiscounted): 5,110
##   ✓ Severe cases (undiscounted): 127
##   ✓ Deaths (undiscounted): 36
##   ✓ DALYs (discounted at 3%): 2,038
##   ✓ Total costs (discounted at 3%): $216,193
## 
## ▶ Executing Combined...
##   ✓ Clinical cases (undiscounted): 3,306
##   ✓ Severe cases (undiscounted): 82
##   ✓ Deaths (undiscounted): 23
##   ✓ DALYs (discounted at 3%): 1,329
##   ✓ Total costs (discounted at 3%): $3,021,914
cat("✓ All simulations completed successfully\n")
## ✓ All simulations completed successfully

Part 4: Cost-Effectiveness Analysis (Discounted)

# ===================================================================
# STEP 10: CALCULATE COST-EFFECTIVENESS METRICS
# ===================================================================
# NOTE: Costs and DALYs are discounted at 3% per annum
# Cases and deaths are undiscounted (reported as is)

ce_table <- data.frame(
  Strategy = strategy_names,
  Clinical_Cases = sapply(results_list, function(x) x$total_clinical),      # UNDISCOUNTED
  Severe_Cases = sapply(results_list, function(x) x$total_severe),          # UNDISCOUNTED
  Deaths = sapply(results_list, function(x) x$total_deaths),                # UNDISCOUNTED
  DALYs = sapply(results_list, function(x) x$total_dalys_discounted),       # DISCOUNTED
  Costs = sapply(results_list, function(x) x$total_costs_discounted)        # DISCOUNTED
)

# Calculate incremental outcomes
baseline_dalys <- ce_table$DALYs[1]
baseline_cost <- ce_table$Costs[1]

ce_table$DALYs_Averted <- baseline_dalys - ce_table$DALYs
ce_table$Inc_Cost <- ce_table$Costs - baseline_cost
ce_table$ICER <- ce_table$Inc_Cost / ce_table$DALYs_Averted
ce_table$ICER[ce_table$DALYs_Averted <= 0 | is.infinite(ce_table$ICER)] <- NA

# Display results
cat("\n\n========================================\n")
## 
## 
## ========================================
cat("COST-EFFECTIVENESS ANALYSIS RESULTS\n")
## COST-EFFECTIVENESS ANALYSIS RESULTS
cat("========================================\n")
## ========================================
cat("Clinical cases, severe cases, and deaths are UNDISCOUNTED\n")
## Clinical cases, severe cases, and deaths are UNDISCOUNTED
cat("Costs and DALYs are DISCOUNTED at 3% per annum\n")
## Costs and DALYs are DISCOUNTED at 3% per annum
cat("========================================\n\n")
## ========================================
# Create display table with proper formatting
display_table <- data.frame(
  Strategy = ce_table$Strategy,
  Clinical_Cases = format_number(ce_table$Clinical_Cases),
  Severe_Cases = format_number(ce_table$Severe_Cases),
  Deaths = format_number(ce_table$Deaths),
  DALYs = format_number(ce_table$DALYs),
  Costs = paste0("$", format_number(ce_table$Costs)),
  DALYs_Averted = format_number(ce_table$DALYs_Averted),
  ICER = ifelse(is.na(ce_table$ICER), "Dominated", paste0("$", format_number(ce_table$ICER)))
)

print(display_table)
##          Strategy Clinical_Cases Severe_Cases Deaths DALYs      Costs
## 1 Baseline (ITNs)          5,791          207     58 3,243 $  113,424
## 2      SMC + ITNs          3,715          133     37 2,099 $2,148,584
## 3  Vaccine + ITNs          5,110          127     36 2,038 $  216,193
## 4        Combined          3,306           82     23 1,329 $3,021,914
##   DALYs_Averted      ICER
## 1             0 Dominated
## 2         1,144    $1,779
## 3         1,205    $   85
## 4         1,914    $1,519
# WHO threshold interpretation
cat("\n💰 INTERPRETATION (WHO threshold: $800 per DALY averted):\n")
## 
## 💰 INTERPRETATION (WHO threshold: $800 per DALY averted):
for(i in 2:nrow(ce_table)) {
  if(!is.na(ce_table$ICER[i])) {
    if(ce_table$ICER[i] < 800) {
      cat(sprintf("  ✓ %s: Cost-effective (ICER = $%s)\n", 
                  ce_table$Strategy[i], format_number(ce_table$ICER[i])))
    } else {
      cat(sprintf("  ✗ %s: Not cost-effective (ICER = $%s)\n", 
                  ce_table$Strategy[i], format_number(ce_table$ICER[i])))
    }
  }
}
##   ✗ SMC + ITNs: Not cost-effective (ICER = $1,779)
##   ✓ Vaccine + ITNs: Cost-effective (ICER = $85)
##   ✗ Combined: Not cost-effective (ICER = $1,519)

Part 5: Geographic Subgroup Analysis (Discounted)

# ===================================================================
# STEP 11: GEOGRAPHIC SUBGROUP ANALYSIS
# ===================================================================

# Aggregate results by geography
geo_summary <- data.frame()

for(i in 1:length(strategies)) {
  geo_temp <- results_list[[i]]$results %>%
    group_by(geography) %>%
    summarise(
      Clinical_Cases = sum(clinical_cases),           # UNDISCOUNTED
      Severe_Cases = sum(severe_cases),               # UNDISCOUNTED
      Deaths = sum(deaths),                           # UNDISCOUNTED
      DALYs_Discounted = sum(dalys_discounted),       # DISCOUNTED
      Treatment_Costs_Discounted = sum(treatment_cost_discounted)  # DISCOUNTED
    ) %>%
    mutate(Strategy = strategy_names[i])
  
  geo_summary <- rbind(geo_summary, geo_temp)
}

# Calculate geographic ICERs
geo_icer_results <- data.frame()

for(geo in c("Urban", "Rural")) {
  # Extract baseline values
  base_dalys <- geo_summary[geo_summary$Strategy == strategy_names[1] & 
                              geo_summary$geography == geo, "DALYs_Discounted", drop = TRUE]
  base_cost <- geo_summary[geo_summary$Strategy == strategy_names[1] & 
                             geo_summary$geography == geo, "Treatment_Costs_Discounted", drop = TRUE]
  
  for(i in 2:length(strategies)) {
    int_dalys <- geo_summary[geo_summary$Strategy == strategy_names[i] & 
                               geo_summary$geography == geo, "DALYs_Discounted", drop = TRUE]
    int_cost <- geo_summary[geo_summary$Strategy == strategy_names[i] & 
                              geo_summary$geography == geo, "Treatment_Costs_Discounted", drop = TRUE]
    
    geo_intervention_cost <- results_list[[i]]$intervention_cost_discounted * 
                             (ifelse(geo == "Urban", 0.3, 0.7))
    
    dalys_averted <- base_dalys - int_dalys
    inc_cost <- (int_cost + geo_intervention_cost) - base_cost
    
    if(length(dalys_averted) == 1 && dalys_averted > 0 && is.finite(dalys_averted)) {
      icer <- inc_cost / dalys_averted
    } else {
      icer <- NA_real_
    }
    
    geo_icer_results <- rbind(geo_icer_results, data.frame(
      Geography = geo,
      Strategy = strategy_names[i],
      DALYs_Averted = as.numeric(dalys_averted),
      ICER = as.numeric(icer)
    ))
  }
}

# Display results
cat("\n\n========================================\n")
## 
## 
## ========================================
cat("GEOGRAPHIC SUBGROUP ANALYSIS\n")
## GEOGRAPHIC SUBGROUP ANALYSIS
cat("========================================\n")
## ========================================
cat("Clinical outcomes are UNDISCOUNTED\n")
## Clinical outcomes are UNDISCOUNTED
cat("Costs and DALYs are DISCOUNTED at 3%\n")
## Costs and DALYs are DISCOUNTED at 3%
cat("========================================\n\n")
## ========================================
for(geo in c("Urban", "Rural")) {
  cat(sprintf("\n📍 %s SETTINGS:\n", toupper(geo)))
  
  # Display clinical outcomes by geography
  geo_clinical <- geo_summary[geo_summary$geography == geo, ]
  cat("\n  Clinical Outcomes (Undiscounted):\n")
  for(j in 1:nrow(geo_clinical)) {
    cat(sprintf("    • %s: %s clinical cases, %s deaths\n", 
                geo_clinical$Strategy[j],
                format_number(geo_clinical$Clinical_Cases[j]),
                format_number(geo_clinical$Deaths[j])))
  }
  
  # Display cost-effectiveness results
  geo_subset <- geo_icer_results[geo_icer_results$Geography == geo, ]
  cat("\n  Cost-Effectiveness (Discounted):\n")
  for(j in 1:nrow(geo_subset)) {
    dalys_text <- if(!is.na(geo_subset$DALYs_Averted[j]) && geo_subset$DALYs_Averted[j] > 0) {
      format_number(geo_subset$DALYs_Averted[j])
    } else {
      "N/A"
    }
    
    icer_text <- if(!is.na(geo_subset$ICER[j]) && is.finite(geo_subset$ICER[j])) {
      paste0("$", format_number(geo_subset$ICER[j]))
    } else {
      "Dominated"
    }
    
    cat(sprintf("    • %s: DALYs averted = %s, ICER = %s\n", 
                geo_subset$Strategy[j], dalys_text, icer_text))
  }
  
  # Best strategy
  valid_idx <- which(!is.na(geo_subset$ICER) & is.finite(geo_subset$ICER) & geo_subset$DALYs_Averted > 0)
  if(length(valid_idx) > 0) {
    best_idx <- valid_idx[which.min(geo_subset$ICER[valid_idx])]
    cat(sprintf("\n  ✓ Most cost-effective: %s (ICER = $%s per DALY averted)\n", 
                geo_subset$Strategy[best_idx], 
                format_number(geo_subset$ICER[best_idx])))
  }
  cat("\n")
}
## 
## 📍 URBAN SETTINGS:
## 
##   Clinical Outcomes (Undiscounted):
##     • Baseline (ITNs): 816 clinical cases, 4 deaths
##     • SMC + ITNs: 507 clinical cases, 2 deaths
##     • Vaccine + ITNs: 687 clinical cases, 2 deaths
##     • Combined: 439 clinical cases, 1 deaths
## 
##   Cost-Effectiveness (Discounted):
##     • SMC + ITNs: DALYs averted = 84, ICER = $7,356
##     • Vaccine + ITNs: DALYs averted = 123, ICER = $274
##     • Combined: DALYs averted = 160, ICER = $5,528
## 
##   ✓ Most cost-effective: Vaccine + ITNs (ICER = $274 per DALY averted)
## 
## 
## 📍 RURAL SETTINGS:
## 
##   Clinical Outcomes (Undiscounted):
##     • Baseline (ITNs): 4,975 clinical cases, 54 deaths
##     • SMC + ITNs: 3,208 clinical cases, 35 deaths
##     • Vaccine + ITNs: 4,423 clinical cases, 34 deaths
##     • Combined: 2,867 clinical cases, 22 deaths
## 
##   Cost-Effectiveness (Discounted):
##     • SMC + ITNs: DALYs averted = 1,060, ICER = $1,337
##     • Vaccine + ITNs: DALYs averted = 1,083, ICER = $64
##     • Combined: DALYs averted = 1,755, ICER = $1,155
## 
##   ✓ Most cost-effective: Vaccine + ITNs (ICER = $64 per DALY averted)

Part 6: Visualizations

# ===================================================================
# STEP 12: CREATE VISUALIZATIONS
# ===================================================================

# Figure 1: Clinical Cases by Strategy (Undiscounted)
p1 <- ggplot(ce_table, aes(x = Strategy, y = Clinical_Cases, fill = Strategy)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_text(aes(label = format_number(Clinical_Cases)), 
            hjust = -0.2, size = 3.5) +
  labs(title = "Figure 1: Clinical Cases by Strategy (Undiscounted)",
       subtitle = "10-year horizon, 10,000 children",
       x = "", y = "Clinical Cases (Undiscounted)") +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size = 14, face = "bold")) +
  coord_flip()

# Figure 2: DALYs Averted (Discounted)
p2 <- ggplot(ce_table[2:4,], aes(x = Strategy, y = DALYs_Averted, fill = Strategy)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_text(aes(label = format_number(DALYs_Averted)), 
            hjust = -0.2, size = 4) +
  labs(title = "Figure 2: DALYs Averted by Strategy (Discounted at 3%)",
       subtitle = "Health outcomes discounted at 3% per annum",
       x = "", y = "Discounted DALYs Averted") +
  scale_y_continuous(labels = scales::comma) +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size = 14, face = "bold")) +
  coord_flip()

# Figure 3: ICER Comparison
ce_plot_data <- ce_table[2:4, ] %>% filter(!is.na(ICER))
if(nrow(ce_plot_data) > 0) {
  p3 <- ggplot(ce_plot_data, aes(x = Strategy, y = ICER, fill = Strategy)) +
    geom_bar(stat = "identity", width = 0.7) +
    geom_hline(yintercept = 800, linetype = "dashed", color = "red", size = 1.2) +
    geom_text(aes(label = sprintf("$%s", format_number(ICER))), 
              hjust = -0.2, size = 4) +
    labs(title = "Figure 3: Cost-Effectiveness Comparison (3% Discounting)",
         subtitle = "WHO threshold: $800 per DALY averted (red dashed line)",
         x = "", y = "ICER (USD per DALY averted)") +
    scale_y_continuous(labels = scales::comma) +
    theme_minimal() +
    theme(legend.position = "none",
          plot.title = element_text(size = 14, face = "bold")) +
    coord_flip()
  
}

# Figure 4: Geographic Distribution of DALYs (Discounted)
p4 <- ggplot(geo_summary, aes(x = geography, y = DALYs_Discounted/1000, fill = Strategy)) +
  geom_bar(stat = "identity", position = "dodge", width = 0.7) +
  labs(title = "Figure 4: Geographic Distribution of DALYs (Discounted)",
       subtitle = "Discounted at 3% per annum",
       x = "", y = "Discounted DALYs (thousands)") +
  scale_fill_brewer(palette = "Set1", name = "Intervention") +
  theme_minimal() +
  theme(plot.title = element_text(size = 14, face = "bold"),
        legend.position = "bottom")

# Figure 5: Age-Specific Burden (Discounted DALYs)
combined_results <- results_list[[4]]$results
p5 <- ggplot(combined_results, aes(x = age_group, y = dalys_discounted, fill = age_group)) +
  geom_bar(stat = "identity", width = 0.7) +
  geom_text(aes(label = format_number(dalys_discounted)), 
            vjust = -0.5, size = 3.5) +
  labs(title = "Figure 5: Age-Specific DALY Burden (Discounted)",
       subtitle = "Combined strategy - DALYs discounted at 3% per annum",
       x = "Age Group", y = "Discounted DALYs") +
  scale_fill_viridis_d() +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_text(size = 14, face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1))

# Display plots
print(p1)

print(p2)

print(p3)

print(p4)

print(p5)

Part 7: Age-Specific Analysis

# ===================================================================
# STEP 13: AGE-SPECIFIC RESULTS
# ===================================================================

cat("\n\n========================================\n")
## 
## 
## ========================================
cat("AGE-SPECIFIC ANALYSIS (Combined Strategy)\n")
## AGE-SPECIFIC ANALYSIS (Combined Strategy)
cat("========================================\n")
## ========================================
cat("Clinical outcomes are UNDISCOUNTED\n")
## Clinical outcomes are UNDISCOUNTED
cat("DALYs are DISCOUNTED at 3% per annum\n")
## DALYs are DISCOUNTED at 3% per annum
cat("========================================\n\n")
## ========================================
age_summary <- combined_results %>%
  group_by(age_group) %>%
  summarise(
    Clinical_Cases = sum(clinical_cases),        # UNDISCOUNTED
    Severe_Cases = sum(severe_cases),            # UNDISCOUNTED
    Deaths = sum(deaths),                        # UNDISCOUNTED
    DALYs_Discounted = sum(dalys_discounted)     # DISCOUNTED
  ) %>%
  mutate(
    Pct_Clinical = Clinical_Cases / sum(Clinical_Cases) * 100,
    Pct_DALYs = DALYs_Discounted / sum(DALYs_Discounted) * 100
  )

# Format for display
age_display <- age_summary
age_display$Clinical_Cases <- format_number(age_display$Clinical_Cases)
age_display$Severe_Cases <- format_number(age_display$Severe_Cases)
age_display$Deaths <- format_number(age_display$Deaths)
age_display$DALYs_Discounted <- format_number(age_display$DALYs_Discounted)
age_display$Pct_Clinical <- round(age_display$Pct_Clinical, 1)
age_display$Pct_DALYs <- round(age_display$Pct_DALYs, 1)

print(age_display)
## # A tibble: 4 × 7
##   age_group Clinical_Cases Severe_Cases Deaths DALYs_Discounted Pct_Clinical
##   <chr>     <chr>          <chr>        <chr>  <chr>                   <dbl>
## 1 0-6m      "  532"        "30"         "12"   "668"                    16.1
## 2 2-5y      "  681"        " 9"         " 1"   " 68"                    20.6
## 3 5+y       "1,157"        " 4"         " 0"   " 13"                    35  
## 4 6-24m     "  936"        "39"         "10"   "579"                    28.3
## # ℹ 1 more variable: Pct_DALYs <dbl>
# Identify highest burden age group
max_clinical_idx <- which.max(age_summary$Pct_Clinical)
max_daly_idx <- which.max(age_summary$Pct_DALYs)

cat("\n📊 KEY FINDINGS:\n")
## 
## 📊 KEY FINDINGS:
cat(sprintf("  • Highest clinical burden: %s (%s%% of cases)\n", 
            age_summary$age_group[max_clinical_idx],
            round(age_summary$Pct_Clinical[max_clinical_idx], 1)))
##   • Highest clinical burden: 5+y (35% of cases)
cat(sprintf("  • Highest DALY burden (discounted): %s (%s%% of DALYs)\n", 
            age_summary$age_group[max_daly_idx],
            round(age_summary$Pct_DALYs[max_daly_idx], 1)))
##   • Highest DALY burden (discounted): 0-6m (50.3% of DALYs)

Part 8: Policy Recommendations

# ===================================================================
# STEP 14: POLICY RECOMMENDATIONS
# ===================================================================

cat("\n\n========================================\n")
## 
## 
## ========================================
cat("POLICY RECOMMENDATIONS\n")
## POLICY RECOMMENDATIONS
cat("========================================\n\n")
## ========================================
# Overall most cost-effective
best_idx <- which.min(ce_table$ICER[2:4])
if(length(best_idx) > 0 && !is.na(ce_table$ICER[best_idx + 1])) {
  best_idx <- best_idx + 1
  cat("1️⃣ MOST COST-EFFECTIVE STRATEGY:\n")
  cat(sprintf("   ✓ %s\n", ce_table$Strategy[best_idx]))
  cat(sprintf("   → ICER: $%s per DALY averted\n", format_number(ce_table$ICER[best_idx])))
  
  if(ce_table$ICER[best_idx] <= 800) {
    cat("   → Cost-effective at WTP threshold (≤ $800/DALY)\n")
  } else {
    cat("   → Not cost-effective at current prices\n")
  }
}
## 1️⃣ MOST COST-EFFECTIVE STRATEGY:
##    ✓ Vaccine + ITNs
##    → ICER: $85 per DALY averted
##    → Cost-effective at WTP threshold (≤ $800/DALY)
# Geographic priorities
cat("\n2️⃣ GEOGRAPHIC PRIORITIES:\n")
## 
## 2️⃣ GEOGRAPHIC PRIORITIES:
for(geo in c("Urban", "Rural")) {
  geo_subset <- geo_icer_results[geo_icer_results$Geography == geo, ]
  valid_idx <- which(!is.na(geo_subset$ICER) & geo_subset$DALYs_Averted > 0)
  
  if(length(valid_idx) > 0) {
    best_geo <- valid_idx[which.min(geo_subset$ICER[valid_idx])]
    cat(sprintf("   • %s: %s (ICER: $%s)\n", 
                geo,
                geo_subset$Strategy[best_geo],
                format_number(geo_subset$ICER[best_geo])))
  } else {
    cat(sprintf("   • %s: No cost-effective strategy identified\n", geo))
  }
}
##    • Urban: Vaccine + ITNs (ICER: $274)
##    • Rural: Vaccine + ITNs (ICER: $64)
# Age targeting
cat("\n3️⃣ AGE TARGETING:\n")
## 
## 3️⃣ AGE TARGETING:
cat("   • Priority: 6-24 months (highest DALY burden)\n")
##    • Priority: 6-24 months (highest DALY burden)
cat("   • Secondary: 2-5 years (SMC candidates)\n")
##    • Secondary: 2-5 years (SMC candidates)
# Budget impact
cat("\n4️⃣ BUDGET IMPACT (National Scale - 5M children):\n")
## 
## 4️⃣ BUDGET IMPACT (National Scale - 5M children):
scale_factor <- 5000000 / 10000
baseline_national <- ce_table$Costs[1] * scale_factor / 1e6
combined_national <- ce_table$Costs[4] * scale_factor / 1e6
cat(sprintf("   • Baseline: $%.1f million (discounted)\n", baseline_national))
##    • Baseline: $56.7 million (discounted)
cat(sprintf("   • Combined strategy: $%.1f million (discounted)\n", combined_national))
##    • Combined strategy: $1511.0 million (discounted)
cat(sprintf("   • Additional funding: $%.1f million\n", combined_national - baseline_national))
##    • Additional funding: $1454.2 million
cat("\n5️⃣ SUMMARY OF DISCOUNTING APPROACH:\n")
## 
## 5️⃣ SUMMARY OF DISCOUNTING APPROACH:
cat("   • Costs: Discounted at 3% per annum\n")
##    • Costs: Discounted at 3% per annum
cat("   • DALYs: Discounted at 3% per annum\n")
##    • DALYs: Discounted at 3% per annum
cat("   • Clinical cases: NOT discounted (reported as raw counts)\n")
##    • Clinical cases: NOT discounted (reported as raw counts)
cat("   • Severe cases: NOT discounted (reported as raw counts)\n")
##    • Severe cases: NOT discounted (reported as raw counts)
cat("   • Deaths: NOT discounted (reported as raw counts)\n")
##    • Deaths: NOT discounted (reported as raw counts)
cat("\n\n========================================\n")
## 
## 
## ========================================
cat("MODEL COMPLETE\n")
## MODEL COMPLETE
cat("========================================\n")
## ========================================