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 use the processed birth estimate data described in Birth Model (birth_estimate.csv). That file provides N (total births, both sexes combined) and prop_male (proportion of births that are male) per timepoint and province; we derive the per-sex counts shown below as N * prop_male (male) and N * (1 - prop_male) (female):

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 and \(H\) is the Heaviside step function. \(t^{(i)}\) is the birth timepoint; since the underlying data is only available at yearly granularity, this is simply the agent’s birth year (e.g. 2008).

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 datetime as dt
import itertools
import statsmodels.api as sm
import statsmodels.formula.api as smf
from leap.utils import date_range, TimeDelta
from leap.data_generation.utils import heaviside
from leap.data_generation.antibiotic_data import load_birth_data, convert_sex_to_numeric, \
    convert_timepoint_to_numeric

# 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 = list(date_range(start=dt.datetime(2000, 1, 1), stop=dt.datetime(2018, 12, 31), step=TimeDelta(years=1)))

# 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", "timepoint"]
)
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()

We can now merge the df_abx and df_birth dataframes:

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

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

fig = px.line(
    df,
    x="timepoint",
    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["sex"].apply(convert_sex_to_numeric)
df["timepoint"] = df["timepoint"].apply(convert_timepoint_to_numeric)

formula = "n_abx ~ timepoint + sex + heaviside(timepoint, 2005.0) * timepoint"

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:                   38
Model:                            GLM   Df Residuals:                       33
Model Family:        NegativeBinomial   Df Model:                            4
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -277.33
Date:                Tue, 08 Sep 2026   Deviance:                       27.142
Time:                        23:20:03   Pearson chi2:                     27.1
No. Iterations:                     6   Pseudo R-squ. (CS):              1.000
Covariance Type:            nonrobust                                         
==========================================================================================================
                                             coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------------------------------
Intercept                                 83.0400     12.249      6.779      0.000      59.032     107.048
timepoint                                 -0.0418      0.006     -6.826      0.000      -0.054      -0.030
sex                                        0.2602      0.012     21.859      0.000       0.237       0.283
heaviside(timepoint, 2005.0)              87.6168     12.855      6.816      0.000      62.422     112.812
heaviside(timepoint, 2005.0):timepoint    -0.0437      0.006     -6.808      0.000      -0.056      -0.031
==========================================================================================================

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

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

Example: Sex = Male, Year = 2008¶

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

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.678742

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.999152, 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 {timepoint}: {n_abx}")

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

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 = {timepoint}",
    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: int, timepoint: float, β_0: float, β_s: float, β_t: float, β_h: float, β_th: float, θ: float
):
    η = (
        β_0 + β_s * sex + β_t * timepoint +
        β_h * heaviside(timepoint, 2005) +
        β_th * heaviside(timepoint, 2005) * timepoint
    )
    µ = 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 = list(date_range(dt.datetime(2000, 1, 1), dt.datetime(2018, 1, 1), step=TimeDelta(years=2)))
df_pred = pd.DataFrame(
    data=list(itertools.product(sexes, years)),
    columns=["sex", "timepoint"]
)
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"], convert_timepoint_to_numeric(x["timepoint"]), β_0, β_s, β_t, β_h, β_th, θ),
    axis=1
)
df_pred["r"] = r
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", "timepoint", "nabx"]
)
df_dist = pd.merge(df_dist, df_pred, on=["timepoint", "sex"])
df_dist.head()

grouped_df = df_dist.groupby(["sex", "timepoint"])
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="timepoint",
    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("timepoint=", "Timepoint: ")))
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)