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.