Asthma Exacerbations Model¶

Data¶

Population Data¶

We use the Statistics Canada population data that was generated and saved as: processed_data/{time_delta_tag}/birth/initial_population.csv. This file is produced by the same data-generation script as birth_estimate.csv — see Birth Model for details on how it is generated from Statistics Canada population data.

Column	Type	Description
`timepoint`	`datetime`	the starting date / time of the time interval that the data applies to
`age`	`int`	the age of the person in years
`province`	`str`	the province of the person (e.g., `AB` = Alberta, `BC` = British Columbia, etc.)
`n_age`	`int`	the number of people in a given age group, time interval, province, and projection scenario
`n_birth`	`int`	the number of births in that time interval, province, and projection scenario
`prop`	`float`	the proportion of the population in that age group, time interval, province, and projection scenario relative to the number of births in that time interval, province, and projection scenario
`prop_male`	`float`	the proportion of the population in a given age group, time interval, province, and projection scenario who are male
`projection_scenario`	`str`	the projection scenario used to generate the data

This dataset has no sex column, so wherever a sex-specific population count \(N(a, s, t)\) is needed elsewhere in this document, it is derived from age \(a\) and timepoint \(t\)’s n_age and prop_male values:

\[\begin{split}N(a, \text{M}, t) &= n_{\text{age}} \cdot p_{\text{male}} \\ N(a, \text{F}, t) &= n_{\text{age}} \cdot (1 - p_{\text{male}})\end{split}\]

where \(n_{\text{age}}\) and \(p_{\text{male}}\) are the n_age and prop_male columns respectively.

Occurrence Data¶

We use the occurrence data that was generated and saved as: processed_data/{time_delta_tag}/asthma_occurrence_predictions.csv

See Occurrence Model 1: Crude Occurrence for more details about this dataset.

Column	Type	Description
`timepoint`	`datetime`	The start of the time interval, e.g. 2024-01-01
`sex`	`str`	`F` = female, `M` = male
`age`	`int`	the age of the patient in years
`incidence`	`float`	the predicted asthma incidence for the given time interval, age, and sex
`prevalence`	`float`	the predicted asthma prevalence for the given time interval, age, and sex

Hospitalization Data¶

The data is from the Hospital Morbidity Database (HMDB) from the Canadian Institute for Health Information (CIHI).

The hospitalization data was collected from patients presenting to a hospital in Canada due to an asthma exacerbation. We will use this data to calibrate the exacerbation model.

Per-province Structure¶

The raw data is saved under original_data/asthma_hosp/{province}/, with one subfolder per province/territory: AB, BC, MB, NB, NL, NS, ON, PE, QC, SK. In addition there are two combined regions:

CA: all of Canada combined.
TR: Nunavut, Northwest Territories, and Yukon combined into a single category, since each of these territories individually has too few cases to report reliably.

Each province subfolder contains 5 files that share the same columns and shape, differing only in what value is reported. Of these, only tab1_rate.csv (the hospitalization rate per 100 000 people) is used to calibrate the exacerbation model (see below); the remaining 4 files (tab1_count.csv, tab1_N.csv, tab1_lower.csv, tab1_upper.csv) and the per-province los.csv file are documented in Raw Data Files in the Developers section.

`tab1_rate.csv` Columns¶

Column	Type	Description
`fiscal_year`	`int`	the year the data was collected
`N`	`float`	the hospitalization rate for all ages and sexes in that year.
`M`	`float`	the rate for all ages who are male in that year.
`F`	`float`	the rate for all ages who are female in that year.
`0`	`float`	the rate for all sexes aged 0 in that year.
`1`	`float`	the rate for all sexes aged 1 in that year.
`...`	`...`	…
`90`	`float`	the rate for all sexes aged 90 in that year.
`90+`	`float`	the rate for all sexes aged over 90 in that year.
`F_0`	`float`	the rate for all females aged 0 in that year.
`...`	`...`	…
`F_90+`	`float`	the rate for all females aged over 90 in that year.
`M_0`	`float`	the rate for all males aged 0 in that year.
`...`	`...`	…
`M_90+`	`float`	the rate for all males aged over 90 in that year.

The hospitalization rate in this table is the hospitalization rate per 100 000 people. For example, in the category F_90+, the value would be the rate for people hospitalized who are female and over 90 during the given year.

Therefore, the observed number of hospitalizations for a given age \(a\), sex \(s\), and timepoint \(t\) can be recovered from the rate:

\[N_{\text{hosp}}(a, s, t) = \dfrac{\text{hospitalization rate}(a, s, t)}{100\,000} \cdot N(a, s, t)\]

Here, \(N(a, s, t)\) is the population count from the Population Data described above.

This is used to calibrate \(\alpha\) below.

The calibration step only uses the per-sex, per-age columns (F_0 … F_90+ and M_0 … M_90+); the sex- and age-aggregated columns (N, M, F, 0 … 90+) are not used.

Province Coverage¶

Although hospitalization data is available for every province/territory listed above, calibration is currently only implemented for BC and CA using the corresponding tab1_rate.csv file (see leap/data_generation/exacerbation_data.py). The resulting exacerbation_calibration.csv used at runtime therefore only contains calibration multipliers for these two regions.

Model¶

The number of exacerbations in a given time interval is modelled using a Poisson distribution. The formula is:

\[N_{\text{exacerbations}}^{(i)} \sim \text{Poisson}(\lambda^{(i)}) = \dfrac{{\lambda^{(i)}}^{\kappa} e^{-\lambda^{(i)}}}{\kappa!}\]

Here \(\lambda^{(i)}\) is the expected number of exacerbations per time interval for agent \(i\), and \(\kappa\) is the number of exacerbations (a non-negative integer) for which we are computing the probability. To obtain \(\lambda^{(i)}\), we must perform a Poisson regression. The Poisson regression assumes that the value we are interested in can be approximated using the following formula:

\[\ln(\lambda^{(i)}) = \ln(\alpha) + \beta_0^{(i)} + \sum_{k=1}^3 \beta_k c_k^{(i)}\]

where:

Variable	Description
\(\alpha\)	the calibration multiplier that adjusts the model to match the hospitalization data
\(\beta_0^{(i)}\)	patient-specific random effect; \(\beta_0^{(i)} \sim \mathcal{N}(0, \sigma^2)\)
\(c_k^{(i)}\)	relative time spent in control level \(k\), given by the probability of control level \(k\) from the Control Model
\(\beta_k\)	control level constant, derived from the EBA and GOAL studies (see Calibration below)

For each agent with asthma, \(\beta_0^{(i)}\) is sampled once and held fixed for their simulated lifetime, representing individual heterogeneity in exacerbation risk beyond what is explained by age, sex, control level, and the population-level calibration \(\alpha\). The mean and variance of this distribution are configuration-dependent; in the current default configuration they are set to approximately \(\mathcal{N}(0, 0)\), effectively disabling this source of individual variation.

This gives us the rate \(\lambda^{(i)}\) of exacerbations of any severity. Once an individual’s total number of exacerbations for a time interval is drawn from this Poisson model, each exacerbation is further assigned a severity level (mild, moderate, severe, or very severe) — see the Asthma Exacerbation Severity Model.

The diagram below summarises how these pieces fit together conceptually, for agent \(i\) with a given age \(a\) and sex \(s\), from calibration through to the final severity-level counts.

flowchart TD
    ALPHA["<b>Calibration multiplier</b>&nbsp;<span style='font-size:1.3em'>$$\alpha$$</span><br/>predicted vs. observed<br/>hospitalizations in CIHI data"]
    BI["<b>Control-level rate constants</b>&nbsp;<span style='font-size:1.3em'>$$\beta_k$$</span><br/>derived from EBA + GOAL studies"]
    CI["<b>Control-level probabilities</b>&nbsp;<span style='font-size:1.3em'>$$c_k^{(i)}$$</span><br/>predicted by the Control Model"]

    RATE["<b>Rate of exacerbations</b>&nbsp;<span style='font-size:1.3em'>$$\lambda^{(i)}$$</span><br/>(any severity)"]

    ALPHA --> RATE
    BI --> RATE
    CI --> RATE

    COUNT["<b>Total exacerbations</b>&nbsp;<span style='font-size:1.3em'>$$N_{\text{exacerbations}}^{(i)}$$</span><br/>this time interval"]
    RATE --> COUNT

    SEV["<b>Severity probabilities</b>&nbsp;<span style='font-size:1.3em'>$$\mathbf{w}^{(i)}$$</span><br/>based on SYGMA II<br/>severity/hospitalization proportions"]
    HIST["<b>History of very severe</b><br/><b>exacerbations</b>&nbsp;<span style='font-size:1.3em'>$$\beta_{\text{prev hosp}}^{(i)}$$</span><br/>(hospitalization)"]
    HIST -->|"increases probability<br/>of very severe"| SEV

    OUTPUT["<b>Exacerbations by</b><br/><b>severity level</b>&nbsp;<span style='font-size:1.3em'>$$(n_{\text{mild}}^{(i)}, \ldots, n_{\text{very severe}}^{(i)})$$</span>"]
    COUNT --> OUTPUT
    SEV --> OUTPUT

    classDef input fill:#fff3e0,stroke:#e65100,color:#3a2400;
    classDef stage fill:#e3f2fd,stroke:#1565c0,color:#0d2b45;
    classDef output fill:#e8f5e9,stroke:#2e7d32,color:#10300f;
    class ALPHA,BI,CI,HIST input;
    class RATE,COUNT,SEV stage;
    class OUTPUT output;

    classDef input fill:#fff3e0,stroke:#e65100,color:#3a2400;
    classDef formula fill:#e3f2fd,stroke:#1565c0,color:#0d2b45;
    classDef sim fill:#e8f5e9,stroke:#2e7d32,color:#10300f;
    class ALPHA,B0,BI,CI input;
    class FORMULA,POISSON formula;
    class W,HIST,MULTI sim;

Calibration¶

We are interested in calculating \(\alpha\). If we rewrite the equation, the meaning of \(\alpha\) becomes more apparent:

\[\lambda^{(i)} = \alpha \cdot e^{\beta_0^{(i)}} \prod_{k=1}^3 e^{\beta_k c_k^{(i)}}\]

How do we obtain \(\alpha\)? We again assume that the mean value has the same form as in a Poisson regression, with the following formula:

\[\ln(\lambda_{C}(a, s)) = \sum_{k=1}^3 \beta_k c_k\]

\(\lambda_C(a, s)\): the predicted mean number of exacerbations per year for a given age \(a\) and sex \(s\) — note this has no timepoint argument, since the Control Model assumes control level probabilities do not vary by calendar year
\(c_k\): the age- and sex-specific probability of control level \(k\), \(P(y^{(i)} = k)\), from the Control Model’s ordinal regression — obtained by subtracting consecutive cumulative probabilities, \(P(y^{(i)} = k) = P(y^{(i)} \leq k) - P(y^{(i)} \leq k-1)\)
\(\beta_k\): control level constant, derived below from the EBA and GOAL studies

The \(\beta_k\) values are derived by combining two literature sources, as described in [Lee, 2026]:

Economic Burden of Asthma (EBA) study (Chen et al. 2013) — a prospective, representative observational study of 618 participants aged 1-85 years (74% aged 18 or older) with self-reported, physician-diagnosed asthma from BC, in which asthma control and the number of exacerbations were measured every 3 months over a year. EBA gives us the overall mean annual exacerbation rate for a person with asthma, \(r = 0.347\), and the overall proportion of time the EBA cohort as a whole spent in each control level over the study period: \(\text{prop}_{\text{wc}} = 0.340\) (well-controlled), \(\text{prop}_{\text{pc}} = 0.474\) (partially-controlled), and \(\text{prop}_{\text{uc}} = 0.186\) (uncontrolled).
GOAL Study (Bateman et al. 2004) — a one-year, randomized, double-blind clinical trial of 3421 participants aged 12-80 years with uncontrolled asthma at study entry, with asthma exacerbations as the primary outcome. An analysis of the GOAL data gives rounded annual exacerbation rates for each control level: well-controlled = 0.1, partially-controlled = 0.2, uncontrolled = 0.3 — i.e. the partially-controlled rate is twice the well-controlled rate, and the uncontrolled rate is three times the well-controlled rate.

Combining these two sources was necessary because the EBA cohort, while representative, did not have enough exacerbation events on its own to robustly estimate a rate for each control level. Instead, we take only the relative rates from GOAL, and solve for an absolute well-controlled rate \(r_{\text{wc}}\) such that the population-weighted average across the three control levels — using the EBA time-in-control proportions together with the GOAL rate ratios — equals the EBA overall rate \(r\):

\[r = \text{prop}_{\text{wc}} \cdot r_{\text{wc}} + \text{prop}_{\text{pc}} \cdot (2 r_{\text{wc}}) + \text{prop}_{\text{uc}} \cdot (3 r_{\text{wc}}) \quad \Longrightarrow \quad r_{\text{wc}} = \dfrac{r}{\text{prop}_{\text{wc}} + 2 \cdot \text{prop}_{\text{pc}} + 3 \cdot \text{prop}_{\text{uc}}}\]

The partially-controlled and uncontrolled rates follow directly from the GOAL ratios: \(r_{\text{pc}} = 2 r_{\text{wc}}\) and \(r_{\text{uc}} = 3 r_{\text{wc}}\). Taking the natural log of each rate gives the \(\beta_k\) values:

\[\begin{split}\beta_1 &:= \ln(r_{\text{wc}}) = \ln(0.1880058) \\ \beta_2 &:= \ln(r_{\text{pc}}) = \ln(0.3760116) \\ \beta_3 &:= \ln(r_{\text{uc}}) = \ln(0.5640174)\end{split}\]

The number of exacerbations predicted by the model is then:

\[\begin{split}N_{\text{asthma}}(a, s, t) &= N(a, s, t) \cdot \eta_{\text{prev}}(a, s, t) \\ N_{\text{exac}}^{\text{(pred)}}(a, s, t) &= \lambda_C(a, s) \cdot N_{\text{asthma}}(a, s, t) \\\end{split}\]

\(N_{\text{asthma}}(a, s, t)\): the number of people of age \(a\), sex \(s\), at timepoint \(t\) with asthma
\(N(a, s, t)\): the number of people of age \(a\), sex \(s\), at timepoint \(t\), from the Population Data described above
\(\eta_{\text{prev}}(a, s, t)\): the prevalence of asthma for age \(a\), sex \(s\), at timepoint \(t\), from Occurrence Model 1: Crude Occurrence
\(\lambda_C(a, s)\): as defined above

and number of hospitalizations is:

\[N_{\text{hosp}}^{\text{(pred)}}(a, s, t) = N_{\text{exac}}^{\text{(pred)}}(a, s, t) \cdot P(\text{hosp})\]

\(N_{\text{exac}}^{\text{(pred)}}(a, s, t)\): the predicted number of exacerbations (of any severity) for age \(a\), sex \(s\), at timepoint \(t\)
\(P(\text{hosp})\): the probability of hospitalization due to asthma given the patient has an asthma exacerbation (a single constant, not stratified by \(a\), \(s\), or \(t\))

As described in the Asthma Exacerbation Severity Model, exacerbations are classified into four severity levels, where very severe is defined as requiring hospital admission. We treat \(P(\text{hosp})\) as the proportion of all exacerbations that are very severe, taken from the Symbicort Given as Needed in Mild Asthma II (SYGMA II) study (Bateman et al. 2018), a double-blind, multi-centre clinical trial of 4176 individuals with mild asthma. SYGMA II reports the distribution of exacerbation severity as 49.5% mild, 19.5% moderate, 28.3% severe, and 2.6% very severe, giving \(P(\text{hosp}) = 0.026\) [Lee, 2026].

Finally, \(\alpha\) can be computed as the ratio of the observed to the predicted number of hospitalizations, for a given age \(a\), sex \(s\), and timepoint \(t\):

\[\alpha(a, s, t) = \dfrac{N_{\text{hosp}}(a, s, t)}{N_{\text{hosp}}^{\text{(pred)}}(a, s, t)}\]

where \(N_{\text{hosp}}(a, s, t)\) is the observed number of hospitalizations, determined from the observed hospitalization rate in CIHI as described in tab1_rate.csv Columns above.

Math: Why α (from Hospitalizations) Applies to λ (All Severities)

Although \(\alpha\) is computed from hospitalizations alone, it is applied to \(\lambda^{(i)}\), the rate of exacerbations of any severity — this is valid because \(P(\text{hosp})\) is treated as a fixed constant, independent of age, sex, province, and timepoint. Substituting \(N_{\text{hosp}}^{\text{(pred)}} = N_{\text{exac}}^{\text{(pred)}} \cdot P(\text{hosp})\), and writing the observed hospitalizations as the true total number of exacerbations times that same constant, \(N_{\text{hosp}} = N_{\text{exac}}^{\text{(true)}} \cdot P(\text{hosp})\), the \(P(\text{hosp})\) terms cancel:

\[\alpha = \dfrac{N_{\text{hosp}}}{N_{\text{hosp}}^{\text{(pred)}}} = \dfrac{N_{\text{exac}}^{\text{(true)}} \cdot P(\text{hosp})} {N_{\text{exac}}^{\text{(pred)}} \cdot P(\text{hosp})} = \dfrac{N_{\text{exac}}^{\text{(true)}}}{N_{\text{exac}}^{\text{(pred)}}}\]

So \(\alpha\) computed from the hospitalization ratio is algebraically identical to the ratio of true to predicted total exacerbations, provided \(P(\text{hosp})\) is indeed constant by age, sex, province, and timepoint (otherwise variation in \(P(\text{hosp})\) would be attributed to \(\alpha\)). Hospitalizations are used to compute it — rather than total exacerbations directly — because they are captured completely in CIHI’s national data, stratified by age/sex/province/year.

\(\alpha\) is computed once per province, age, sex, and timepoint as part of data generation, and saved as: processed_data/{time_delta_tag}/exacerbation_calibration.csv. This file is looked up at simulation runtime — by province, timepoint, sex, and age — to calibrate each agent’s \(\lambda^{(i)}\).

Processed Data¶

Column	Type	Description
`timepoint`	`datetime`	The start of the time interval, e.g. 2024-01-01
`sex`	`str`	`F` = female, `M` = male
`age`	`int`	the age of the patient in years
`province`	`str`	the 2-letter province abbreviation (currently only `BC` and `CA`, see Province Coverage above)
`calibrator_multiplier`	`float`	the calibration multiplier \(\alpha(a, s, t)\) for the given age, sex, timepoint, and province