Blog

The School of Mathematics Understanding the Relationship between Hospital Prescription Rates & Antimicrobial Resistance by Ravi Patel, S1889112 July 2019 Supervised by Dr. Gail Robertson, Dr. Meghan Perry, Dr. Serveh Sharifi Far & Dr. Vanda Ina′cio de Carvalho Own Work Declaration All work contained within is my own, except where otherwise referenced. Contents 1 Introduction 3 1.1 Background & Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Software and Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.5 Terminology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 EDA 4 2.1 Single Variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Multiple Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 EDA Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3 Model 9 3.1 Ordered Logit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Resistance Score as the Latent Variable . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Bayesian Specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 4 Bayesian Workflow 12 4.1 Prior Predictive . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4.2 Model Fit and Algorithm Diagnostics . . . . . . . . . . . . . . . . . . . . . . . . 13 4.3 Posterior Predictive (Hypothesis Model) . . . . . . . . . . . . . . . . . . . . . . . 14 4.4 Model Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.5 Posterior Predictive (mc41) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 5 Results 23 5.1 DDD and AMR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.1.1 Route . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2 Drug-Specific Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5.3 Class Group Heterogeneity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 5.4 Comments on Other Fitted Models . . . . . . . . . . . . . . . . . . . . . . . . . . 26 6 Model Discussion 27 6.1 Concerns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 6.2 Alternative Model: Zero-one inflated Beta (ZOIB) . . . . . . . . . . . . . . . . . 27 7 Conclusion 28 Appendices 32 A Data Preparation 32 B Models Compared 33 C R Session Information 34 D Code & Data Access 35 D.1 Stan Code Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 1 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. 6 Figure 6: DDD-Proportion Relationship Figure 7: Proportion over time by Class Group Figure 8: TRAK-Age Relationship 7 2.3 EDA Summary We have 1. Transformed the proportion to a set of ordered groups. 2. Grouped drug classes based on their level of resistance for initial modelling. 3. Standardised the average age, and logarithmised the DDD/1000 OBD. 4. Decided to include log DDD as a population level effect, with an interaction on the route. 5. Included age as a population-level effect 6. Not included week due to a lack of any visible relationship in the data. 7. Not included TRAK due to its strong relationship with age. 8 3 Model We use a Bayesian ordered logit with varying intercepts. This section provides an overview of ordered logit, discusses the latent variable, and specifies our priors. 3.1 Ordered Logit This subsection gives an overview of ordered logit. Options for a more comprehensive treatment include Walker & Duncan (1967), McCullagh (1980), and Greene (2012) (ordered probit). In ordered logit, our response variable is the binned proportion (y). Ordinal regression supposes there is a latent variable (y ∈ R) which leads to our bins by the process below. We have thresholds cj , j ∈ {0, 1, . . .K = 7}, where c0 = ∞, and c7 =∞. y = 0 if y ∈ ( ∞, c1] (0, 0.25] if y ∈ (c1, c2] (0.25, 0.5) if y ∈ (c2, c3] 0.5 if y ∈ (c3, c4] (0.5, 0.75] if y ∈ (c4, c5] (0.75, 0.99) if y ∈ (c5, c6] 1 if y ∈ (c6,∞) For the cumulative logit, we assume the latent variable follows a logistic distribution with scale 1 (variance pi 2 3 ) and mean η (the usual linear predictor x ′β, excluding the intercept). We define the cumulative distribution function (CDF) of the standard logistic distribution as Λ. y ～ Logistic(η, 1) Λ(x) = exp(x) 1 + exp(x) Figure 9: Latent Variable Categorisation The probability y is in category k given η is Pr(y = k) = Pr(y ∈ (ck 1, ck)) = Λ(ck η) Λ(ck 1 η). This concept is illustrated in Figure 9. For example, the yellow region is the probability the latent variable is between c2 and c3, i.e. the probability of the observed category being category 3. We therefore define the pdf and likelihood as: 9 f(y) = K∏ j=1 [Λ(cj η) Λ(cj 1 η)]I(y=j) (3.1) L = n∏ i=1 K∏ j=1 [Λ(cj ηi) Λ(cj 1 ηi)]I(yi=j) (3.2) Where ηi refers to covariates changing by observation. Define γj = Pr(y ≤ cj). We then have γj = Pr(y ≤ cj) = Λ(cj η). Applying Λ 1: Λ 1(γj) = logit(γj) = log ( γj 1 γj ) = cj η (3.3) We can interpret the coefficients in η in relation to the log odds. Suppose β1 > 0 is the coefficient on x1. A one unit increase in x1 increases the log-odds of a higher category by β1. This interpretation is independent of the category, which is where “proportional odds” comes from. 3.2 Resistance Score as the Latent Variable It may be concerning that we know the bins were created using the proportion, which is non-logistic. The application of ordered logit therefore requires a level of abstraction. We posit an underlying “resistance score” that is impossible to measure, but can be partially captured by the proportion. Our latent variable is mapped to the proportion, which we transform to our bins. We suppose there is a function f which maps from the reals to the closed interval [0,1]. Additionally, there is a function g which maps from the interval (c1, cK 1) to the open interval (0,1). We define f below, along with a mathematical statement of the function mappings. f :R 7→ [0, 1] g :(c1, cK 1) 7→ (0, 1) f(y ) = 0 if y ≤ c1 g(y ) if y ∈ (c1, cK 1) 1 if y ≥ cK 1 g′(y ) > 0, to ensure higher scores lead to higher proportions. The latent variable is 0 if below c1, and 1 if above cK 1. Otherwise, it passes through some function g to an intermediate proportion. We divide this intermediate proportion into 5 regions to get our binned observations. There also exists a discrimination parameter equal to the inverse standard deviation, which allows heterogeneity in the latent variable. It is modelled on the log-scale, since it must be positive. We allow for heterogeneity by each class group, using the “A” group as the reference category. For more information, see Bu¨rkner & Vuorre (2018). 3.3 Bayesian Specification An alternative parameterisation of the latent variable equation accounting for group-level effects is: y = x′β + z′u+ ～ Logistic(0, 1) (3.4) Hence, we can redefine (3.3) as: logit(γj) = cj x′β z′u (3.5) Where z defines the groups, and u defines the group-level parameters. Our priors and hyperpriors are as follows: 1. Population-level parameters βj are N(0, 3 2) priors to be weakly informative. Recall that the standard deviation of the standard logistic is pi√ 3 ≈ 1.8, we increase this value to 3 to add some extra uncertainty to reduce the chance of underfitting. 10 2. The thresholds cj have N(0, 3 2) priors for the same reason. 3. The standard deviation of group-level parameters were modelled as Exponential(1), following the usage by McElreath (2019). 4. brms uses the non-centered parameterisation, so mean hyperpriors were not explicitly set. 5. For the discrimination, we allow heterogeneity over each class grouping, and place a N(0,3) prior on each class group intercept. 6. On the occasions we model varying slopes and intercepts, we use the LKJ(1) (Lewandowski et al., 2009) prior on the correlations. 11 4 Bayesian Workflow The Bayesian workflow will be presented based on our model hypothesised following section 2. In general, we follow the steps of analysis from Gabry et al. (2019). This model is presented in brms syntax for conciseness. Proportion Gr ～ Age Standardised + log DDDp1kOBD*Route + (1|Class Gr) disc ～ 0 + Class Gr The first equation can be thought of in latent variable form as: y = β1Age+ β2lDDD + β3NonInvasive+ β4lDDD NonInvasive+ α[Class Group] + (4.1) We have population level effects on the standardised age and the log DDD, with a group-level intercept on the grouped class. We allow heterogeneity by each class grouping. The group intercept on the grouped class is subject to change to the ungrouped class and the drug name. Howevver disc (discrimination) will always be dependent on the grouped class due to sparsity issues with the ungrouped class and drug name. 4.1 Prior Predictive Using proper priors allows us to generate samples from the prior predictive. Our goal is to ensure that the observed data is feasible under our prior. Figure 10 shows our prior is not urealistic given our key test statistics (number of 0s, 1s, 0.5s), with the sampled data (yrep) capturing the test statistic within the bounds of the histogram. Similarly, the empirical CDF (ECDF) shows the data is unlikely but not impossible under the prior. Figure 10: Prior Predictive Plots 12 4.2 Model Fit and Algorithm Diagnostics We ran 3 chains, with warmup of 2000, and 5000 total iterations per chain, meaning 9000 post warmup samples with combined chains. This may seem small for those familiar with Gibbs sampling, however the NUTS (Hoffman & Gelman, 2014) allows for convergence with fewer samples, as the effective size per iteration is much higher (Bu¨rkner, 2017). The model had 0 divergent transitions1. The split R should be below 1.1 for all parameters (Gelman et al., 2013), and our model is below this threshold as per Figure 11. The minimum effective size ratio is above 0.15, meaning our effective sizes are all above 1000 which is regarded as usually being enough for stable estimates (Bu¨rkner & Vuorre, 2018). Trace plots display no mixing issues. Figure 11: Convergence Diagnostics and ESS 1Divergent transitions are an indicator of potential bias in the HMC algorithm, and are a more sensitive diagnostic than the usual R measure. For more information, see the online case study by Betancourt (2017). 13 4.3 Posterior Predictive (Hypothesis Model) This section explores whether the fitted model would be able to generate a dataset similar to what we observed in reality. We do this primarily by computing test statistics, and checking these against the true data. We examine the test statistics split by the class group, to get a more precise look at the data. Posterior predictive checks can be thought of as the Bayesian analogue to goodness-of-fit tests. Figure 12 shows the count of each proportion group segmented by the class group. We see that the counts within each class largely match up to the observed values. The 2nd class in group F however looks suspect, with a relatively low level of uncertainty, yet being above the true count indicated by the bar. There is a substantial amount of variation in groups E to G, suggesting there is a lot of uncertainty around these posterior frequencies. Figure 12: Frequency Plot by Grouped Class Figure 13 shows the density of the overall dataset appears to be well fitted by the simulations from the posterior predictive. Figure 13: Posterior ECDF 14 Figure 14 shows the distribution of the counts of 0s, 0.5s and 1s within each class group. In general, the model captures the test statistics well. However the number of 1s in B is slightly underestimated. The number of 0.5s is not unreasonable for any groups, however B, E, and F are slightly off center. On the whole, these test statistics show an adequate fit to the data. Figure 14: Test Statistics 15 Figure 15 shows a similar frequency plot to Figure 12, but now the groups are the true observations. In terms of direct predictions the model does poorly, proposing similar patterns within each bin. That is each bin has a peak at 0 and 1, with a wide normal shape on (0,1). Our goal may be inference, however this lack of predictive ability is concerning as it suggests that the data is too variable to predict, meaning we should be sceptical of our coefficients. Examining the predicted linear component (i.e. the values on the latent scale) in Figure 16 we see that, as desired, higher bins have a higher value of the latent variable. This breaks down for the 7th group (proportions of 1), but the general trend is promising. This suggests there is some ability to distinguish between proportion groups, but that there is too much variance for specific predictions to be made. Figure 15: Frequencies vs True Values Figure 16: Linear Predictor vs. True Group 16 4.4 Model Comparison We compare our hypothesised model to alternative formulations to establish room for improvement, using Pareto-smoothed importance sampling leave-one-out cross-validation (PSIS-LOO) (Vehtari et al., 2017a). PSIS-LOO has been established as preferable to the widely applicable information criterion (WAIC) developed by Watanabe (2013) – another common model comparison tool. Both aim to approximate the expected log pointwise predictive density (elpd), however WAIC has a larger bias, and the bias-variance tradeoff is harder to control (Vehtari et al., 2017b). In general, we interpret a higher elpd as better, and a model as clearly better if the difference in elpd is greater than 2 standard errors. However this is only a general heuristic, since standard error estimates have high variance (Piironen & Vehtari, 2017). Details of the models compared are appended. The idea is to use a different group for the group-level effect (drug class, drug name), allow varying slopes, and remove the interaction between DDD and the route. We also use an overfit model to get a rough estimate of the maximum predictive potential on the data under the ordered logit, by plotting the equivalent of Figure 15 in Figure 17. Our hypothesised model is mp1, while comparison models begin with mc. The overfit model is mc1, varying slopes is mc2, mc3x models use the ungrouped class as the varying intercept, and mc4x models use the drug name. No divergent transitions were reported, R was below 1.1 for all models, effective sample sizes were all above 1000, and k values (a diagnostic for PSIS-LOO) were all less than 0.7 except the overfit model, which had 6 “bad” values (in (0.7, 1]) and 1 “very bad” value (above 1). Model elpd diff se diff elpd loo se elpd loo p loo se p loo looic se looic mc41 0.00 0.00 -3438.12 51.90 38.24 1.75 6876.24 103.80 mc43 -0.64 1.71 -3438.76 51.92 37.79 1.72 6877.51 103.84 mc44 -0.89 1.64 -3439.01 51.93 37.33 1.74 6878.02 103.86 mc42 -1.23 1.89 -3439.34 51.92 39.16 1.76 6878.69 103.85 mc1 -9.90 5.08 -3448.02 52.11 59.48 3.81 6896.05 104.22 mp1 -23.73 9.80 -3461.85 50.91 20.36 0.68 6923.69 101.82 mc2 -29.85 10.30 -3467.97 50.94 22.37 0.70 6935.93 101.87 mc31 -30.03 9.45 -3468.15 51.02 27.24 0.91 6936.31 102.04 mc34 -36.42 10.07 -3474.54 51.03 26.27 0.88 6949.08 102.06 mc33 -36.70 10.15 -3474.82 51.01 26.25 0.88 6949.64 102.02 mc32 -37.17 9.98 -3475.29 51.03 28.94 0.92 6950.58 102.06 Table 2: Model Comparison Table 2 implies a varying intercept for the class group is poor relative to a name-varying intercept, but better than the raw class. In general, population-level effects for the age with an interaction between the route and DDD performed better than alternatives. This can be seen as mc41 and mc31 are better than others within their group, and mp1 is better than mc2. The overfit model performs better than our hypothesised model, however this elpd approximation is likely very biased due to the high k values. The prediction barplot for the overfit model can be seen in Figure 17. Again, the raw predictions themselves do poorly, suggesting there is too much variance in the dataset to generate good predictions. Similarly, Figure 18 has the same pattern as previously. Given the similarity of the “best” model (mc41) to our hypothesis (the only difference is the group-level intercept), we briefly review the posterior predictive for mc41, and proceed to results. 17 Figure 17: Frequencies vs True Values (Overfit) Figure 18: Linear Predictor vs True (Overfit) 18 4.5 Posterior Predictive (mc41) First, consider Figure 19. Though at first glance there are many “wrong” estimates (e.g. the top left panel), this is just a result of shrinkage to the mean that comes with using group-level effects. Since there are not many obserations in certain groups, their intercept is pushed towards the mean as a result of the common distribution shared by each group. This is how varying intercepts/slopes prevent overfitting in hierarchical models. McElreath (2019) refers to this as adaptive regularisation. See also Gelman et al. (2013), or Gelman & Hill (2007). The ECDF in Figure 20 looks to fit the data well. The test statistics for 0s, 1s, and 0.5s in Figures 21, 22, and 22 largely seem appropriate. Some exceptions are the number of 1s for vancomycin, and possibly flucloxacillin, and the numbers of 0.5s for vancomycin, cefalexin, ceftriaxone, and ciprofloxacin. The barplot of frequencies against true values (Figure 24), and the latent variable plot (Figure 25) are the same story as previously. Poor actual predictions, but a mostly desirable pattern in the latent variable. Figure 19: Frequency by Name 19 Figure 20: ECDF Figure 21: Test Statistic by Name (Number of 0s) 20 Figure 22: Test Statistic by Name (Number of 1s) Figure 23: Test Statistic by Name (Number of 0.5s) 21 Figure 24: Predictions vs Truth Figure 25: Linear Predictor Value by True Value 22 5 Results 5.1 DDD and AMR The results for the population-level effects are in Table 3. The distributions for the key coefficients (i.e. involving DDD and route) are shown in Figure 26. Under our model, there is a 95% probability of a positive association between resistance and the DDD for invasively administered drugs. Specifically, we estimate a unit increase in the log DDD for invasive treatments increases the log-odds of a higher proportion group by 0.24. It may be more intuitive to look at the latent scale, where Table 4 will assist us. We are still looking at invasively administered drugs. Looking at the 2nd column, we can say a unit increase in the log DDD closes the distance between the first and second threshold by 26%. Inverting this number, we see that to transition the whole gap, the log DDD needs to increase by 3.81. From this, we can calculate transition ratios – the multiplier needed on raw DDD to push the linear predictor from the bottom of one group to the bottom of the next group up. This is the exponent of the values in the Difference/βlDDD column. For example, to transition from the bottom of group 2 to the bottom of group 3 given a specific value of DDD x1, the difference x2 x1 needs to be such that x2x1 is 45.2 This suggests that although the effect of DDD on resistance is likely positive, it is not overly important materially – particularly considering the middle transition ratios. Recalling Figure 1 – we can see that there is a low proportion of the higher values of raw DDD. The proportion of DDD values larger than 400 is only 6%, suggesting that DDD changes are unlikely to be able to discriminate between groups, since such large values are so rare. However, it appears that transitioning from group 6 to 7 is more strongly determined by DDD (relatively), with a multiplicative effect of 24 for a bottom-to-bottom transition. In general, the importance of DDD in group transitions appears to be higher at more extreme proportions, in terms of the ability to change the proportion group. We might hypothesise that there is some medical tipping-point phenomenon near the more extreme proportions, however this would require more domain expertise. Parameter Mean SD 2.5% 5% 95% 97.5% Pr(>0) Pr(<0) c1 1.94 0.81 0.42 0.66 3.30 3.58 0.99 0.01 c2 2.85 0.85 1.32 1.55 4.33 4.64 1.00 0.00 c3 4.11 0.96 2.45 2.66 5.79 6.18 1.00 0.00 c4 5.82 1.17 3.76 4.05 7.86 8.35 1.00 0.00 c5 6.96 1.34 4.60 4.93 9.32 9.88 1.00 0.00 c6 7.72 1.46 5.14 5.50 10.28 10.89 1.00 0.00 βAge 0.14 0.11 -0.07 -0.03 0.32 0.37 0.90 0.10 βlDDD 0.24 0.15 -0.04 0.00 0.51 0.57 0.95 0.05 βNonInvasive 1.75 0.97 -0.01 0.25 3.45 3.84 0.97 0.03 βlDDD NonInvasive -0.36 0.21 -0.81 -0.72 -0.04 0.02 0.03 0.97 βlDDD + βlDDD NonInvasive -0.12 0.15 -0.42 -0.37 0.12 0.17 0.20 0.80 Table 3: Results Table (Population Effects) Difference βlDDD/Difference Difference/βlDDD Transition Ratios c2 c1 0.91 0.26 3.81 45.04 c3 c2 1.26 0.19 5.27 193.99 c4 c3 1.71 0.14 7.15 1269.58 c5 c4 1.15 0.21 4.79 120.90 c6 c5 0.76 0.32 3.16 23.65 Table 4: Group Transition Table 2log(x2) log(x1) = 3.81 x2x1 = exp(3.81) ≈ 45 23 For non-invasive administration, there is a 97% probability the association between DDD and proportion is smaller than for invasive administration (βlDDD NonInvasive). The overall ceteris paribus3 effect of an increase in log DDD for a non-invasively administered drug is, on average, a reduction of 0.12 on the latent scale, with an 80% probability the association is negative, hence a 20% probability it is positive. Therefore, there is lots of uncertainty around the association of DDD for non-invasive drugs, so we should not commit to any conclusions. Summarising, we can note that the association between DDD and the proportion: 1. Is positive for invasively administered drugs (95%, βlDDD) 2. Is smaller for non-invasively administered drugs (97%, βlDDD NonInvasive) 3. May be negative for non-invasive drugs (80%, βlDDD + βlDDD NonInvasive) Figure 26: Posterior Area Plots (DDD and Route) 5.1.1 Route Ceteris paribus, non-invasive administration has a much higher chance of falling into a higher bin, with a mean of 1.75, and a 97% probability of non-invasive administration being associated with higher proportions. Reviewing the thresholds, it appears that the administration route being non-invasive is enough, on average, to achieve a bottom to bottom transition for each group – suggesting it is an important predictor of AMR. 3All else being equal 24 5.2 Drug-Specific Effects Posterior interval plots of the varying intercepts are shown in Figure 27, with lines at -3 and 3. It appears a large part of the hospital resistance is being driven by a few drugs, with only 5 drugs having 90% symmetric credible intervals not containing 0. The highest magnitudes come from amoxicillin, doxycycline, and trimethoprim (minocycline is excluded due to its wide intervals). On the other side, linezolid, meropenem, and metronidazole all have a strong negative relationship with resistance. Figure 27: Posterior Interval Plots (Drug Name) 25 5.3 Clas