Asthma Occurrence Model¶

This document describes the asthma occurrence model, which is used to predict the incidence and prevalence of asthma in British Columbia. The model is divided into two parts:

1. Occurrence Model 1: Crude Occurrence: A Generalized Linear Model (GLM) that predicts asthma incidence and prevalence based on age, sex, and timepoint.

2. Occurrence Model 2: Risk Factors: A model that incorporates risk factors such as family history and antibiotic use during infancy to predict asthma incidence and prevalence, along with the results from the first model.

Occurrence Model 1: Crude Occurrence¶

In the first model, we will use data collected from the BC Ministry of Health on the incidence and prevalence of asthma in British Columbia. We will use this data to fit a Generalized Linear Model (GLM) to predict the incidence and prevalence of asthma based on the age, sex, and timepoint. However, asthma occurrence doesn’t just depend on someone’s age or sex, but it also depends on risk factors such as family history and antibiotic use during infancy. We will address these in the second model: Occurrence Model 2: Risk Factors.

Datasets¶

The BC Ministry of Health Administrative Dataset contains asthma incidence and prevalence data for the years 2000-2019, in 5-year age intervals.

The data is formatted as follows:

Model: Generalized Linear Model - Poisson¶

Since our model projects into the future, we would like to be able to extend this data beyond 2019. Our model also makes predictions at customized age intervals, not 5-year age intervals. To obtain these projections, we use a Generalized Linear Model (GLM) with a Poisson distribution and log link function. Incidence and prevalence are counts of people diagnosed with or living with asthma in a given time interval, making the Poisson distribution a natural choice. See Statistical Background for more information on GLMs, including the Poisson distribution and log link function.

Formula¶

Now that we have our distribution and link function, we need to decide on a formula for the linear predictor. The log link means we are modelling \(\log(\bar{p}_{\text{inc}})\) for incidence and \(\log(\bar{p}_{\text{prev}})\) for prevalence. The scaling of the data (rates per 100 people) is absorbed into the intercept \(\beta_0\). We are permitted to use linear combinations of functions of the features in our dataset.

Let’s start with incidence. We want a formula using age, sex, and timepoint. Since asthma depends on factors such as pollution and antibiotic use, and these factors change over time, it follows that asthma incidence should depend on the timepoint. Antibiotic use also depends on age, so we should include age in our formula. Finally, there is a sex difference in asthma incidence, so we should include sex in our formula.

where:

And \(a_i\) is the age, \(t_i\) is the timepoint, \(s_i\) is the sex of individual \(i\).

There are \(4 + 5 + 5 = 14\) coefficients in the incidence model.

Next we have the prevalence. We again want a formula using age, sex, and timepoint. Since asthma prevalence depends on the number of people who have asthma, and this number changes over time, we should include timepoint in our formula. Asthma prevalence also depends on age, so we should include age in our formula. Finally, there is a sex difference in asthma incidence and hence prevalence, so we should include sex in our formula.

where:

and \(a_i\) is the age, \(t_i\) is the timepoint, \(s_i\) is the sex of individual \(i\).

There are \((1 + 1 + 5 + 5) + (2 + 2 + 10 + 10) = 36\) coefficients in the prevalence model.

Assumptions¶

The following assumptions are applied when generating predictions from the fitted model:

• Age: The training data contains 5-year age bands up to a maximum midpoint of 63 years. For ages greater than 63, incidence and prevalence are assumed to remain constant at the rates predicted for age 63.

• Timepoint: Incidence and prevalence trends are predicted for all timepoints using the fitted GLM. At runtime, the timepoint is capped at 2026-01-01 — the maximum timepoint in asthma_occurrence_correction.csv — so that incidence and prevalence remain constant at the 2026 rates for all subsequent timepoints.

Processed Data¶

The processed data produced by this model is stored in asthma_occurrence_predictions.csv (under the time_delta_<days> directory, where <days> is the number of days in the simulation’s time step — e.g. time_delta_365 for yearly intervals). The data contains predicted asthma incidence and prevalence at 1-year age intervals, for each timepoint and sex. The variables are:

Occurrence Model 2: Risk Factors¶

Model 1 produces age-, sex-, and timepoint-specific rates for the general population, but it treats everyone in a stratum identically. In reality, individuals differ in ways that affect their asthma risk — most notably whether a parent has asthma, and whether they received antibiotics in early life. Model 2 builds on Model 1 by incorporating these risk factors, so that the simulation can assign each agent an individualised probability of asthma incidence or prevalence rather than a population average.

This is done in two phases:

• Data generation (run once, before simulation): for every combination of age, sex, and timepoint, a calibration term \(\alpha\) is computed that ensures the population-weighted average of the risk-factor-adjusted probabilities still matches the target rates \(\bar{p}_{\text{prev}}\) and \(\bar{p}_{\text{inc}}\) from Model 1. The results are saved to asthma_occurrence_correction.csv.

• Simulation (at runtime): each agent’s individual risk factors are combined with \(\alpha\) from the lookup table to produce a personalised asthma probability on every simulated timepoint of life.

Within the data generation phase, there is a fixed order of operations across three steps. Computing \(\alpha\) requires knowing \(\log(\omega_{\text{fhx}})\) and \(\log(\omega_{\text{abx}})\) for each risk factor combination \(\lambda\) at each (age, sex, timepoint) stratum, since these determine the individual-level probabilities that must be population-weighted to match the Model 1 targets \(\bar{p}_{\text{prev}}\) and \(\bar{p}_{\text{inc}}\). The log-ORs depend on age-dependent slope parameters \(\beta_{\lambda,\text{age}}\), which are known from the literature for prevalence but must be estimated for incidence — which creates the following sequence:

1. Prevalence calibration. The age-dependent OR slope coefficients (\(\beta_{\text{fhx,age}}\), \(\beta_{\text{abx,age}}\)) for prevalence are derived directly from the literature, so \(\log(\omega_{\text{fhx}})\) and \(\log(\omega_{\text{abx}})\) are fully determined at each age. Broyden-Fletcher-Goldfarb-Shanno (BFGS) then solves for \(\alpha_{\text{prev}}\) per stratum as the value that makes the population-weighted average of the individual prevalence probabilities match the Model 1 target \(\bar{p}_{\text{prev}}\).

2. Incidence \(\beta_{\lambda,\text{age}}\) estimation. Because no published studies provide age-dependent OR slopes for incidence, these are estimated by optimisation under two constraining assumptions: the intercept \(\beta_0\) and antibiotic dose coefficient \(\beta_{\text{abx,dose}}\) in the incidence OR equations are inherited directly from the prevalence values and held fixed — only the age slopes \(\beta_{\text{fhx,age}}\) and \(\beta_{\text{abx,age}}\) are free to vary; and the prevalence and incidence OR equations are equal at age 3, which is guaranteed by the shared intercepts. The optimiser is initialised from the prevalence slope values and finds the age slopes that simultaneously satisfy two conditions: (i) the average incidence across risk factor combinations \(\lambda\), weighted by their population proportions \(\text{prop}(\lambda)\) within each (age, sex, timepoint) stratum, matches \(\bar{p}_{\text{inc}}\) from Model 1; and (ii) the ORs implied by the contingency tables simulated forward one timepoint from the calibrated prevalence distribution at age \(t-1\) match the literature-derived prevalence ORs across age groups. The converged slopes are saved to occurrence_calibration_parameters.json.

3. Incidence calibration. With the estimated \(\beta_{\lambda,\text{age}}\) slopes from step 2, \(\log(\omega_{\text{fhx}})\) and \(\log(\omega_{\text{abx}})\) for incidence are fully determined. BFGS then solves for \(\alpha_{\text{inc}}\) per stratum: the value that makes the population-weighted average of the at-risk incidence probabilities match the Model 1 target \(\bar{p}_{\text{inc}}\).

The diagram below summarises how the components relate, from the Model 1 targets through the offline calibration steps to the personalised probabilities used at runtime.

flowchart TD
    M1["<b>Model 1 population targets</b><br/>$$\overline{\vphantom{b}p}_{\text{prev}} \text{ and } \overline{\vphantom{b}p}_{\text{inc}}$$"]
    LIT["<b>Literature ORs (prevalence)</b>"]

    subgraph OFFLINE["Data generation (run once)"]
        direction TB
        S1["<b>Step 1 — Prevalence calibration</b><br/>
        $$\sum_\lambda \text{prop}(\lambda) \cdot p_{\text{prev},\lambda} = \overline{\vphantom{b}p}_{\text{prev}}$$"]
        S2["<b>Step 2 — Incidence β_age estimation</b><br/>optimise $$\quad \beta_{fhx_{age}}, \beta_{abx_{age}}$$"]
        S3["<b>Step 3 — Incidence calibration</b><br/>
        $$\sum_\lambda \text{prop}(\lambda) \cdot p_{\text{inc},\lambda} = \overline{\vphantom{b}p}_{\text{inc}}$$"]
        S1 -->|"calibrated prevalence at age $$~ t-1$$"| S2
        S2 -->|"$$\beta_{fhx_{age}}, \beta_{abx_{age}}$$"| S3
    end

    SIM["<b>Simulation (runtime)</b><br/>
    $$\text{logit}(p) = \text{logit}(\overline{\vphantom{b}p}) + \log(\omega_{\text{fhx}}) + \log(\omega_{\text{abx}}) - \alpha$$
    "]

    M1 -->|"$$\overline{\vphantom{b}p}_{\text{prev}}$$"| S1
    LIT --> S1
    M1 -->|"$$\overline{\vphantom{b}p}_{\text{inc}}$$ target + literature prevalence ORs"| S2
    M1 -->|"$$\overline{\vphantom{b}p}_{\text{inc}}$$"| S3
    S1 -->|"$$\alpha_{\text{prev}}$$"| SIM
    S3 -->|"$$\alpha_{\text{inc}}$$"| SIM
    M1 -->|"$$\overline{\vphantom{b}p}_{\text{prev}}, \overline{\vphantom{b}p}_{\text{inc}}$$"| SIM

    classDef target fill:#e3f2fd,stroke:#1565c0,color:#0d2b45;
    classDef lit fill:#fff3e0,stroke:#e65100,color:#3a2400;
    classDef sim fill:#e8f5e9,stroke:#2e7d32,color:#10300f;
    class M1 target;
    class LIT lit;
    class SIM sim;

Model: Risk Factors¶

We want to incorporate the effects of family history and antibiotic use on asthma incidence and prevalence.

Formula¶

Before we begin, let us define some terms. We have two risk factors we are interested in: family history and antibiotic use. There are \(2 * 4 = 8\) possible combinations of these two risk factors:

The effect of each risk factor on asthma risk is expressed as an odds ratio (OR) sourced from external published studies. An odds ratio of \(\omega\) means that a person with the risk factor has \(\omega\) times the odds of having asthma compared to a person without it.

We work in log-odds (logit) space rather than probability space because log-odds are additive: the combined effect of independent risk factors is simply the sum of their individual log-odds contributions. Because independent ORs are multiplicative, their logarithms are additive:

Prevalence¶

Applying individual risk factor ORs directly to the Model 1 log-odds would shift the population-weighted average probability away from the target \(\bar{p}_{\text{prev}}\). The calibration term \(\alpha\) corrects for this: it is a single scalar per (age, sex, timepoint) stratum that shifts the baseline log-odds so that the population-weighted average of \(p_{\text{prev}}\) across all 8 risk factor combinations matches \(\bar{p}_{\text{prev}}\). It plays the same role as an intercept correction in a regression model.

The predicted prevalence for an individual agent is:

Solving for the Correction Term¶

The correction term \(\alpha\) is solved once, prior to runtime, using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm to find the value that minimises the difference between the population-weighted average asthma prevalence and the Model 1 target asthma prevalence:

where \(\text{prop}(\lambda)\) is the proportion of the population with risk factor combination \(\lambda\).

Age-dependent Odds Ratios for Prevalence¶

The log-ORs for each risk factor are age-dependent — their functional forms and fixed coefficients are derived from external published studies.

Antibiotic Exposure¶

The antibiotic log-OR was derived by Lee et al. [Lee, 2024a] using a random effects meta-regression model. It applies only to children aged 7 or under who received antibiotics in infancy, and is zero otherwise:

where \(a_i\) is the agent’s age and \(d_i\) is the number of courses of antibiotics taken during the first year of life (capped at 3). All three coefficients are sourced directly from the Lee et al. meta-regression, which estimated age as an explicit covariate:

• \(\beta_{\text{abx}_0} = 1.826\)

• \(\beta_{\text{abx}_\text{age}} = 0.225\)

• \(\beta_{\text{abx}_\text{dose}} = 0.053\)

Family History¶

The family history log-OR was derived from the CHILD Study by Patrick et al. [Patrick, 2020] using logistic regression. See Family History Model for how the family history risk factor is assigned to agents and how \(p = 0.2927\) is estimated from the same study. The log-OR applies to all ages but the age-dependent component plateaus at age 5:

where \(f_i = 1\) if at least one parent has asthma, 0 otherwise. Both coefficients are derived from the two empirical ORs reported by Patrick et al. at ages 3 and 5 — OR = 1.13 at age 3 and OR = 2.40 at age 5. The age-dependent slope is the linear interpolation between those two points on the log-OR scale:

• \(\beta_{\text{fhx}_0} = \log(1.13) = 0.122\)

• \(\beta_{\text{fhx}_\text{age}} = \dfrac{\log(2.40) - \log(1.13)}{2} = 0.377\)

For ages above 5, the OR is held constant at the age-5 value.

For prevalence, all OR coefficients — including the age-dependent slopes — are fully determined by the literature before calibration runs. The only quantity calibrated for prevalence is the scalar correction term \(\alpha\), which shifts the population-weighted average onto the Model 1 target \(\bar{p}_{\text{prev}}\).

Incidence¶

The incidence formula has the same structure as prevalence, but applies only to the at-risk population — agents who do not currently have asthma. The predicted incidence for an individual agent is:

Solving for the Correction Term¶

The correction term \(\alpha\) is solved once, prior to runtime, using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm to find the value that minimises the difference between the population-weighted average asthma incidence and the Model 1 target asthma incidence:

where \(\text{prop}_{\text{no asthma},\lambda}(t-1)\) is the proportion of the population who are asthma-free at \(t-1\) and have risk factor combination \(\lambda\). Only asthma-free agents are included because incidence counts new diagnoses only.

Calibrating Age-Dependent Odds Ratios¶

Unlike prevalence, there are no published studies that directly estimate the age-dependent slopes of \(\log(\omega_{\text{fhx}})\) and \(\log(\omega_{\text{abx}})\) for incidence. Instead, the slopes are estimated by optimisation under two constraining assumptions:

• The intercept \(\beta_0\) and antibiotic dose coefficient \(\beta_{\text{abx,dose}}\) in the incidence OR equations are inherited from prevalence and held fixed. Only the age slopes \(\beta_{\text{fhx,age}}\) and \(\beta_{\text{abx,age}}\) are free to vary.

• The prevalence and incidence OR equations are equal at age 3, which is guaranteed by sharing the same intercept.

The optimiser is initialised from the corresponding prevalence slope values and must simultaneously satisfy two conditions:

1. The average incidence across risk factor combinations \(\lambda\), weighted by their population proportions \(\text{prop}(\lambda)\) within each (age, sex, timepoint) stratum, matches the Model 1 target \(\bar{p}_{\text{inc}}\).

2. The ORs implied by the contingency tables — simulated forward one timepoint from the calibrated prevalence distribution at age \(t-1\) — match the literature-derived prevalence ORs across age groups.

Condition 2 ensures that introducing individual-level incidence risk does not distort the aggregate OR structure established by the prevalence literature. See Contingency Tables for an introduction to contingency tables and worked examples.

To evaluate condition 2, we track a cohort across one timepoint using contingency tables of risk factor combination \(\lambda\) against asthma diagnosis. The labels \(t=0\) and \(t=1\) denote two consecutive timepoints. The procedure is:

1. Build the baseline table at \(t=0\) so that its odds ratio equals the (literature-anchored) prevalence OR at the previous timepoint.

2. Step forward one timepoint by applying reassessment to those who already had asthma and new incidence — computed with the candidate age-slopes \(\beta_{\lambda,\text{age}}\) — to those who did not, producing the combined table at \(t=1\).

3. Compare the odds ratio of that combined table (the prevalence OR the model actually produces at age \(a\)) against the target prevalence OR, and adjust \(\beta_{\lambda,\text{age}}\) until they match.

The four tables below build this up step by step, and the diagram summarises the flow.

flowchart LR
    BASE["<b>Baseline table (t=0)</b><br/>built from literature-based<br/>prevalence OR"]
    NEW["<b>New diagnoses (t=1)</b><br/>incidence with candidate $$~\beta_{age}$$"]
    EX["<b>Existing diagnoses (t=1)</b><br/>reassessment with prob $$~\rho$$"]
    COMB["<b>Combined table (t=1)</b><br/>prevalence at age a:<br/> evaluate prevalence OR vs target prevalence OR"]

    BASE -->|"no asthma at t=0"| NEW
    BASE -->|"had asthma at t=0"| EX
    NEW --> COMB
    EX --> COMB
    COMB -.->|"adjust $$~\beta_{age}$$"| NEW

    classDef base fill:#e3f2fd,stroke:#1565c0,color:#0d2b45;
    classDef mid fill:#fff3e0,stroke:#e65100,color:#3a2400;
    classDef comb fill:#e8f5e9,stroke:#2e7d32,color:#10300f;
    class BASE base;
    class NEW,EX mid;
    class COMB comb;

Baseline Contingency Table (t=0)¶

Since the table is 2×2, we compare one risk factor combination at a time against the no-risk-factor baseline (\(\lambda = 0\), i.e. no family history and no antibiotic exposure). This produces 7 separate tables — one for each non-baseline combination. The non-binary nature of antibiotic dose is handled implicitly through the \(\lambda\) indexing: dose levels 1, 2, and 3 each appear as distinct combinations and are each compared independently against the baseline rather than against each other.

For each comparison, let \(N\) be the total population across all risk factor combinations (a hypothetical size, e.g. 100,000 — the odds ratio is scale-invariant). The row and column totals are then:

where \(n_1\) is the number of people with risk factor \(\lambda\), \(n_0\) is the number with no risk factors, \(n\) is the grand total of the 2×2 table, and \(n_2\) is the total number with asthma across both groups.

Given \(n_1\), \(n_0\), \(n_2\), \(n\), and the odds ratio \(\omega_\lambda\), we solve for the cell count \(a_0\) (people with both risk factor \(\lambda\) and asthma) such that the implied odds ratio of the table matches \(\omega_\lambda\). This is a non-trivial solve because all four cells are simultaneously constrained by the marginal totals and the odds ratio — we use the method from Di Pietrantonj (2006) [C., 2006]. The remaining cells follow directly from \(a_0\) and the marginal totals.

Existing Diagnoses: Reassessment at t=1¶

In our model, an asthma diagnosis is not permanent — a person may be reassessed and lose their diagnosis from one timepoint to the next. This table tracks what happens at \(t=1\) to the \(a_0 + c_0\) people who already had asthma at \(t=0\). Let \(\rho\) be the probability of retaining a diagnosis. Then each cell is simply the corresponding baseline cell scaled by \(\rho\) (retained) or \(1 - \rho\) (lost):

where \(\rho\) is the probability of retaining an asthma diagnosis from one timepoint to the next, applied equally regardless of risk factor status.

New Diagnoses at t=1¶

This table tracks what happens at \(t=1\) to the \(b_0 + d_0\) people who were asthma-free at \(t=0\). Each person may receive a new diagnosis based on their incidence probability: \(p_{\text{inc},\lambda}\) for those with risk factor \(\lambda\), and \(p_{\text{inc},0}\) for those with no risk factors.

Combined Contingency Table: Existing and New Diagnoses (t=1)¶

The combined table at \(t=1\) sums the reassessment and new diagnosis components for each cell, giving the full joint distribution of risk factor status and asthma diagnosis at the end of the time interval:

where:

From these values, we can compute the odds ratio:

Optimization¶

We want to find the age-dependent slope values (\(\beta_{\lambda, \text{age}}\)) that minimize the mean absolute difference between \(\log(\hat{\omega})\) and \(\log(\omega_{\lambda})\) — the fixed, age-dependent log-ORs for prevalence sourced from external studies — averaged across all age groups and all 7 non-baseline risk factor combinations.

Once optimised, these slopes (\(\beta_{\lambda, \text{age}}\)) are used to compute \(\log(\omega_{\text{fhx}})\) and \(\log(\omega_{\text{abx}})\) for each (age, sex, timepoint) stratum. BFGS then uses those stratum-specific log-ORs to solve for the calibration term \(\alpha\) per stratum, which is stored in asthma_occurrence_correction.csv as described in the Processed Data section below.

Processed Data¶

The calibration terms produced by Model 2 are stored in asthma_occurrence_correction.csv (under the time_delta_<days> directory matching the simulation’s time step). Each row gives the value of \(\alpha\) for a specific age, sex, timepoint, and outcome type (prevalence or incidence). This file is used at runtime by the simulation to look up the correction for each agent at each timepoint of life.