Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
between Hospital Prescription Rates
& Antimicrobial Resistance
S1889112
Executive summary
Antimicrobial resistance (AMR) is a big threat, with a low estimate placing the potential
death toll due to AMR in 2050 at 10,000,000 (O’Neill, 2016). Studies at the individual
(Costelloe et al., 2010) and population (Goossens et al., 2005) level show a positive association
between drug use and AMR – our goal is to investigate this association in a hospital setting.
Using a Bayesian ordered logit model with drug-varying intercepts, we conclude there is a 95%
probability of a positive association between drug use and resistance – specifically for
invasively administered drugs. Additionally, there is a 97% probability that this association is
smaller when considering non-invasive administration. Both these effects are very small in
magnitude, and poor predictive performance means these results should be looked upon
sceptically. We also propose a potential idea to expand the zero-one-inflated beta model that
could account for ordering, and may be a good model for similar such studies in future.
2
1 Introduction
1.1 Background & Motivation
Antimicrobial resistance (AMR) is similar to antibiotic resistance, also including non-bacterial
microbes such as viruses, parasite, and fungi (World Health Organisation (WHO), 2018).
Their use is widespread in common medical procedures such as C-sections, joint replacements,
chemotherapy, and organ transplants (Roope et al. (2019), WHO (2018)), making the rise of
resistance a serious threat.
Cassini et al. (2019) claim that in 2015, resistant infections accounted for between 28,480 and
38,430 deaths, and between 768,837 and 989,068 disability-adjusted life years lost (i.e. a year
lost in healthy life) at 95% confidence. O’Neill (2014, 2016) estimates there will be 10,000,000
deaths from AMR in 2015, and suggests this may even be an underestimate.
Furthermore, we would expect an economic consequence due to both morbidity and mortality.
The World Bank (2017) estimate in a pessimistic (but not unrealistic) scenario, that the
economic impact of AMR could be similar to that of the 2008 financial crisis.
1.2 Objective
This paper aims to provide a small contribution to the literature on AMR by examining the
effect of prescription rates in a hospital setting. Previous research has focused on individual
level resistance, such as Costelloe et al. (2010) who conducted a meta-analysis finding bacterial
resistance developed to administered antibiotics in primary care patients. Similarly, at a
population level, Goossens et al. (2005) showed an association between high-consuming
countries (measured in DDD per 1000 inhabitants per day), and antibiotic resistance.
1.3 Data
Data was collected over 3 months from the Western General Hospital (WGH) in Edinburgh.
Each row represents an antibiotic prescription case within a ward, within a week (a
ward-week). Variables of interest are the drug name (e.g. amoxicillin), the class the drug falls
into (e.g. Broad penetration β-lactams), the route of drug administration, the TRAK code (a
ward identifier), the average age of patients in a ward-week (henceforth just called age), the
defined daily dose of drugs administered in a ward week per 1000 occupied bed days
(DDD/1000 OBD), and the proportion of clinical isolates tested in a ward week that showed
resistance to the prescribed antibiotic (resistance proportion). DDD/1000 OBD can be viewed
as a measure of the drug strength in a ward-week for an antibiotic that accounts for the
number of patients in the ward. Our interest is the association between DDD/1000 OBD and
the resistance proportion. After data preparation (see appendix), we were left with 2548 rows.
1.4 Software and Algorithms
Models were fit with brms (Bu¨rkner, 2018), an R (R Core Team, 2019)) interface to the Stan
language (Stan Development Team, 2018) which implements Hamiltonian Monte Carlo (HMC)
(Duane et al., 1987; Neal, 2012) and an improvement – the No U-Turn Sampler (NUTS)
(Hoffman & Gelman, 2014). Plots were created with ggplot2 (Wickham, 2016), cowplot
(Wilke, 2019), and bayesplot (Gabry & Mahr, 2019). Data manipulation was performed
using dplyr (Wickham et al., 2019) and data.table (Dowle & Srinivasan, 2019), while tables
were created with xtable (Dahl et al., 2019). Session information for R is in the appendix.
1.5 Terminology
I avoid the terms ‘fixed effects’ and ‘random effects’ following the advice of Gelman & Hill
(2007). I instead follow brms terminology and refer to population-level effects, and group-level
effects. These are also sometimes referred to as ‘constant’ and ‘varying’ (respectively) for ease
of writing, as per Gelman (2005).
3
2 EDA
2.1 Single Variable
Figure 1 shows the age is approximately normal (though with a second mode in the 80s), while
the DDD is right-skewed due to extreme values. Both have large standard deviations, which
can cause issues in sampling and complicate prior specification. Consequently, we standardise
the age, and take the logarithm of the DDD/1000 OBD.
Figure 1: Continuous Covariate Histograms
Figure 2 illustrates the relationship of proportion with factor variables, and the week. Red
points signify the standard deviations, and annotated numbers refer to sample sizes. We see
that non-invasive administration is associated with higher resistance, but no visible time trend.
Some wards appear to have less resistance than others, however there is no natural grouping in
the wards as all of the error bars (1 standard error) contain adjacent wards.
Drug classes show separation into different groups based on the proportion. We group based
on both observations, and proportion differences. For instance we group the top 2 classes as
their means and errorbars overlap. However we do not include “BL - Broad Pen” in the group
as it has enough observations to be its own group. The groupings are shown in Table 1. This
enables cleaner visualisation, and will be used as a starting point for modelling.
An important feature is the volatility in the data, with standard deviations generally above
30% for subgroups. This makes it difficult to generate good predictions on the data, as there is
too much noise relative to the small dataset.
Group Drug Classes
A Anti-TB; Chloramphenicol; Nitroimidazole; Oxalidinone
B β-Lactam - Carbapenem; Polymixin; Glycopeptides/Lipopeptides
C Other; Aminoglycoside
D β-Lactams (Extended Pen; Narrow Pen; Cephalosporin); Quinolones
E Macrolide/Lincosamide
F β-Lactam - Broad Pen
G Tetracycline; Trimethoprim and sulphur based drugs
Table 1: Drug Class Groupings
4
Figure 2: Week and Discrete Covariates vs Proportion
Figure 3 visualises the proportion in 3 ways. The left panel exhibits a large number of 0s, 1s,
and possibly 0.5s. This is clearer in the middle panel, which shows a frequency barplot for each
individual proportion. There are spikes at very specific values, such as 0, 1, 0.5, 0.25 and 0.33.
As such, we can capture most of the shape in the proportion by grouping them into relevant
intervals. This yields the barplot in the right panel, motivating the use of ordinal regression.
Figure 3: Proportion Distribution
2.2 Multiple Variables
Recall that our interest is the association of DDD with resistance, so we include it in our
model. We will also include some sort of drug-specific effect – the grouped class, raw class, or
drug name – based on the class-level proportion differences seen earlier.
Figure 4 shows the estimated linear relationship between the proportion and the DDD,
segmented by route. Invasively administered drugs seemingly have a lower resistance
proportion for low levels of DDD, but increases in DDD have a stronger bearing on resistance
than for non-invasively administered drugs. As such, we will include the route in our model,
and interact it with age to account for the different slopes.
Figure 5 suggests a negative relationship between age and DDD, which may be expected due
5
Figure 4: Proportion-DDD Relationship by Route
to the increased sensitivity of the body to certain drugs in elderly patients as suggested in
Hughes (2002) and Turnheim (1998). This relationship slightly varies by class, but not
enought to motivate varying slopes. Adjusted R2 values range from 0 to 0.12 for each class,
with all coefficients significant at the 5% level except for group E. All slope coefficients were
negative, ranging from -0.05 to -2.99.
We include age as a potential confounder, as it may have a direct influence on the proportion
due to human-human resistance transmission (Holmes et al., 2016) in the elderly as a result of
nursing homes. Age is therefore considered a potential “fork” confounder as discussed in
McElreath (2019) and Pearl (2009).
Figure 5: Age-DDD Relationship by Class Grouping
Figure 6 shows a slight positive slope for groups B, C, and D but flat slopes for E, F, and G.
However due to the level of variation in the confidence intervals and the small differences, we
posit population-level effects, as any group-level effects will be too subtle with our noisy data.
Figure 7 suggests that the week has little influence on the proportion. There is a downward
slope for group E, however there is a large degree of variation. Therefore we consider week to
have no effect on the proportion.
Figure 8 points to a strong relationship between the age and the ward, with the ward
seemingly a good predictor of the age due to the low age variation within most wards. A linear
regression of age on TRAK yields an adjusted R2 of 0.9329, suggesting including both TRAK
and age would be redundant. We select age for theoretical reasons, and to avoid including too
many parameters.