Usage data for 26 European countries for the years 2008–2018 is drawn from IQVIA MIDAS database. IQVIA reports the volume of sales of antibiotic molecules used in human medicine based on national surveys. We convert antibiotic sales into Defined Daily Doses (DDDs) using ATC/DDD Index 2020. The World Health Organization defines DDD as the assumed average maintenance dose per day for a drug used for its main indication in adults. We adjust for the population using population estimates from the World Bank DataBank to obtain DDDs per 1000 inhabitants (DpTI).

Antibiotic molecules are active chemical compounds that can be broadly categorized into classes of antibiotics based on their mechanisms of action with bacteria. Using the ATC/DDD index we aggregate DpTI for antibiotics into 12 antibiotic classes, some of which contain only a single antibiotic (Supplementary Table S1). Sales are thus identified by antibiotic class, year, and country. We categorize the 26 European countries into Northern, Southern, Eastern and Western Europe, as detailed in Supplementary Table S2.

Usage levels vary significantly across classes and countries, as shown in Supplementary Figs. S1–S3. To help make the usage data comparable, we convert DpTI into z-scores. A z-score is defined as the difference between an observation and the sample mean divided by the standard deviation for a class in each country over the years 2008–2018.

The European Antimicrobial Resistance Surveillance System (EARS–Net) collects data on antimicrobial resistance in 8 different bacterial pathogens to 12 antibiotic classes as reported in the Surveillance Atlas of Infectious Diseases. The bacteria-antibiotic combinations covered are detailed in Supplementary Table S3. According to EARS–Net data documentation, all main geographical regions are covered, and, on average, data are considered as representative of the national epidemiology. However, the percentage of population coverage varies among reporting countries and the population under surveillance changes over time.

The data record resistance of specific bacteria—in percentage—by specific antibiotic class, year, and country. We categorize the 8 different bacterial pathogens into Gram-positives and Gram-negatives, as detailed in Supplementary Table S4. Our resistance data therefore vary by country, bacteria-antibiotic class combination, and year from 2008 to 2018, with a range between 0 and 100%, as shown in Supplementary Figs. S4–S7.

We combine the two datasets using year, country, and antibiotic class to obtain resistance and usage data, for 26 European countries, 11 years, and 26 bacteria-class combinations. We construct a panel on the universe of all units, where a unit is a combination given by country, year, and bacteria-class. Since some countries did not report resistance for every bacteria-class in all years, our panel is an unbalanced one with 6586 observations.

Our main dependent variable is the prevalence of resistance, as given by the percentage of resistant isolates identified by EARS–Net. We use usage z-scores or change in z-scores, depending on the specification employed, as our main explanatory variable. We select a 6-year panel of the resistance data covering the years 2012–2017 for our analysis. This choice helps set up the event study design with an event horizon running from 1 year prior to 4 years after usage, so that we require usage data from 2008–2018.

We use three related empirical models to estimate the causal effect of usage on resistance: (1) a Distributed–Lag (D–L) model with fixed effects, (2) a D–L model with first differences and (3) an Event–Study (E–S) model with binned endpoints. Binning refers to the practice of treating the last lag (lead) as an open interval capturing all known changes that have happened (will happen) in the past (future). Our outcome variable is the resistance for a bacteria-class i, in country c, in year t, \(R_}}}}}}},t}\). The explanatory variable of interest in model (1) is the z-score of usage, happening j periods away for a bacteria-class i, in country c, in year t, \(U_}}}}}}},t - j}\). The explanatory variable of interest in models (2) and (3) is the change in z-score of usage, happening j periods away for a bacteria-class i, in country c, in year t, \(\Delta U_}}}}}}},t - j}\).

Models (1) and (2) are specified as:

With fixed effects:

$$R_}}}}}}},t} = \mathop \nolimits_^4 }}}}}}},t - j} + \mu _i + \mu _c + \theta _t + \varepsilon _}}}}}}},t}}$$

(1)

With first differences:

$$\Delta R_}}}}}}},t} = \mathop \nolimits_^4 }U_}}}}}}},t - j} + \theta _t + \varepsilon _}}}}}}},t}}$$

(2)

The event window spans the period from 1 year before usage (or to a change in usage) to 4 years after. We include bacteria-class (µi) and country (µc) fixed effects in Model (1), which accounts for unobserved bacteria, class and country effects which are constant across time. In Model (2), first–differencing resistance and usage controls for these effects. For both models, we include year fixed effects (θt). The error term is given by εi,c,t. The coefficients γj and λj denote the marginal effects of usage on resistance, measuring the slope of these effects from one year to the next. The estimated values of γj and λj in Models (1) and (2) are expected to differ significantly only if the effect of usage on resistance continues to unfold beyond the 4-year window [25].

Model (3) is specified as:

$$R_}}}}}}},t} = \mathop \nolimits_^4 _}}}}}}},t - j} + \mu _i + \mu _c + \theta _t + \varepsilon _}}}}}}},t}}$$

(3)

where the binned variables \(v_}}}}}}},c,t}^\) are defined as:

$$v_}}}}}}},c,t}^ = \left\} \nolimits_^ \quad}}}}}}}\;j = - 2} } \\ \quad\quad\quad\quad}}}}}}} - 2 \, < \, j \, < \, 4} \\ \nolimits_^\infty }}}}}}},t - k}\quad\quad\quad}}}}}}}\;j = 4} } \end} \right.$$

(4)

Model (3) is a regression of levels on binned changes. The coefficients βj are the treatment effects, j time periods before or after usage, dynamically unfolding over time and are expressed relative to a reference period. The coefficient for the reference period is normalized to zero. Binning the upper and lower endpoints is equivalent to assuming that \(\gamma _j = 0\) for all j > 4 and for all \(j \le - 2\). In other words, this model assumes effects of usage on resistance stay constant for all j > 4 and for all \(j \le - 2\). Due to the nature and construction of binned variables, the E-S model requires data 2 years prior to the event. These assumptions on the effect window acknowledge that dynamic effects cannot be estimated using the infinite past and the infinite future due to limited data availability and sample restrictions [25].

Model 3 provides readily interpretable coefficients βj whereas the coefficients γj and λj from Models 1 and 2 must be linearly transformed to derive the dynamic effects. Statistical properties such as consistency and asymptotic normality are preserved during this linear transformation of D–L model estimates, and variances and covariances of the estimated parameters γj and λj can be used to recover standard errors for the estimated βj using a standard linear combination formula [25]. These estimated dynamic treatment effects are unbiased under linear and additive assumptions. Event study models with binned endpoints and D–L models are expected to yield similar treatment effects [25].

The reported dynamic treatment effects estimated using all three models are interpreted as the impact on resistance due to a 1 standard deviation increase in usage. We can observe this impact in the years prior to and after usage. Identification of dynamic treatment effects requires that there be no statistically significant impacts on resistance before any usage happens. Moreover, identification is achieved within bacteria-class and country over time.

We also estimate all the models using segmented data. Specifically, we segment the data by Gram-positive and Gram-negative data to assess whether they respond differently. We also segment by countries, dividing the countries into four regions, West, North, East, and South, to assess whether the same effects are observed in all geographies.

We further segment the data according to whether changes in usage are positive or negative. Using our D–L specifications, with a slight modification, we test for effects on resistance due to positive and negative changes in usage, separately. Our hypothesis is that a positive change in usage will increase resistance whereas a negative change in usage will decrease resistance. Our empirical model to assess this hypothesis is given by

$$R_}}}}}}},t} = \mathop \nolimits_^4 }}}}}}},t - j} + \mu _i + \mu _c + \theta _t + \varepsilon _}}}}}}},t}}$$

(5)

where \(D_}}}}}}},t - j}\) is an indicator variable that indicates whether there was a positive or negative change above a specific threshold in the usage of antibiotics treating bacteria-class i, in country c, j periods away.