Patterns of antibiotic cross-resistance by bacterial sample source: a retrospective cohort study ================================================================================================ * Stacey S. Cherny * Michal Chowers * Uri Obolski ## Abstract **Background** Antimicrobial resistance is a major healthcare burden, exasperated when it extends to multiple drugs. While cross-resistance across multiple drugs is well-studied experimentally, it is not the case in clinical settings, and especially not while considering confounding. In addition, bacteria from different sample sources may have undergone different evolutionary trajectories, therefore examining cross-resistance across sources is desirable. **Methods** We employed additive Bayesian network (ABN) modelling to examine antibiotic cross-resistance in five major bacterial species, obtained from different sources (urine, wound, blood, and sputum) in a clinical setting, collected in a large hospital in Israel over a 4-year period. ABN modelling allowed for examination of the relationship between resistance to different drugs while controlling for major confounding variables. **Findings** Patterns of cross-resistance differ across sample sources. Importantly, even when positive cross-resistance exists between two antibiotics in different sample sources, its magnitude may vary substantially. **Interpretation** Our results highlight the importance of considering sample sources when assessing likelihood of antibiotic cross-resistance and determining antibiotic treatment regimens and policies. ## Introduction Antimicrobial resistance is a substantial healthcare burden, causing increased use of hospital resources, changes in treatment protocols, and excess morbidity and mortality.1 When bacteria are resistant to multiple drugs, the problem is exacerbated. As a result, studying bacterial resistance is a major area of inquiry, both by experimental models and clinical data. Many predictors of resistance in a clinical setting have been identified, including age, patients’ independence status, and previous antibiotic usage.2,3 Additionally, the interrelationship of resistance between different antibiotics is of primary clinical relevance in antibiotic selection. Cross-resistance of bacteria to different drugs is a well-established phenomenon.4–6 If a bacterial isolate is susceptible (or resistant) to a particular drug, it will often be susceptible (or resistance) to a different drug of the same drug family, although differences in the drugs may lead to imperfect associations. Discordant resistance, where susceptibility to one drug is associated with resistance to another (and vice versa), is repeatedly demonstrated in experimental studies.7–10 The paradigm typically involves exposure of bacterial cultures to increased doses of a single antibiotic over many generations of selection and then testing the adapted population for resistance to other antibiotics. Under such conditions, both cross-resistance and discordant resistance are found, with evolutionary mechanisms able to explain both. However, the situation doesn’t necessarily parallel bacteria living in human populations. As a result, it is perhaps not surprising that discordant resistance is rarely observed in clinical settings.4,11 Nonetheless, a recent paper proposed a method of analysis of clinical MIC data that is somewhat analogous to the experimental approach, that not only examines the existence of cross or discordant resistance, but attempts to establish direction of effect.11 They found some of the limited discordant (and cross) resistance present was unidirectional. Revealing which cross-resistance patterns occur in clinical settings is not trivial. Observational data are subject to confounding, which can substantially mask the underlying cross-resistance patterns. We have previously employed additive Bayesian network (ABN) modelling to deal with such confounding in clinical data.5 ABN models allow uncovering the association structure among a set of variables, while controlling for relevant confounders. However, the differences of cross-resistance patterns between sample sources have not been explored in depth previously, while adjusting for confounders. Bacteria from different sample sources may have undergone divergent evolutionary trajectories as a result of different selective forces, perhaps leading to distinct patterns of resistance. Such differences can have clinical implications on short- and long-term antibiotic treatment selection. In this study, we use a large clinical dataset obtained from an Israeli hospital over a 4-year period, to study cross-resistance in different sample sources. We employ ABN modelling to control for potential confounding variables and explore the resistance network structures by bacterial species and sources of the bacterial culture. This allows for unbiased examination of whether relationships between drug resistance differs by culture source. ## Materials and methods ### Data We obtained data pertaining to all positive bacterial cultures drawn in Meir Medical Center, Israel, a 740 bed hospital with approximately 60,000 admissions per year, from 2016-01-02 to 2019-12-31. The corresponding medical history, demographics, previous hospitalizations, and previous in-hospital antibiotic usage in the year prior to the infection, of patients from whom the cultures were drawn, were also available. Bacterial cultures were tested for antibiotic resistance for an array of antibiotics, and results of non-susceptibility and resistance were combined into a ‘resistant’ category. All antibiotics tested were systemic, with the exception of mupirocin, which is topical. Nonetheless, given we had data on bacterial resistance to mupirocin, it was still included in the analyses. Bacterial infections were considered nosocomial if cultures were drawn >48 hours after admission. A summary of these variables is presented in Tables 1 and 2. For our analyses, we selected the five bacterial species with the largest sample sizes available in the dataset: *Escherichia coli, Staphylococcus aureus, Pseudomonas aeruginosa, Klebsiella pneumoniae*, and *Proteus mirabilis*, in order of decreasing frequency. View this table: [Table 1:](http://medrxiv.org/content/early/2022/04/01/2022.03.31.22273223/T1) Table 1: Descriptive statistics of patients for *E coli* and *K pneumoniae* bacterial isolates, by source of sample View this table: [Table 2:](http://medrxiv.org/content/early/2022/04/01/2022.03.31.22273223/T2) Table 2: Descriptive statistics of patients for *P aeruginosa, P mirabilis*, and *S aureus* bacterial isolates, by source of sample ### Statistical analysis We selected which antibiotics to include in the analysis by keeping only those with minimal missing data and which did not reduce the number of complete cases appreciably (<10% loss). We performed some variable selection to assure stable statistical models with no perfect or near-perfect separation, by not including perfectly or near perfectly correlated antibiotics and selecting only antibiotics which contained a minimum of 3% resistance in each bacterial subsample. Antibiotics excluded from analysis due to high collinearity with included variables are presented in Supplementary Table S1, along with their range of tetrachoric correlations with the relevant included variables, across the various sample sources. Excluded antibiotics were usually, but not always, of the same drug family. For example, while ceftazidime, ceftriaxone, cefuroxime, and cefalexin are all cephalosporins and we often removed three of these four, ampicillin, a beta-lactam, was removed due to high correlation with the cephalosporins. This resulted in analysis of between three and five antibiotics for the five bacterial species cultured from the various sources, each analyzed separately. When constructing the ABN, the following covariates were included, in addition to the antibiotic resistance tests: demographic variables (*age, sex*, and *days hospitalized* in the previous year), presence of five medical conditions (*immunosuppression, dementia, diabetes*, chronic obstructive pulmonary disorder (*COPD*), chronic renal failure (*CRF*), and *obesity* (BMI over 30), binary culture type variables (*nosocomial* and *polymicrobial*), and a binary variable for antibiotic use in hospital in the previous year. Antibiotics used were first grouped into drug families and the six most frequent families across the entire dataset were included, along with whether any other antibiotics were taken (which did not belong to the six largest families). The six families were (in descending frequency of use) cephalosporins, beta-lactams (that are not cephalosporins), beta-lactamase inhibitors, nitroimidazoles, aminoglycosides, and fluoroquinolones, plus a seventh category including all other antibiotics. Table 1 presents summary statistics for *E coli* and *K pneumoniae*, for all variables used in our models, and Table 2 presents summary statistics for *P aeruginosa, P mirabilis, and S aureus*. For some models, certain covariates needed to be excluded from the analysis in order for the ABN models to converge due to insufficient sample size. In addition, for the *E coli* urine analysis, to allow the model to run within a reasonable amount of time (days rather than months), two variables were omitted which did not have a direct connection with resistance in the initial search. Excluded variables are noted in Tables 1 and 2. Data were analyzed using ABN modelling,12,13 with version 2.5.0 of the R package abn14 on an R 4.1.0 installation.15 Briefly, ABN modelling is a purely data-driven, exploratory approach, originating in machine learning, that is often used for hypothesis generation for causation among a set of variables. By searching likely networks linking a set of variables, a model for the dependency of the variables in the data can be inferred without making strong prior assumptions. This model shows which variables are directly connected, or linked, via arcs, and the coefficients of these arcs are directly analogous to the adjusted odds ratios obtained from multiple logistic regression analysis.12 While no assumptions of causal relationships are required, we restricted the model space to disallow causal paths that made no sense for our study, details of which are included in the Supplementary Methods. We followed a multi-step procedure: first identifying the most-likely network structure, using the abn package; then performing a parametric bootstrap, using JAGS to simulate 1000 datasets per model; then using ABN to reanalyze each simulated dataset, to correct for overfitting; and finally, re-estimating parameters using ABN, only allowing for connections seen in at least 50% of the analyses of the simulated datasets, and computing 95% Bayesian credible intervals (CIs) for each of the parameter estimates. The procedure we employed has been previously described and applied to antibiotic resistance data by us.5 One difference in the method of analysis in the present paper is that in order to speed up the parametric bootstrap, we dropped all variables that had no connection to any other variable in the initial network and reran the full procedure on the reduced dataset. This would have no effect on the results presented here. In addition, we limited the search for the maximum number of parents to eight, due to prohibitive computational time needed to conduct the search allowing more parents. Supplementary Table S2 contains the number of arcs present in the model both before and after bootstrapping, with the one model limited to 8 without reaching a plateau noted. To test whether arcs present in different sources of the same bacterial species were significantly different in magnitude, we randomly sampled from the posterior distributions of a given arc and computed the difference between a random pair posterior ln(OR) estimates, which yielded a distribution of ln(OR) differences. From that, we present the 95% CIs of the differences and examine whether zero was inside the intervals, indicating no significant difference in arc parameter estimates. The posterior estimates were computed at a granularity of 10,000 per parameter and 100,000 pairs were sampled with replacement to generate the empirical distribution of ln(OR) differences. ### Ethics The study was approved by the Institutional Review Board (Helsinki) Committee of Meir Medical Center. Since this was a retrospective study, using archived medical records, an exemption from informed consent was granted by the Helsinki Committee. ## Results In Figure 1, we present an example of an ABN estimated for *P aeruginosa* in a wound sample. The antibiotic resistance tests had direct links to *sex, polymicrobial infection*, and previous use of fluoroquinolones, beta-lactam inhibitors, and less-common antibiotics (*other)*. However, resistance to MEM was not linked to any variable. Thus, the network yields less biased associations between the different antibiotics by controlling for various patient covariates. ABNs for all bacterial species investigated are presented in Supplementary Figures S1-S14, with parameter estimates and 95% CIs presented in Supplementary Tables S3-S16. ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/04/01/2022.03.31.22273223/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2022/04/01/2022.03.31.22273223/F1) Figure 1: Additive Bayesian networks for *P aeruginosa* in wound samples. Variables in rectangles are binary and the one in a circle is continuous. Drug family variables refer to a particular drug taken in the prior year. Antibiotic resistance tests are shown in shaded boxes: CAZ, ceftazidime; GEN, gentamicin; MEM, meropenem, TZP, piperacillin/tazobactam. To better summarize our results, we present consensus graphs (Figure 2). These graphs present the relationships among drug resistances from all bacterial sources for a given bacterial species. Figure 2 contains 5 panels, A-E, representing the relationships among the antibiotic resistance tests in *E coli, K pneumoniae, P aeruginosa, P mirabilis, and S aureus*, respectively, obtained from fitting the ABN models (as the one presented in Figure 1; see Supplementary). The nodes in the consensus graphs are denoted with both the antibiotic resistance tests included in the ABN models and corresponding antibiotics excluded due to very high correlations with the included antibiotic, as seen from examination of the pairwise tetrachoric correlations (Supplementary Table 1). ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/04/01/2022.03.31.22273223/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2022/04/01/2022.03.31.22273223/F2) Figure 2: Consensus network graphs for the five bacterial species: *E coli* (A),*K pneumoniae* (B), *P aeruginosa* (C), *P mirabilis* (D), and *S aureus* (E). Line color depicts in which bacterial source model the link is present, with line width proportional to the ln(OR) (or mean ln(OR) for black lines). Orange=urine; red=wound; blue=aerobic blood; black=urine, blood, and wound, except in *P aeruginosa*, where black denotes urine, blood, and sputum. AMC, amoxicillin/clavulanate; AMP, ampicillin; CAZ, ceftazidime; CIP, ciprofloxacin; CRO, ceftriaxone; CXM, cefuroxime; ERY, erythromycin; GEN, gentamicin; IPM, imipenem; LEX, cefalexin; MEM, meropenem; MUP, mupirocin; OFX, ofloxacin; OXA, oxacillin; PIP, piperacillin; TZP, piperacillin/tazobactam; SXT, sulfamethoxazole/trimethoprim. In Figure 2, arcs are color-coded by sample source: orange arcs are in urine, red in wound, blue in aerobic blood, and black for arcs present in all three sources for a particular bacterial species. In all four of the species where ciprofloxacin and gentamicin were tested, they were linked, with the exception of *P mirabilis*, where they were linked in two of the three sample sources. Ceftazidime and gentamicin were also linked in four species: in *P aeruginosa*, in all three sources and in two of three sources in *E coli* and *P mirabilis*. Piperacillin/tazobactam was directly linked to ciprofloxacin in two of the three species for which it was tested and in two sample sources of *K pneumoniae*. Importantly, all direct connections were consistently positive throughout all antibiotics and all bacterial species. Finally, we sought to analyse whether the magnitude of the cross-resistance links differed between sample sources. Parameter estimates of the cross-resistance links, along with their 95% CIs, are presented in Figure 3, for those parameters estimated in more than one sample source of a particular bacterial species. The median differences between these parameter estimates in different sample sources, along with 95% CIs for the difference, within bacterial species, are presented in Supplementary Table S17. While many arcs were present in more than one culture source for a given species, only in three of 18 instances were the magnitudes of those parameter estimates not significantly different between sources compared. In *E coli*, the link between ofloxacin and ceftazidime was not significantly different in aerobic blood than in urine, as was the link between ceftazidime and gentamicin in urine versus wound. The third pair of estimates not significantly different was cefuroxime and ofloxacin in the *P mirabilis* cultures from urine versus wound. All other links present in two different culture sources were found significantly different between sources, as determined by the 95% CIs of the differences not including zero. ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2022/04/01/2022.03.31.22273223/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2022/04/01/2022.03.31.22273223/F3) Figure 3: Parameter estimates (ln(OR)) and their 95% credible intervals, for arcs between pairs of antibiotics, by bacteria and sample source, where both antibiotics appear in more than one sample source for a given bacteria. NS near a pair denotes that the two estimates were not significantly different, with * denoting within .01 from the null. The 95% CI of the estimate of the difference between all other pairs tested within a bacterial species did not contain zero. Line color depicts in which bacterial source the estimate was obtained: orange=urine; blue=aerobic blood; red=wound; green=sputum. AMC, amoxicillin/clavulanate; CAZ, ceftazidime; CIP, ciprofloxacin; CXM, cefuroxime; GEN, gentamicin; OFX, ofloxacin; TZP, piperacillin/tazobactam. *The 95% credible interval contains zero to one decimal place. ## Discussion In this work, we examined patterns of cross-resistance across bacterial culture sources, while adjusting for potential confounders, using ABNs. We found that the patterns and magnitude of cross-resistance among pairs of antibiotics may vary significantly between sample sources. Some antibiotics were found to have a nearly perfect association in our data and were retained irrespectively to any patient covariates. Such a grouping was formed by the cephalosporins ceftazidime, ceftriaxone, cefuroxime, and cefalexin in the three (*E coli, K pneumoniae, and P mirabilis*) bacterial species in which these antibiotics were used, for all tested sample sources. However, this is somewhat inconsistent with a recent study by Beckley and Wright, which found ceftazidime, ceftriaxone, and cefuroxime to be concurrently resistant, but the relationships were far from perfect.4 We note that Beckley and Wright used the mutual information score as their measure of association and hence the comparisons to our results are only qualitative. Furthermore, while we did find near perfect correlations among cephalosporins, there are instances in our data where resistance to an older generation of cephalosporin does not necessitate resistance to a newer generation, but almost never vice versa. Because those instances are relatively rare, the tetrachoric correlations were still very high. Other discrepancies observed between our study and Beckley and Wright4 could be due to several major differences between the analyses. Firstly, we controlled for an array of relevant covariates using ABNs. Secondly, we examined relationships among resistances stratified by sample sources. Finally, our data may be more homogenous, since they originated from a single center, in contrast with data used in Beckley and Wright, which originated from over 35 separate hospitals across the US.4 While all links identified in our study were positive or non-existent, 14 of 18 pairwise links were significantly different in their magnitude between different sample sources. This suggests that the extent of cross resistance found might depend on the selective pressures exerted by the environment the bacteria were sampled from. In several cases, these differences were dramatic. For instance, in *K pneumoniae*, the OR between gentamicin and ofloxacin ranged from around 6 in urine to nearly 40 in blood. In *E coli*, the same antibiotics varied in their ORs from 3 in urine to 11 in blood. Such differences can distinguish between a treatment being merely unlikely to succeed, given one type of resistance, to being almost futile. A limitation of our analyses is that sample size differed greatly both in bacterial species and sample sources. This translates into differing levels of statistical power across ABN analyses. For example, we observe the most arcs in the network for *E coli* obtained from urine, potentially because the sample size was largest in that analysis. For the three other bacterial species cultured from urine, it is apparent that there are more arcs in the models connecting antibiotic resistance tests as compared to other culture sources, again potentially due to urine having the largest sample sizes. Nonetheless, while sample sizes were modest in some groups, ample data was available for others, and the study examined all bacterial cultures obtained in a large hospital over a four year period, making it quite comprehensive. Furthermore, the patterns found in our analyses are not necessarily generalizable to other countries or even hospitals. Resistance patterns are determined by prescription patterns, patient demographics and a host of other factors, and are expected to differ across time and space.3,16 However, our results of different magnitudes of cross-resistance between sample sources have a biological rationale, and we therefore expect them to be more generalizable. To conclude, using a large clinical dataset obtained from admissions to a major hospital in Israel, over a 4-year period, we applied ABN modelling to control for potential confounding variables and explored resistance network structures stratified by bacterial species and source of the bacterial culture. This allowed for less biased examination of the relationships between drug resistance and culture sources. We found similar patterns of the existence of cross-resistance in different culture sources within bacterial species; however, the magnitude of those relationships differed substantially between sample sources. This highlights the importance of considering sample sources when assessing the likelihood of antibiotic cross-resistance. Future antibiotic prescription policies aiming to minimize collateral resistance should therefore be differentially determined by the source of infection. ## Supporting information Supplementary Material [[supplements/273223_file02.pdf]](pending:yes) ## Data Availability Data are proprietary but can be made available upon reasonable request from the authors. ## Contributors ### Funding This study was supported by the Israel Science Foundation (ISF 1286/21). ### Transparency declarations None to declare. ### Data sharing Data are proprietary but can be made available upon reasonable request from the authors. * Received March 31, 2022. * Revision received March 31, 2022. * Accepted April 1, 2022. * © 2022, Posted by Cold Spring Harbor Laboratory The copyright holder for this pre-print is the author. All rights reserved. The material may not be redistributed, re-used or adapted without the author's permission. ## References 1. Friedman ND, Temkin E, Carmeli Y. The negative impact of antibiotic resistance. Clin Microbiol Infect 2016; 22: 416–22. 2. Lewin-Epstein O, Baruch S, Hadany L, Stein GY, Obolski U. Predicting Antibiotic Resistance in Hospitalized Patients by Applying Machine Learning to Electronic Medical Records. Clin Infect Dis 2021; 72: e848–55. 3. Chatterjee A, Modarai M, Naylor NR, et al. Quantifying drivers of antibiotic resistance in humans: a systematic review. Lancet Infect Dis 2018; 18: e368–78. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-099(18)30296-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F04%2F01%2F2022.03.31.22273223.atom) 4. Beckley AM, Wright ES. Identification of antibiotic pairs that evade concurrent resistance via a retrospective analysis of antimicrobial susceptibility test results. Lancet Microbe 2021; 0. DOI:10.1016/S2666-5247(21)00118-X. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S2666-5247(21)00118-X&link_type=DOI) 5. Cherny SS, Nevo D, Baraz A, et al. Revealing antibiotic cross-resistance patterns in hospitalized patients through Bayesian network modelling. J Antimicrob Chemother 2021; 76: 239–48. 6. Obolski U, Dellus-Gur E, Stein GY, Hadany L. Antibiotic cross-resistance in the lab and resistance co-occurrence in the clinic: Discrepancies and implications in E.coli. Infect Genet Evol J Mol Epidemiol Evol Genet Infect Dis 2016; 40: 155–61. 7. Imamovic L, Sommer MOA. Use of Collateral Sensitivity Networks to Design Drug Cycling Protocols That Avoid Resistance Development. Sci Transl Med 2013; 5: 204ra132–204ra132. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTE6InNjaXRyYW5zbWVkIjtzOjU6InJlc2lkIjtzOjE0OiI1LzIwNC8yMDRyYTEzMiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIyLzA0LzAxLzIwMjIuMDMuMzEuMjIyNzMyMjMuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 8. Pál C, Papp B, Lázár V. Collateral sensitivity of antibiotic-resistant microbes. Trends Microbiol 2015; 23: 401–7. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.tim.2015.02.009&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25818802&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F04%2F01%2F2022.03.31.22273223.atom) 9. Lázár V, Nagy I, Spohn R, et al. Genome-wide analysis captures the determinants of the antibiotic cross-resistance interaction network. Nat Commun 2014; 5: 4352. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ncomms5352&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25000950&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F04%2F01%2F2022.03.31.22273223.atom) 10. Munck C, Gumpert HK, Wallin AIN, Wang HH, Sommer MOA. Prediction of resistance development against drug combinations by collateral responses to component drugs. Sci Transl Med 2014; 6: 262ra156–262ra156. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTE6InNjaXRyYW5zbWVkIjtzOjU6InJlc2lkIjtzOjE0OiI2LzI2Mi8yNjJyYTE1NiI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIyLzA0LzAxLzIwMjIuMDMuMzEuMjIyNzMyMjMuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 11. Zwep LB, Haakman Y, Duisters KLW, Meulman JJ, Liakopoulos A, van Hasselt JGC. Identification of antibiotic collateral sensitivity and resistance interactions in population surveillance data. JAC-Antimicrob Resist 2021; 3: dlab175. 12. Kratzer G, Lewis FI, Comin A, Pittavino M, Furrer R. Additive Bayesian Network Modelling with the R Package abn. ArXiv Prepr ArXiv191109006 2019. 13. Lewis FI, Ward MP. Improving epidemiologic data analyses through multivariate regression modelling. Emerg Themes Epidemiol 2013; 10: 4. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/1742-7622-10-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23683753&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2022%2F04%2F01%2F2022.03.31.22273223.atom) 14. Kratzer G, Pittavino M, Lewis FI, Furrer R. abn: an R package for modelling multivariate data using additive Bayesian networks. 2019 [https://CRAN.R-project.org/package=abn](https://CRAN.R-project.org/package=abn). 15. R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing, 2020 [https://www.R-project.org/](https://www.R-project.org/). 16. Baraz A, Chowers M, Nevo D, Obolski U. Stable temporal relationships as a first step towards causal inference: an application to antibiotic resistance. 2022; : 2022.01.31.22270156.