Antibiotic Exposure Model¶

In this section, we will describe the model used to predict the number of antibiotics prescribed to infants in their first year of life. This model will be incorporated into the risk factors for developing asthma later in life.

Input Datasets¶

We obtained data from the BC Ministry of Health on the number of antibiotics prescribed to infants in their first year of life. The data is available for the years 2000 to 2018. The data is formatted as follows:

Since the n_abx column gives us the total number of antibiotics prescribed, we need to use population data to convert this to a per infant value. We obtained population data from Table 17-10-00005-01 from Statistics Canada, the same past-data source used in the Birth Model:

Merged Dataset: InfantAbxBC.csv¶

The two input datasets are merged on year and sex to produce InfantAbxBC.csv, saved in the processed_data directory. This combined dataset is what the GLM is fitted on:

Model: Generalized Linear Model - Negative Binomial¶

Since our model projects into the future, we would like to be able to extend this data beyond 2018. To obtain these projections, we use a Generalized Linear Model (GLM). A GLM is a type of regression analysis which is a generalized form of linear regression.

The response variable — number of antibiotic courses prescribed during the first year of life — is overdispersed count data, so we use the Negative Binomial distribution with a log link function. See Example 3: Negative Binomial Distribution with Log Link in Statistical Background for the full distributional derivation and motivation.

The mean parameter has a lower bound to prevent unrealistically small extrapolations into future years:

In other words, the predicted mean number of antibiotic courses per infant is at least 0.05.

Formula¶

The birth data enters the model as a log offset on n_birth, so the GLM predicts a per-capita rate directly rather than a raw count.

Since prescribing practices change over time, and since infections requiring antibiotic prescriptions also change over time, birth year is included in the formula. Sex is also included, since there are known sex differences in antibiotic prescriptions.

There is an additional factor specific to BC regulations. In 2005, the BC government introduced an antibiotic conservation program, which reduced the number of antibiotics prescribed [Mamun, 2019]. To account for this structural break, the formula includes a Heaviside step function \(H\), which is 0 for years before 2005 and 1 for years from 2005 onward.

where:

And \(s^{(i)}\) is the sex, \(t^{(i)}\) is the birth year, and \(H\) is the Heaviside step function.

Usage in Simulation¶

Once fitted, the \(\beta\) coefficients and \(\theta\) (the Negative Binomial overdispersion parameter; see Example 3: Negative Binomial Distribution with Log Link) are stored in config.json and used at runtime. When an agent is initialised at birth, the simulation draws their antibiotic exposure count directly from the Negative Binomial distribution. For an agent with sex \(s^{(i)}\) and birth year \(t^{(i)}\):

1. Compute the linear predictor \(\eta^{(i)}\) using the formula above.

2. Convert to the mean: \(\mu^{(i)} = \max(\exp(\eta^{(i)}),\; 0.05)\).

3. Convert to the Negative Binomial success probability: \(p^{(i)} = \theta / (\theta + \mu^{(i)})\).

4. Draw the exposure count: \(n_{\text{abx}}^{(i)} \sim \text{NegBin}(\theta,\, p^{(i)})\).

This count is fixed for the agent’s lifetime and is capped at 3 courses when computing the antibiotic exposure odds ratio \(\omega_{\text{abx}}\) in the Asthma Occurrence Model. See Simulation for how antibiotic exposure is assigned during agent initialisation.

The worked example below demonstrates steps 1–4 for a specific agent (sex = male, birth year = 2008), and then visualises the resulting distribution across multiple birth years and sexes.

Example¶

Setup¶

%load_ext autoreload
%autoreload 2

import plotly.io as pio
pio.renderers.default = "sphinx_gallery"
import plotly.express as px
import pandas as pd
import numpy as np
import itertools
import statsmodels.api as sm
import statsmodels.formula.api as smf
from leap.data_generation.utils import heaviside
from leap.data_generation.antibiotic_data import load_birth_data

# Plot Settings
pio.templates.default = "plotly_white"
config = {
    'toImageButtonOptions': {
        'format': 'png',
        'scale': 2
    }
}

Data¶

Since the antibiotic dosing data is private, we will create here a mock dataset. In this dataset, we will have the total number of courses of antibiotics prescribed to all infants in BC, stratified by sex and year.

# Data was collected for the years 2000-2018
years = np.arange(2000, 2018)

# Mock data for antibiotic courses prescribed to infants using an approximate linear trend
n_abx_female = [-640 * i + 16000 for i, _ in enumerate(years)]
n_abx_male = [-880 * i + 22000 for i, _ in enumerate(years)]
n_abx = n_abx_female + n_abx_male

# Create a DataFrame for antibiotic courses prescribed to infants
df_abx = pd.DataFrame(
    data=list(itertools.product(["F", "M"], years)),
    columns=["sex", "year"]
)
df_abx["n_abx"] = n_abx
df_abx.head()

Next, since the antibiotic data measures the number of antibiotics prescribed at a population level, we need to use population data from Statistics Canada to determine the number of antibiotics prescribed per capita:

df_birth = load_birth_data()
df_birth.head()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/pandas/core/indexes/base.py:3641, in Index.get_loc(self, key)
   3640 try:
-> 3641     return self._engine.get_loc(casted_key)
   3642 except KeyError as err:

File pandas/_libs/index.pyx:168, in pandas._libs.index.IndexEngine.get_loc()
--> 168 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/index.pyx:197, in pandas._libs.index.IndexEngine.get_loc()
--> 197 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/hashtable_class_helper.pxi:7668, in pandas._libs.hashtable.PyObjectHashTable.get_item()
-> 7668 'Could not get source, probably due dynamically evaluated source code.'

File pandas/_libs/hashtable_class_helper.pxi:7676, in pandas._libs.hashtable.PyObjectHashTable.get_item()
-> 7676 'Could not get source, probably due dynamically evaluated source code.'

KeyError: 'year'

The above exception was the direct cause of the following exception:

KeyError                                  Traceback (most recent call last)
Cell In[4], line 1
----> 1 df_birth = load_birth_data()
      2 df_birth.head()

File /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/leap/data_generation/antibiotic_data.py:88, in load_birth_data(province, min_year, max_year)
     81 df.rename(
     82     columns={"REF_DATE": "year", "GEO": "province", "SEX": "sex", "VALUE": "n_birth"},
     83     inplace=True
     84 )
     86 # select only the age = 0 age group and the years where min_year <= year <= max_year
     87 df = df.loc[
---> 88     (df["year"] >= min_year) & 
     89     (df["year"] <= max_year) & 
     90     (df["AGE_GROUP"] == "0 years")
     91 ]
     93 # select only the columns we need
     94 df = df[["year", "province", "sex", "n_birth"]]

File /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/pandas/core/frame.py:4378, in DataFrame.__getitem__(self, key)
   4374 
   4375         if is_single_key:
   4376             if self.columns.nlevels > 1:
   4377                 return self._getitem_multilevel(key)
-> 4378             indexer = self.columns.get_loc(key)
   4379             if is_integer(indexer):
   4380                 indexer = [indexer]
   4381         else:

File /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/pandas/core/indexes/base.py:3648, in Index.get_loc(self, key)
   3643     if isinstance(casted_key, slice) or (
   3644         isinstance(casted_key, abc.Iterable)
   3645         and any(isinstance(x, slice) for x in casted_key)
   3646     ):
   3647         raise InvalidIndexError(key) from err
-> 3648     raise KeyError(key) from err
   3649 except TypeError:
   3650     # If we have a listlike key, _check_indexing_error will raise
   3651     #  InvalidIndexError. Otherwise we fall through and re-raise
   3652     #  the TypeError.
   3653     self._check_indexing_error(key)

KeyError: 'year'

We can now merge the df_abx and df_birth dataframes:

df = pd.merge(df_abx, df_birth, on=["year", "sex"], how="left")
df.head()

color_map={
    "M": "#09bfc4",
    "F": "#f39c12",
}

fig = px.line(
    df,
    x="year",
    y="n_abx",
    facet_col="sex",
    color="sex",
    labels={"n_abx": "Total Number of Antibiotic Courses Prescribed to Infant Population"},
    color_discrete_map=color_map
)
fig.update_xaxes(title_text="Birth Year", nticks=10)
fig.update_yaxes(title_font=dict(size=12))

fig.show()

Model¶

df["sex"] = df.apply(lambda x: 1 if x["sex"] == "F" else 2, axis=1)

formula = "n_abx ~ year + sex + heaviside(year, 2005) * year"

Since we are using a simulated dataset, we set $\theta = 800$.

θ = 800
alpha = 1 / θ

# Fit the GLM model
model = smf.glm(
    formula=formula,
    data=df,
    family=sm.families.NegativeBinomial(alpha=alpha),
    offset=np.log(df["n_birth"])
)
results = model.fit(maxiter=1000)

print(results.summary())

                 Generalized Linear Model Regression Results                  
==============================================================================
Dep. Variable:                  n_abx   No. Observations:                   36
Model:                            GLM   Df Residuals:                       31
Model Family:        NegativeBinomial   Df Model:                            4
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -260.35
Date:                Tue, 12 Aug 2025   Deviance:                       18.039
Time:                        13:53:48   Pearson chi2:                     18.0
No. Iterations:                     5   Pseudo R-squ. (CS):              1.000
Covariance Type:            nonrobust                                         
==============================================================================================
                                 coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------------------
Intercept                     74.9737     16.193      4.630      0.000      43.236     106.711
year                          -0.0377      0.008     -4.664      0.000      -0.054      -0.022
sex                            0.2595      0.012     21.258      0.000       0.236       0.283
heaviside(year, 2005)         87.6702     16.652      5.265      0.000      55.033     120.308
heaviside(year, 2005):year    -0.0437      0.008     -5.259      0.000      -0.060      -0.027
==============================================================================================

Now that the model has been trained, we can obtain the $\beta$ parameters:

β_0 = results.params["Intercept"]
β_t = results.params["year"]
β_s = results.params["sex"]
β_h = results.params["heaviside(year, 2005)"]
β_th = results.params["heaviside(year, 2005):year"]

Example: Sex = Male, Year = 2008¶

sex = 2
year = 2008
η = β_0 + β_s * sex + β_t * year + β_h * heaviside(year, 2005) + β_th * heaviside(year, 2005) * year

Now, $\eta$ is the number of antibiotics prescribed per number of infants; what we want to know, is given an age, birth year, and sex, how many antibiotics were prescribed to that person during infancy? To do so, we need to use the Negative Binomial distribution. First, we compute $\mu$, using the link function.

µ = np.exp(η)

Next, we make sure that $\mu > 0.05$:

µ = max(µ, 0.05)
print(f"µ = {µ:.6f}")

µ = 0.674037

The standard parametrization of the Negative Binomial function uses $p$ and $r$, not $\mu$. So, we need to convert $\mu$ to $p$:

p = θ / (θ + µ)
r = θ
print(f"p = {p:.6f}, r = {r:.4f}")

p = 0.999158, r = 800.0000

The Negative Binomial distribution gives us the probability of an infant being prescribed $k$ courses of antibiotics during infancy:

$$ \begin{aligned} P(Y = k; r, p) := \binom{k+r-1}{k}(1-p)^k p^r \end{aligned} $$

n_abx = np.random.negative_binomial(r, p)
sex = "female" if sex == 1 else "male"
print(f"Number of courses of antibiotics prescribed during infancy to a {sex} person born in {year}: {n_abx}")

Number of courses of antibiotics prescribed during infancy to a male person born in 2008: 1

The number of antibiotics obtained is just a sampling from the distribution; if we repeated the selection multiple times, we could get different values. To visualize this, we can look at the distribution:

n_abx = np.random.negative_binomial(r, p, size=100000)
fig = px.histogram(
    n_abx,
    nbins=10,
    title=f"Distribution of Antibiotic Courses Prescribed in Infancy<br>Sex = {sex}, Birth Year = {year}",
    color_discrete_sequence=["#09bfc4"],
    width=600
)
fig.update_xaxes(title_text="Number of Courses of Antibiotics", nticks=10)
fig.update_yaxes(title_text="Count")
fig.update_layout(
    legend_title_text="",
    title_x=0.5,
    title_y=0.95,
    margin=dict(t=50), # Adjust top margin for title
    showlegend=False
)
fig.show(config=config)

Example: Visualize Multiple Years / Sexes¶

The example above was just the distribution for a single year and sex; we can visualize multiple years / sexes together:

def compute_abx_distribution_parameters(sex, year, β_0, β_s, β_t, β_h, β_th, θ):
    η = (
        β_0 + β_s * sex + β_t * year +
        β_h * heaviside(year, 2005) +
        β_th * heaviside(year, 2005) * year
    )
    µ = np.exp(η)
    µ = max(µ, 0.05)  # Ensure µ is greater than 0.05
    p = θ / (θ + µ)
    return p
    

Let’s create a dataframe with even years from 2000 to 2018:

sexes = [1, 2]
years = np.arange(2000, 2018, step=2)
df_pred = pd.DataFrame(
    data=list(itertools.product(sexes, years)),
    columns=["sex", "year"],
    dtype='Int64'
)
df_pred.head()

Next, let’s find $p$ and $r$ for each year / sex combination:

df_pred["p"] = df_pred.apply(
    lambda x: compute_abx_distribution_parameters(x["sex"], x["year"], β_0, β_s, β_t, β_h, β_th, θ),
    axis=1
)
df_pred["r"] = df_pred.apply(
    lambda x: r,
    axis=1
)
df_pred.head()

We now create a second dataframe, which we will use to sample the Negative Binomial distribution.

df_dist = pd.DataFrame(
    data=list(itertools.product(sexes, years, np.zeros(1000))),
    columns=["sex", "year", "nabx"],
    dtype=int
)
df_dist = pd.merge(df_dist, df_pred, on=["year", "sex"])
df_dist.head()

grouped_df = df_dist.groupby(["sex", "year"])
df_dist["nabx"] = grouped_df.apply(
    lambda x: np.random.negative_binomial(x["r"], x["p"], size=len(x)), include_groups=False
).explode().astype(int).values
df_dist.head()

Now we can visualize our results:

color_map={
    2: "#09bfc4",
    1: "#f39c12",
}

fig = px.histogram(
    df_dist,
    x="nabx",
    facet_col="year",
    facet_col_wrap=3,
    color="sex",
    nbins=10,
    title=f"Distribution of Antibiotic Courses Prescribed in Infancy",
    height=1000,
    width=850,
    facet_row_spacing=0.1,
    color_discrete_map=color_map,
    barmode="group"
)
fig.for_each_trace(lambda t: t.update(name="female" if t.name == "1" else "male"))
fig.for_each_annotation(lambda a: a.update(text=a.text.replace("year=", "Year: ")))
fig.update_xaxes(
    title_text="Number of Courses of Antibiotics",
    nticks=10,
    row=1,
    title=dict(font=dict(size=12))
)
fig.update_xaxes(showticklabels=True, nticks=10)
fig.update_yaxes(title_text="Count", col=1)
fig.update_layout(
    legend_title_text="",
    title_x=0.5,
    title_y=0.95,
    margin=dict(t=150), # Adjust top margin for title
    showlegend=True
)
fig.show(config=config)