---
title: "Inputs and Outputs Structure"
output:
  html_document:
    toc: true
    toc_float: true
    toc_depth: 5
vignette: >
  %\VignetteIndexEntry{InputsOutputsStructure}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

oldopt <- options(rmarkdown.html_vignette.check_title = FALSE)
on.exit(options(oldopt))

```


## Understanding Input Parameters

The EPIC model uses configuration files to define all input parameters. Default input values are loaded from JSON configuration files (e.g., `config_canada.json` or `config_us.json`). You can retrieve these values using `get_input()`.

### Retrieving Input Parameters

```{r, eval = TRUE, echo = TRUE}
# Get default inputs for Canada (default)
library(epicR)
input <- get_input()

# Get inputs for US
input_us <- get_input(jurisdiction = "us")

# The function returns a nested list with three components:
# - values: The actual parameter values used in simulation
# - help: Descriptions of what each parameter means
# - ref: Literature references for each parameter
```

### Input Structure Overview

The `input$values` list contains the following main categories:

| Category | Description |
|----------|-------------|
| `global_parameters` | Time horizon, age ranges, discount rates, cohort settings |
| `agent` | Population demographics, mortality, height/weight distributions |
| `smoking` | Smoking initiation, cessation, relapse rates, mortality factors |
| `COPD` | COPD prevalence and incidence equations |
| `lung_function` | FEV1 baseline, decline, and prediction equations |
| `exacerbation` | Exacerbation rates, severity, duration, and mortality |
| `symptoms` | Cough, phlegm, dyspnea, wheeze probabilities |
| `outpatient` | GP visit rates |
| `diagnosis` | Diagnosis probabilities, case detection parameters |
| `medication` | Medication adherence, effectiveness, costs, utilities |
| `cost` | Healthcare costs (exacerbations, maintenance, visits) |
| `utility` | Quality of life utilities by disease stage |
| `manual` | Advanced calibration parameters |

---

## Detailed Input Categories

### Global Parameters (`global_parameters`)

Controls overall simulation settings.

| Parameter | Description | Default (Canada) |
|-----------|-------------|------------------|
| `age0` | Starting age in the model | 40 |
| `time_horizon` | Model time horizon in years | 20 |
| `discount_cost` | Annual discount rate for costs | 0 |
| `discount_qaly` | Annual discount rate for QALYs | 0.03 |
| `closed_cohort` | 0 = open population (new agents enter), 1 = closed cohort | 0 |

```{r, eval = TRUE, echo = TRUE}
# Example: Modify global parameters
input <- get_input()
input$values$global_parameters$time_horizon <- 10
input$values$global_parameters$discount_cost <- 0.015
input$values$global_parameters$discount_qaly <- 0.015

results <- simulate(input = input$values)
```

---

### Agent Demographics (`agent`)

Defines the simulated population characteristics.

| Parameter | Description |
|-----------|-------------|
| `p_female` | Proportion of females in the population (default: 0.5) |
| `height_0_betas` | Regression coefficients for baseline height (meters) |
| `height_0_sd` | Standard deviation of baseline height |
| `weight_0_betas` | Regression coefficients for baseline weight (kg) |
| `weight_0_sd` | Standard deviation of baseline weight |
| `height_weight_rho` | Correlation between height and weight |
| `p_prevalence_age` | Age distribution at baseline (111 values, ages 0-110) |
| `p_incidence_age` | Age distribution for new population entering model |
| `p_bgd_by_sex` | Life table - annual probability of background death by age and sex |
| `l_inc_betas` | Log-rate of new population incidence |
| `ln_h_bgd_betas` | Modifiers for background mortality (time trends, etc.) |

```{r, eval = TRUE, echo = TRUE}
# Example: Change population sex distribution
input <- get_input()
input$values$agent$p_female <- 0.52  # 52% female

results <- simulate(input = input$values)
```

---

### Smoking Parameters (`smoking`)

Controls smoking behavior and its health effects.

| Parameter | Description |
|-----------|-------------|
| `logit_p_current_smoker_0_betas` | Probability of being current smoker at creation |
| `logit_p_never_smoker_con_not_current_0_betas` | Probability of being never-smoker (given not current) |
| `minimum_smoking_prevalence` | Floor for smoking prevalence |
| `mortality_factor_current` | Mortality multipliers for current smokers by age group |
| `mortality_factor_former` | Mortality multipliers for former smokers by age group |
| `pack_years_0_betas` | Pack-years regression at baseline |
| `pack_years_0_sd` | Standard deviation of pack-years |
| `ln_h_inc_betas` | Log-hazard of smoking initiation/relapse |
| `ln_h_ces_betas` | Log-hazard of smoking cessation |
| `smoking_ces_coefficient` | Decay rate for cessation treatment effect |
| `smoking_cessation_adherence` | Proportion adhering to cessation treatment |

```{r, eval = TRUE, echo = TRUE}
# Example: Increase smoking cessation rate
input <- get_input()
# The first coefficient is the intercept - increasing it increases cessation rate
input$values$smoking$ln_h_ces_betas[1] <- input$values$smoking$ln_h_ces_betas[1] + 0.5

results <- simulate(input = input$values)
```

---

### COPD Parameters (`COPD`)

Defines COPD prevalence and incidence.

| Parameter | Description |
|-----------|-------------|
| `logit_p_COPD_betas_by_sex` | Logit-probability of having COPD at creation (by sex) |
| `ln_h_COPD_betas_by_sex` | Log-hazard of developing COPD over time (by sex) |

The beta coefficients typically include:

- Intercept
- Age effect
- Age-squared effect
- Pack-years effect
- Current smoking effect
- Calendar time effect

---

### Lung Function Parameters (`lung_function`)

Controls FEV1 modeling.

| Parameter | Description |
|-----------|-------------|
| `fev1_0_prev_betas_by_sex` | FEV1 at creation for prevalent COPD cases |
| `fev1_0_prev_sd_by_sex` | SD of FEV1 at creation (prevalent) |
| `fev1_0_inc_betas_by_sex` | FEV1 at COPD onset for incident cases |
| `fev1_0_inc_sd_by_sex` | SD of FEV1 at COPD onset |
| `fev1_0_ZafarCMAJ_by_sex` | Coefficients for Zafari CMAJ FEV1 decline model |
| `pred_fev1_betas_by_sex` | Predicted FEV1 calculation coefficients |
| `fev1_betas_by_sex` | Mixed-effects model for FEV1 decline |
| `dfev1_sigmas` | Random effect variances for FEV1 decline |
| `dfev1_re_rho` | Correlation between random effects |

---

### Exacerbation Parameters (`exacerbation`)

Controls exacerbation occurrence and outcomes.

| Parameter | Description |
|-----------|-------------|
| `ln_rate_betas` | Log-rate of exacerbation (10 coefficients) |
| `ln_rate_intercept_sd` | Random intercept SD for exacerbation rate |
| `logit_severity_betas` | Proportional odds model for severity (9 coefficients) |
| `logit_severity_intercept_sd` | Random intercept SD for severity |
| `rate_severity_intercept_rho` | Correlation between rate and severity random effects |
| `exac_end_rate` | Rate of exacerbation resolution by severity (4 values) |
| `logit_p_death_by_sex` | Probability of death during exacerbation |

**Exacerbation severity levels:**

| Level | Severity |
|-------|----------|
| 1 | Mild |
| 2 | Moderate |
| 3 | Severe |
| 4 | Very severe |

```{r, eval = TRUE, echo = TRUE}
# Example: Reduce exacerbation rate
input <- get_input()
# Reduce the baseline rate (first coefficient is intercept)
input$values$exacerbation$ln_rate_betas[1] <- input$values$exacerbation$ln_rate_betas[1] - 0.2

results <- simulate(input = input$values)
```

---

### Symptom Parameters (`symptoms`)

Controls respiratory symptoms.

| Parameter | Description |
|-----------|-------------|
| `logit_p_cough_COPD_by_sex` | Cough probability for COPD patients |
| `logit_p_cough_nonCOPD_by_sex` | Cough probability for non-COPD |
| `logit_p_phlegm_COPD_by_sex` | Phlegm probability for COPD patients |
| `logit_p_phlegm_nonCOPD_by_sex` | Phlegm probability for non-COPD |
| `logit_p_dyspnea_COPD_by_sex` | Dyspnea probability for COPD patients |
| `logit_p_dyspnea_nonCOPD_by_sex` | Dyspnea probability for non-COPD |
| `logit_p_wheeze_COPD_by_sex` | Wheeze probability for COPD patients |
| `logit_p_wheeze_nonCOPD_by_sex` | Wheeze probability for non-COPD |
| `covariance_COPD` | 4x4 covariance matrix for symptom random effects (COPD) |
| `covariance_nonCOPD` | 4x4 covariance matrix for symptom random effects (non-COPD) |

---

### Outpatient Parameters (`outpatient`)

Controls GP visit frequency.

| Parameter | Description |
|-----------|-------------|
| `ln_rate_gpvisits_COPD_by_sex` | GP visit rate for COPD patients |
| `dispersion_gpvisits_COPD` | Dispersion parameter for COPD GP visits |
| `ln_rate_gpvisits_nonCOPD_by_sex` | GP visit rate for non-COPD |
| `dispersion_gpvisits_nonCOPD` | Dispersion parameter for non-COPD GP visits |

---

### Diagnosis Parameters (`diagnosis`)

Controls COPD diagnosis and case detection.

| Parameter | Description |
|-----------|-------------|
| `p_case_detection` | Vector of case detection probabilities by year |
| `case_detection_start_end_yrs` | Years when case detection program is active |
| `years_btw_case_detection` | Interval between case detection opportunities |
| `min_cd_age` | Minimum age for case detection |
| `min_cd_pack_years` | Minimum pack-years for case detection eligibility |
| `min_cd_symptoms` | Require symptoms for case detection (0 or 1) |
| `case_detection_methods` | Sensitivity, specificity, cost for each method |
| `case_detection_methods_eversmokers` | Methods for ever-smoker population |
| `case_detection_methods_symptomatic` | Methods for symptomatic population |
| `logit_p_prevalent_diagnosis_by_sex` | Probability of diagnosis at baseline |
| `logit_p_diagnosis_by_sex` | Probability of diagnosis during follow-up |
| `p_hosp_diagnosis` | Probability of diagnosis during hospitalization |
| `logit_p_overdiagnosis_by_sex` | Probability of false COPD diagnosis |
| `p_correct_overdiagnosis` | Probability of correcting overdiagnosis |

```{r, eval = TRUE, echo = TRUE}
# Example: Enable case detection program
input <- get_input()
input$values$diagnosis$case_detection_start_end_yrs <- c(0, 20)  # Active years 0-20
input$values$diagnosis$p_case_detection <- rep(0.1, 20)  # 10% annual probability

results <- simulate(input = input$values)
```

---

### Medication Parameters (`medication`)

Controls medication use and effects.

| Parameter | Description |
|-----------|-------------|
| `medication_adherence` | Proportion adhering to medication (default: 0.7) |
| `medication_ln_hr_exac` | Log hazard ratio for exacerbation reduction by treatment |
| `medication_costs` | Annual costs by medication class |
| `medication_utility` | Utility gain by medication class |

**Medication Classes:**

| Class | Code | Description |
|-------|------|-------------|
| `none` | 0 | No medication |
| `SABA` | 1 | Short-acting beta-agonist |
| `LAMA` | 4 | Long-acting muscarinic antagonist |
| `LAMA/LABA` | 6 | Long-acting muscarinic antagonist/Long-acting beta-agonist |
| `ICS/LAMA/LABA` | 14 | Triple Therapy |


Combinations are represented as sums (e.g., ICS+LAMA+LABA = 8+4+2 = 14).

---

### Cost Parameters (`cost`)

All costs are in local currency (CAD for Canada, USD for US).

| Parameter | Description |
|-----------|-------------|
| `bg_cost_by_stage` | Annual maintenance costs by GOLD stage (0=no COPD, 1-4=GOLD stages) |
| `exac_dcost` | Incremental cost per exacerbation by severity (4 values) |
| `cost_case_detection` | Cost of case detection procedure |
| `cost_outpatient_diagnosis` | Cost of outpatient diagnosis |
| `cost_gp_visit` | Cost per GP visit |
| `cost_smoking_cessation` | Cost of smoking cessation treatment |

```{r, eval = TRUE, echo = TRUE}
# Example: View and modify costs
input <- get_input()

# View default exacerbation costs
print(input$values$cost$exac_dcost)
# Columns: mild, moderate, severe, very_severe

# Double the cost of severe exacerbations
input$values$cost$exac_dcost[3] <- input$values$cost$exac_dcost[3] * 2

results <- simulate(input = input$values)
```

---

### Utility Parameters (`utility`)

Quality of life values (0-1 scale, where 1 = perfect health).

| Parameter | Description |
|-----------|-------------|
| `bg_util_by_stage` | Background utility by GOLD stage (5 values: no COPD + 4 stages) |
| `exac_dutil` | Utility decrement during exacerbation (4x4 matrix: severity x GOLD) |

```{r, eval = TRUE, echo = TRUE}
# Example: View utility values
input <- get_input()

# Background utilities by stage
print(input$values$utility$bg_util_by_stage)
# Typical values: No COPD=0.86, GOLD1=0.81, GOLD2=0.72, GOLD3=0.68, GOLD4=0.58

# Exacerbation disutility matrix
print(input$values$utility$exac_dutil)
```

---

### Manual/Advanced Parameters (`manual`)

Calibration and advanced settings.

| Parameter | Description |
|-----------|-------------|
| `MORT_COEFF` | Global mortality coefficient multiplier |
| `smoking$intercept_k` | Smoking intercept adjustment |
| `explicit_mortality_by_age_sex` | COPD-attributable mortality to subtract from background |

---

## Accessing Parameter Help and References

```{r, eval = TRUE, echo = TRUE}
input <- get_input()

# Get help text for a parameter
input$help$exacerbation$ln_rate_betas
# Returns: "Regression coefficients for the random-effects log-hazard model..."

# Get literature reference
input$ref$smoking$mortality_factor_current
# Returns: "Meta-analysis. doi:10.1001/archinternmed.2012.1397"

# List all help for a category
names(input$help$cost)
```

---

## Running Simulations with Custom Inputs

```{r, eval = TRUE, echo = TRUE}
# Complete workflow example
input <- get_input(jurisdiction = "canada")

# Modify multiple parameters
input$values$global_parameters$time_horizon <- 10
input$values$agent$p_female <- 0.52
input$values$smoking$ln_h_ces_betas[1] <- input$values$smoking$ln_h_ces_betas[1] + 0.3

# Run simulation
results <- simulate(input = input$values, n_agents = 50000)

# Access results
print(results$basic)
```

---

## Output Structure

The `simulate()` function returns a list with simulation results.

### Basic Output (`results$basic`)

Summary statistics from the simulation.

| Output | Description |
|--------|-------------|
| `n_agents` | Total number of agents simulated |
| `cumul_time` | Cumulative person-years of follow-up |
| `n_deaths` | Total number of deaths |
| `n_COPD` | Number of agents with COPD at end of simulation |
| `total_exac` | Vector of total exacerbations by severity (4 values) |
| `total_exac_time` | Cumulative time spent in exacerbation by severity |
| `total_pack_years` | Total pack-years across all agents |
| `total_cost` | Total discounted costs |
| `total_qaly` | Total discounted QALYs |
| `total_diagnosed_time` | Person-years with COPD diagnosis |

```{r, eval = TRUE, echo = TRUE}
results <- simulate(n_agents = 50000)

# Access basic results
basic <- results$basic
cat("Total agents:", basic$n_agents, "\n")
cat("Total deaths:", basic$n_deaths, "\n")
cat("Mild exacerbations:", basic$total_exac[1], "\n")
cat("Moderate exacerbations:", basic$total_exac[2], "\n")
cat("Severe exacerbations:", basic$total_exac[3], "\n")
cat("Very severe exacerbations:", basic$total_exac[4], "\n")
cat("Total cost:", basic$total_cost, "\n")
cat("Total QALYs:", basic$total_qaly, "\n")
```

---

### Extended Output (`results$extended`)

Detailed breakdowns by time, age, sex, and disease stage. Extended results are included by default.
Note: `ctime = calendar time`, `ltime = local time`

#### Population Counts by Time

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_alive_by_ctime_sex` | time_horizon x 2 | Alive agents by calendar time and sex |
| `n_alive_by_ctime_age` | time_horizon x 111 | Alive agents by calendar time and age |
| `n_alive_by_age_sex` | 111 x 2 | Alive agents by age and sex (total) |

#### Smoking Status by Time

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_smoking_status_by_ctime` | time_horizon x 3 | Agents by smoking status (never/current/former) |
| `n_current_smoker_by_ctime_sex` | time_horizon x 2 | Current smokers by time and sex |
| `cumul_time_by_smoking_status` | 3 | Person-years by smoking status |

#### COPD Outcomes by Time


| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_COPD_by_ctime_sex` | time_horizon x 2 | COPD cases by time and sex |
| `n_COPD_by_ctime_age` | time_horizon x 111 | COPD cases by time and age |
| `n_COPD_by_ctime_severity` | time_horizon x 5 | COPD cases by time and GOLD stage (0-4) |
| `n_COPD_by_age_sex` | 111 x 2 | COPD cases by age and sex (total) |
| `n_inc_COPD_by_ctime_age` | time_horizon x 111 | Incident COPD cases by time and age |
| `cumul_non_COPD_time` | scalar | Total person-years without COPD |
| `cumul_time_by_ctime_GOLD` | time_horizon x 5 | Person-years by time and GOLD stage. GOLD0 = no COPD |

#### Diagnosis Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_Diagnosed_by_ctime_sex` | time_horizon x 2 | Diagnosed COPD cases by time and sex |
| `n_Diagnosed_by_ctime_severity` | time_horizon x 4 | Diagnosed cases by time and GOLD stage |
| `n_Overdiagnosed_by_ctime_sex` | time_horizon x 2 | False diagnoses by time and sex |
| `n_case_detection_by_ctime` | time_horizon x 4 | Case detection events by time and method |
| `n_case_detection_eligible` | scalar | Total eligible for case detection |
| `n_diagnosed_true_CD` | scalar | True positives from case detection |
| `n_agents_CD` | scalar | Agents who received case detection |

#### Exacerbation Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_exac_by_ctime_age` | time_horizon x 111 | Exacerbations by time and age |
| `n_exac_by_ctime_severity` | time_horizon x 4 | Exacerbations by time and severity |
| `n_exac_by_ctime_sex` | time_horizon x 2 | Exacerbations by time and sex |
| `n_exac_by_ctime_GOLD` | time_horizon x 4 | Exacerbations by time and GOLD stage |
| `n_exac_by_gold_severity` | 4 x 4 | Exacerbations by GOLD stage and severity |
| `n_exac_by_gold_severity_diagnosed` | 4 x 4 | Exacerbations among diagnosed patients |
| `n_severep_exac_by_ctime_age` | time_horizon x 111 | Severe and very severe exacerbations by time and age |
| `n_exac_by_ctime_severity_undiagnosed` | time_horizon x 4 | Exacerbations in undiagnosed patients |
| `n_exac_by_ctime_severity_diagnosed` | time_horizon x 4 | Exacerbations in diagnosed patients |
| `n_exac_by_ctime_severity_female` | time_horizon x 4 | Exacerbations in females by severity |

#### Death Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_death_by_age_sex` | 111 x 2 | Deaths by age and sex |
| `n_exac_death_by_ctime_age` | time_horizon x 111 | Exacerbation deaths by time and age |
| `n_exac_death_by_ctime_severity` | time_horizon x 4 | Exacerbation deaths by severity |
| `n_exac_death_by_age_sex` | 111 x 2 | Exacerbation deaths by age and sex |

#### Cost and QALY Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `annual_cost_ctime` | time_horizon | Annual costs by calendar time |
| `cumul_cost_ctime` | time_horizon | Cumulative costs by calendar time |
| `cumul_cost_gold_ctime` | time_horizon x 5 | Cumulative costs by time and GOLD stage |
| `cumul_qaly_ctime` | time_horizon | Cumulative QALYs by calendar time |
| `cumul_qaly_gold_ctime` | time_horizon x 5 | Cumulative QALYs by time and GOLD stage |

#### Additional Summary Statistics

| Output | Dimensions | Description |
|--------|------------|-------------|
| `sum_fev1_ltime` | time_horizon | Sum of FEV1 values by local time |
| `sum_p_COPD_by_ctime_sex` | time_horizon x 2 | Sum of COPD probabilities by time/sex |
| `sum_pack_years_by_ctime_sex` | time_horizon x 2 | Sum of pack-years by time and sex |
| `sum_age_by_ctime_sex` | time_horizon x 2 | Sum of ages by time and sex |
| `sum_time_by_ctime_sex` | time_horizon x 2 | Person-time by calendar time and sex |
| `sum_time_by_age_sex` | 111 x 2 | Person-time by age and sex |
| `sum_weight_by_ctime_sex` | time_horizon x 2 | Sum of weights by time and sex |

#### Symptom and GP Visit Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `n_GPvisits_by_ctime_sex` | time_horizon x 2 | GP visits by time and sex |
| `n_GPvisits_by_ctime_severity` | time_horizon x 5 | GP visits by time and GOLD stage |
| `n_GPvisits_by_ctime_diagnosis` | time_horizon x 2 | GP visits by diagnosis status |
| `n_cough_by_ctime_severity` | time_horizon x 5  | Number with cough by time and GOLD |
| `n_phlegm_by_ctime_severity` | time_horizon x 5 | Number with phlegm by time and GOLD |
| `n_wheeze_by_ctime_severity` | time_horizon x 5 | Number with wheeze by time and GOLD |
| `n_dyspnea_by_ctime_severity` | time_horizon x 5 | Number with dyspnea by time and GOLD |

#### Medication Outcomes

| Output | Dimensions | Description |
|--------|------------|-------------|
| `medication_time_by_class` | 5 | Person-time on each medication class |
| `medication_time_by_ctime_class` | time_horizon x 5 | Medication time by year and class |
| `n_smoking_cessation_by_ctime` | time_horizon | Smoking cessation treatments by time |
| `n_exac_by_medication_class` | 5 x 4 | Exacerbations by medication and severity |

```{r, eval = TRUE, echo = TRUE}
results <- simulate(n_agents = 50000)

# Access extended results
extended <- results$extended

# Calculate COPD prevalence over time
copd_prevalence <- rowSums(extended$n_COPD_by_ctime_sex) /
                   rowSums(extended$n_alive_by_ctime_sex)
plot(copd_prevalence, type = "l",
     xlab = "Year", ylab = "COPD Prevalence",
     main = "COPD Prevalence Over Time")

# Calculate annual exacerbation rates by severity
annual_exac_rates <- extended$n_exac_by_ctime_severity /
                     rowSums(extended$n_COPD_by_ctime_sex)

# Analyze costs by GOLD stage
costs_by_gold <- colSums(extended$cumul_cost_gold_ctime)
names(costs_by_gold) <- c("No COPD", "GOLD 1", "GOLD 2", "GOLD 3", "GOLD 4")
barplot(costs_by_gold, main = "Total Costs by GOLD Stage")
```

---

## Event History (`results$events`)

For detailed agent-level analysis, you can capture the complete event history by setting `return_events = TRUE`.

**Note:** Event recording requires significant memory. Use with caution for large simulations.

```{r, eval = TRUE, echo = TRUE}
# Run simulation with event recording
results <- simulate(
  n_agents = 10000,
  return_events = TRUE
)

# Access events data frame
events <- results$events
head(events)
```

### Event Matrix Columns (33 columns)

| Column | Description |
|--------|-------------|
| `id` | Unique agent identifier |
| `local_time` | Time since agent creation (years) |
| `female` | Sex (0=male, 1=female) |
| `time_at_creation` | Calendar time when agent entered model |
| `age_at_creation` | Age when agent entered model |
| `pack_years` | Cumulative pack-years of smoking |
| `gold` | GOLD stage (0=no COPD, 1-4=GOLD stages) |
| `event` | Event type code (see table below) |
| `FEV1` | Current FEV1 value (liters) |
| `FEV1_slope` | FEV1 decline slope |
| `FEV1_acceleration` | FEV1 decline acceleration |
| `pred_FEV1` | Predicted FEV1 based on age/height/sex |
| `smoking_status` | Current smoking status (0=not smoking, 1=smoking) |
| `localtime_at_COPD` | Local time when COPD developed |
| `age_at_COPD` | Age when COPD developed |
| `weight_at_COPD` | Weight when COPD developed |
| `height` | Height (meters) |
| `followup_after_COPD` | Time since COPD onset |
| `FEV1_baseline` | FEV1 at COPD onset |
| `exac_status` | Current exacerbation status (0=none, 1-4=severity) |
| `cough` | Cough symptom present (0/1) |
| `phlegm` | Phlegm symptom present (0/1) |
| `wheeze` | Wheeze symptom present (0/1) |
| `dyspnea` | Dyspnea symptom present (0/1) |
| `gpvisits` | Annual number of GP visits |
| `diagnosis` | Diagnosis status (0=undiagnosed, 1=diagnosed) |
| `medication_status` | Medication class bitfield |
| `case_detection` | Case detection received (0/1) |
| `cumul_cost` | Cumulative cost for this agent |
| `cumul_qaly` | Cumulative QALYs for this agent |
| `time_at_diagnosis` | Local time when diagnosed |
| `exac_history_n_moderate` | Cumulative count of moderate exacerbations |
| `exac_history_n_severe_plus` | Cumulative count of severe and very severe exacerbations |

### Event Type Codes

| Event | Code | Description |
|-------|------|-------------|
| Start | 0 | Agent creation |
| Annual | 1 | Annual update (fixed time event) |
| Birthday | 2 | Birthday event (currently disabled) |
| Smoking change | 3 | Smoking status change (start or quit) |
| COPD incidence | 4 | New COPD diagnosis |
| Exacerbation | 5 | Exacerbation onset |
| Exacerbation end | 6 | Exacerbation resolution |
| Death by exacerbation | 7 | Death during exacerbation |
| Doctor visit | 8 | Doctor visit (not implemented) |
| Medication change | 9 | Medication change (not implemented) |
| Background death | 13 | Death from background causes |
| End | 14 | End of follow-up |

```{r, eval = TRUE, echo = TRUE}
# Analyze events for a specific agent
results <- simulate(n_agents = 10000, return_events = TRUE)
events <- results$events

# Get all events for agent #5
agent_5_events <- events[events$id == 5, ]
print(agent_5_events[, c("local_time", "event", "gold", "exac_status")])

# Count events by type
event_counts <- table(events$event)
names(event_counts) <- c("Start", "Annual", "Birthday", "Smoking",
                         "COPD", "Exac", "Exac_end", "Exac_death",
                         "Doctor", "Med_change", "", "", "",
                         "Bgd_death", "End")[as.numeric(names(event_counts)) + 1]
print(event_counts)

# Find all exacerbation events
exac_events <- events[events$event == 5, ]
cat("Total exacerbations:", nrow(exac_events), "\n")
cat("By severity:\n")
print(table(exac_events$exac_status))
```

---

## Record Modes

The `record_mode` setting controls what level of detail is captured.

| Mode | Value | Description | Memory Usage |
|------|-------|-------------|--------------|
| `record_mode_none` | 0 | Aggregate output only | Minimal |
| `record_mode_agent` | 1 | Agent-level summaries | Low |
| `record_mode_event` | 2 | Complete event history | High |
| `record_mode_some_event` | 3 | Selected event types only | Medium |

```{r, eval = TRUE, echo = TRUE}
# Using settings directly
settings <- get_default_settings()
settings$record_mode <- 2  # Full event recording
settings$n_base_agents <- 10000

results <- simulate(settings = settings, return_events = TRUE)

# Or simply use return_events = TRUE (automatically sets record_mode = 2)
results <- simulate(n_agents = 10000, return_events = TRUE)
```

---

## Complete Example: Cost-Effectiveness Analysis

```{r, eval = TRUE, echo = TRUE}
# Baseline scenario
input_baseline <- get_input(jurisdiction = "canada")
results_baseline <- simulate(
  input = input_baseline$values,
  n_agents = 50000,
  seed = 12345  # For reproducibility
)

# Intervention scenario: Enhanced smoking cessation
input_intervention <- get_input(jurisdiction = "canada")
input_intervention$values$smoking$ln_h_ces_betas[1] <-
  input_intervention$values$smoking$ln_h_ces_betas[1] + 0.5

results_intervention <- simulate(
  input = input_intervention$values,
  n_agents = 50000,
  seed = 12345  # Same seed for comparability
)

# Calculate incremental cost-effectiveness
delta_cost <- results_intervention$basic$total_cost - results_baseline$basic$total_cost
delta_qaly <- results_intervention$basic$total_qaly - results_baseline$basic$total_qaly
icer <- delta_cost / delta_qaly

cat("Baseline total cost:", results_baseline$basic$total_cost, "\n")
cat("Intervention total cost:", results_intervention$basic$total_cost, "\n")
cat("Baseline total QALYs:", results_baseline$basic$total_qaly, "\n")
cat("Intervention total QALYs:", results_intervention$basic$total_qaly, "\n")
cat("Incremental cost:", delta_cost, "\n")
cat("Incremental QALYs:", delta_qaly, "\n")
cat("ICER:", icer, "per QALY\n")
```

---

## Important Notes

- **Starting Age:** EPIC only simulates patients starting at age 40
- **Comorbidities:** Cardiovascular comorbidities (MI, stroke, heart failure) are not implemented in the current version
- **Memory:** Event recording can require substantial memory for large simulations
- **Reproducibility:** Use the `seed` parameter for reproducible results
- **Discounting:** Costs and QALYs are discounted at the rates specified in `global_parameters`