## Abstract

Methicillin-resistant *Staphylococcus aureus* (MRSA) colonizes multiple body sites, and carriage is an important risk factor for MRSA infection. Successful decolonization reduces disease incidence; however, decolonization protocols vary in the number of body sites targeted, and the impact of site-specific treatments is not well understood. Here, we used data from a randomized controlled trial (RCT) of MRSA decolonization using chlorhexidine body and mouth wash and nasal mupirocin to quantify the contribution of each treatment component to the success of the protocol. We estimated mouthwash as the least effective treatment component and the combined effect of MRSA clearance at the nares, skin, and wound as 93% (90% credible interval 85%-99%) of the full decolonization. Our model can estimate the effectiveness of hypothetical treatments *in silico* and shows enhancing MRSA clearance at nares will achieve the largest gains. This study demonstrates the use of machine learning to go beyond what is typically achieved by RCTs, facilitating evidence-based decision-making to streamline clinical protocols.

## 1 Introduction

Methicillin-resistant *Staphylococcus aureus* is a common antimicrobial-resistant pathogen in community and healthcare settings[1, 2], causing an estimated 320,000 infections in hospitalized patients and over 10,000 deaths in the United States in 2017[3]. Progress in reducing invasive MRSA infections has slowed, underscoring the importance of continued innovation and effort to prevent disease[4]. As MRSA carriage is a major risk factor for invasive disease, efforts at prevention center on the promotion of decolonization and body hygiene as well as environmental cleaning[5]. The most common *S. aureus* carriage site is the anterior nares, and it also colonizes the perineum and groin, the axilla, the pharynx, as well as other body sites[6, 7]. While the anterior nares have been identified as a key reservoir for transmission and a risk factor for invasive disease[8], the extent of interaction among colonization sites and the importance of addressing each site in decolonization protocols remain unclear. Decolonization protocols primarily target the nares and can include other sites as well, but it would be ideal to simplify, where possible, without significantly impacting effectiveness. Achieving this goal requires a detailed understanding of the dynamic relationship among colonization sites and the extent to which colonization at one site affects colonization at other sites.

The CLEAR (Changing Lives by Eradicating Antibiotic Resistance) Trial demonstrated that the use of a post-discharge decolonization protocol in MRSA carriers reduces infection and hospitalization rates[9]. In the trial, 2,121 study participants were randomized into two groups to test the impact of the decolonization protocol: the *education* group (n=1,063) received and reviewed an educational binder on hygiene, cleanliness, and MRSA transmission; the *decolonization* group (n=1,058) received this information and as well underwent decolonization for five days twice monthly for six months, with the protocol consisting of nasal mupirocin and chlorhexidine body and mouth wash. During these six months, swabs were collected from the participants at discharge from the hospitalization and three follow-up visits, which approximately took place at months 1, 3, and 6 after the discharge. Samples were taken from the nares, skin (axilla/groin), throat, and, if present, any wound. Participants had different numbers of observations because of trial exit or skipped visits. Overall, 20,506 samples were received in the first 6 months of follow up. Additionally, the dataset includes reported adherence to the protocol, enabling the assessment of both real-world uptake and consideration of the ideal scenario of full compliance.

To model the interactions of MRSA colonization between body sites, we used an extension of a Hidden Markov Model (HMM, a well-known machine learning algorithm[10]), coupled HMM (cHMM)[11, 12], where the probability of colonization at a given site depends on the colonization of all body sites in the previous time step[13, 14] (Figure 1). We considered two formulations: *Additive-cHMM* and *Or-cHMM*. In the former, the probability of colonization at a given site is an additive function of colonization at the other sites. In the latter formulation, other body sites are considered jointly, and the probability of colonization in the body site of interest increases if *any* of the other sites is colonized. According to a predictive model selection criterion (leave-one-out cross-validation, LOO-CV)[15], the *Additive-cHMM* had superior accuracy, and consequently, results for this model are discussed in the main text (see Supplementary Materials for additional results and details). We provide an R-package that implements these models with an efficient collapsed Metropolis-within-Gibbs Markov Chain Monte Carlo (MCMC) algorithm, which yields Bayesian credible intervals (CI) for all model parameters (see Methods).

## 2 Results

### Nasal carriage of MRSA is the most persistent and most affected by the decolonization protocol

For each body site, we estimated the persistence of colonization by calculating the probability that the site would be colonized in the next time step when it was colonized in the previous time step, assuming other sites were not colonized. We found that in the education group, the estimated persistence was highest in the nares (Figure 2a). In the decolonization group, application of the protocol decreased persistence at each site, with the largest decrease in the nares[16] (Figure 2a). We also estimated the relapse probability, i.e., the probability of a site getting colonized when it was not colonized in the previous time step, assuming other sites were not colonized. However, this did not seem to be significantly affected by the decolonization protocol (Figure 2b). The number of samples from wounds was relatively small in the analysis, which caused the estimated parameters for the wound to have considerably wider CIs than estimates for the other sites. Repeating the analysis only on patients who reported adherence to the protocol led to almost identical results since most (70%) of the data were collected from adherent patients (Supplementary Materials).

### The decolonization protocol reduces the transmission of MRSA between body sites

We quantified the strength of dependencies of MRSA colonization between body sites using the increase in the log-odds of colonization in the next time step caused by colonized sites in the previous time step (Figure 3). First, as expected, the most important determinant of colonization at any site was whether the site itself was colonized in the previous time step. Second, the nares was the only site linked to all the other sites in the education group, highlighting it as the likely source from which MRSA spreads to the other sites. Third, the decolonization protocol decreased or even removed the impact the nares had on any other site (Figure 3). In contrast, other dependencies between body sites were not significantly affected by the treatment. These findings indicate two mechanisms responsible for the efficiency of the decolonization protocol: 1) the protocol decreased persistence in the nares, which is a central hub for colonization (Figure 2a), and 2) the protocol weakened links between body sites, reducing self-inoculation (Figure 3).

### The cHMM accurately predicts the decrease in MRSA carriage over the study period

The key metric of success for decolonization is the extent to which it reduces MRSA carriage. To examine the cHMM predictive accuracy, we simulated patient trajectories using parameters from the estimated models and compared the decrease in colonization in our model-based simulations to the observed decrease from the CLEAR trial data. The predicted reduction of the carriage was from 67% to 48% (90% CI is 45%-51%) in the education group and down to 27% (24%-29%) in the decolonization group (Figure 4). These predictions accurately match the reductions observed in the trial. However, a discrepancy between predictions and observations was seen in intermediate visits; especially in the decolonization group, the observed carriage rate decreases more rapidly than predicted in the beginning and then more slowly. This outcome may reflect heterogeneity in treatment response among the study population: some patients responded very fast to the treatment, others more slowly.

### Treatments on nares, skin, and wound account for 93% of the effect of full decolonization

The importance of the nares has previously been noted in data from an intensive care unit[17]. However, the added value of decolonization efforts that focus on body sites other than the nares is not well understood. To quantify the impact of decolonization efforts at each individual site, we predicted the efficiency of hypothetical simplified treatments, consisting of only a subset of the decolonization protocol (Figure 4). Nasal mupirocin decreased the estimated carriage from 67% to 34% (31%-37%), compared to 27% (24%-29%) of the full decolonization (Figure 4). Other single treatments were far less effective, reducing the carriage to 46% (43%-48%) with chlorhexidine mouthwash (CHG Oral) or 44% (41%-46%) with chlorhexidine bodywash (CHG Skin), a minor improvement compared to the 48% (45%-51%) of education alone. The best combination of two treatment components includes mupirocin and CHG Skin, i.e., leaving out CHG Oral, and this decreased carriage to 30% (27%-33%), which corresponds to 93% (85%-99%) of the impact of the full decolonization protocol (Figure 4). In contrast, the combination of CHG Skin + CHG Oral and mupirocin + CHG Skin decreased the carriage only to 41%(38%-44%) and 31%(29%-34%), respectively (Suplementary Materials).

### Enhancing clearance at the nares is predicted to achieve the largest gain in overall decolonization success

We used our model to predict the impact of future improvements to components of the decolonization protocol. We assumed optimal treatment components, defined as treatments that resulted in immediate clearance of a specific site, could be applied on top of the existing protocol. The results were consistent with the previous finding of the nares’ central role in colonization and its clearance. In particular, the estimated final carriage with an optimized clearance of MRSA nasal carriage was 13% (11%-16%), while optimizing the throat or skin/wound treatments decreased the carriage to 17% (15%-20%) or 18% (15%-20%), respectively (Figure 4). Optimizing treatments both on nares + throat has the potential to reduce carriage further down to 5% (4%-7%) compared to 6% (5%-7%) of the nares + skin/wound, and 10% (9%-12%) of the throat + skin/wound (Suplementary Materials).

## 3 Discussion

Our analysis made several simplifying assumptions. First, we only considered the presence or absence of MRSA, ignoring the type and the burden of MRSA[18]. However, heterogeneity in treatment responses among patients and strain type of MRSA might affect the colonization duration. Additionally, the analysis focused on individual MRSA clearance, but the decolonization protocol may have other benefits, i.e., reduction in transmission between patients[19, 20]. Second, our model assumed that patients came to the follow-up visits at the same time intervals, while in practice, there were some variations. Third, we assumed that whether an observation was missing did not depend on the colonization status of the body site, which is likely not always correct; for example, if a previously colonized wound is healed and not colonized anymore, no additional samples would be taken from the wound. Finally, our analysis is based on a Bayesian approach and required a specification of prior probabilities on model parameters to express what is believed about their values before seeing the data. We used relatively noninformative priors, and a comparison included in the Supplementary Materials shows that the model fit was not sensitive to the exact values of the priors. Despite these limitations, the model successfully estimated interactions between body sites, quantified the impact of individual treatment components, and predicted the decrease in MRSA carriage over the study period. Additionally, because of its Markovian property, it allows us to estimate the sensitivity and specificity of MRSA colonization detection as a byproduct (Suplementary Materials).

We presented a model for the dynamics and clearance of MRSA colonization. Our model estimated interactions between multiple body sites, disentangled the impact of individual treatment components on the decolonization protocol’s overall efficiency, and accurately predicted the decrease in carriage percentage over the study period. We found that not all treatment components contribute equally to successful decolonization. Among the components of the CLEAR Trial, mupirocin was estimated as the most effective treatment component, and our results indicated that omission of chlorhexidine mouthwash still yields a highly effective decolonization regimen. More broadly, this analysis indicates that the impact of individual interacting components of complex clinical protocols and RCT data can be rigorously deconvoluted and assessed *in silico* by machine learning tools like those used here, thereby enabling the design of interventions that are more efficient, easier to adhere to, and thus more likely to be successful.

## 4 Methods

### Isolate Collection

MRSA isolates were collected as part of the CLEAR (Changing Lives by Eradicating Antibiotic Resistance) Trial. The trial was designed to compare the impact of a repeated decolonization protocol plus education on general hygiene and environmental cleaning to education alone on MRSA infection and hospitalization [3]. Study subjects in the trial were recruited from hospitalized patients based on an MRSA positive culture or surveillance swabs. After recruitment, swabs were obtained from different body parts of subjects (nares, skin, throat, and wound) around the time of hospital discharge (ENRL) and at 1,3,6, and 9 months (V1-V4, respectively) following the initial visit. The decolonization lasted only for 6 months, and consequently in this study we only model visits until V3. Note that some enrolled study subjects, despite a history of MRSA colonization or infection, did not have swabs positive for MRSA at the first time point (ENRL). More details about sample collection and data preprocessing can be found in the [3].

### Hidden Markov Model

Markov models are a well-known tool for time series analysis. Hidden Markov Models (HMM) are a special case where Markov process is observed indirectly through noisy observations. Discrete-time discrete-state Hidden Markov model with a latent sequence ** π**,

*π*

_{t}∈ {1, …, K} and observations

**x**,

*x*

_{t}∈ {1, …, L} for

*t*∈ {1, …,

*τ*}, is defined as where K is number of latent states,

*τ*is the number of time steps, and

*L*the number of possible observed states. Parameters of the model are the initial state probability

*π*

_{0}, transition probability

*T*, and emission probability

*E*, defined as

Note that the rows of *T* and *E* are probability distributions by definition. We will jointly denote the model parameters as *θ* = {*π*_{0}, *T, E*}. If *θ* is known, then we can estimate the latent states *π* conditional on parameters *θ* and observations **x** using the forward-filtering backward-sampling algorithm.

### Forward-Filtering Backward-Sampling

Filtering is the estimation of the current hidden state by using all observations so far, *p* (*π*_{t} | *x*_{1:t}) [12]. This is known as *forward-filtering* and it can be represented as follows:

We define *α* (*π*_{t}) = *p* (*π*_{t} |*x*_{1:t}) and substitute that into Equation 7, which yields the following recursive equation known as *α-recursion* or *forward recursion*;

The recursion starts with *α* (*π*_{1}) = *p* (*π*_{1} | *x*_{1}) = *p* (*x*_{1} | *π*_{1}) *p* (*π*_{1}) */p* (*x*_{1}). This formula shows that the filtered distribution *α*(.) is propagated through to the next time-step, where it acts like a **‘prior’**. In other words, at each time step the calculated posterior becomes the new prior for the next time-step [21].

To estimate the posterior distribution, we sample from the joint distribution of the latent sequence, *p*(*π*_{1:τ} | *x*_{1:τ}):

According to Equation 10, to be able to compute the posterior distribution, we need **‘time-reversed’** transitions *p*(*π*_{t} | *π*_{t+1},, *x*_{1:t}), which are obtained using *α-recursion*s from the forward filtering. The simulation starts from the end of the sequence and proceeds backwards recursively. Starting point is

For each time step *t*, the unnormalized sampling probability can be calculated using the previously sampled latent state , and sampled from using

This procedure is known as forward-filtering backward-sampling[21].

### Learning Model Parameters

Conditionally on the latent states and observations the parameters *θ* = {*π*_{0}, *T, E*} of the model can be estimated. The posterior distribution is given by

Factorizing the likelihood term in Equation 14 as in Equation 1, we can get the following result:

Since, *π*_{0}, rows of the *T*, and rows of the *E* are probability distribution, [.]’s inside the Equation 16 corresponds to Multinomial distribution with parameters and , where corresponds to the initial state counts, *T* ^{counts} corresponds to the state to state transition counts, and *E*^{counts} corresponds to the state to observation emission counts. Since the Dirichlet is conjugate prior for the Multinomial distribution, we set a Dirichlet prior for *θ* = {*π*_{0}, *T, E*} with hyperparameters , and , respectively. Because of the conjugacy, the resulting posterior in Equation 14 is also a Dirichlet distribution:

### Inference of HMM with MCMC Sampling

There are several ways to estimate the parameters of HMM when both hidden states and model parameters *θ* are unknown, for example the expectation-maximization algorithm, variational Bayes, or MCMC sampling. In this paper, we will use MCMC. The goal of MCMC is to produce draws from the posterior *p*(** π**,

*θ*|

**x**). The Gibbs sampler is an MCMC algorithm suitable for high dimensional problems, such as the HMM, and it samples parameters one-by-one using their distributions conditional on the other parameters. In the HMM, it will alternate between sampling the model parameters

*θ*conditional on the latent sequence

**from**

*π**p*(

*θ*|

**), and the latent sequence**

*π*, x**conditional on the model parameters**

*π**θ*from

*p*(

**|**

*π**θ*,

**x**). The algorithm iterates the following steps

*N*times, which gives

*N*draws from the posterior:

Sample the model parameters from Equation 17, Equation 18, and Equation 19, respectively. (Note that they are independent of each other given latent sequence.)

Sample the latent sequence of HMM using the updated model parameters

*θ*using the Forward-FilteringBackward-Sampling algorithm.

Before calculating the posteriors, the first *N*_{0} samples are discarded as a convergence (warm-up) period of the Markov chain. The plate diagram of the HMM is given in the Figure 5a and the algorithm is described in Algorithm 1.

### Coupled Hidden Markov Model

In the ordinary HMM transition parameters do not change. However, in the CHMM, transitions in one chain are affected by other chains. In principle, it would be possible to define a single HMM where the latent state represents jointly the latent states of all individual chains, and adapt the solution for the HMM described above. However, this is inefficient where there are many chains because the number of states would grow as *O*(*K*^{C}), where *K* is the number hidden states in a chain and *C* the number of chains. On the high level, our key idea is that the transition matrix of each chain is modeled conditionally on the states of the other chains (and not as a single large joint transition matrix). A similar formulation is considered by [14], however the details are different.

Note that a transition matrix *T* for a specific chain will not be time-independent anymore, but changes at each time step depending on the states of the other chains. The graphical representation of the model (showing just two chains) is given in Figure 1, and the plate diagram of the CHMM is given in Figure 5b. In this formulation there are *O*(*KC*) parameters, which is much more efficient many chains.

Defining the transition matrix is a critical part of the algorithm, and it is essential to carry as much information about the other chains as possible. We model the dependencies between the chains with parameters ** β**. Assume there are

*C*chains and let C denote the set of chains. Further assume that each chain can be in one of

*K*possible states, and let K denote the set of states. Finally, assume that there is a baseline state such that if a chain is in that state, it does not affect other chains (in our application this state corresponds to the absence of MRSA colonization in the respective body site). The transition matrix for the chain

*ĉ*, at time

*t*, denoted as, is defined by;

The transition matrix is obtained in Equation 21 by applying a row-wise softmax operator *σ*_{row} to the unnormalized transition matrix . Parameter corresponds to an intercept matrix and it specifies the transition matrix of the target chain *ĉ* when all other chains are in the baseline state . Parameter represents the impact of chain *c* on the target chain *ĉ* and it is added to whenever chain *c* was in state in the previous time step. We will denote all the parameters for target chain *ĉ* as *β*^{[ĉ ←.]}.

In Equation 20, the unnormalized transition probability *U* is an additive function of the latent states of the other chains, hence we call it the *Additive-CHMM*. We also implemented formulation that was simple than the Additive-CHMM, which we call *Or-CHMM*. In the Or-CHMM, the target chain is not affected by the other chains individually; instead, there the target chain is affected whenever *any* of the other chains is in some state different from the baseline state, and the unnormalized transition matrix is given by:

In our application we have *K* = 2, corresponding to the presence and absence of colonization in a given body site. In this case, Equation 22 has a simple form: If all other chains are in state , the output is only *β*_{0}, and if any of the other chains is in a state other , then the output is the sum of *β*_{0} and *β*_{k}.

### Design and Interpretation of the β parameters

In the CHMM, the transition matrix *T* is modeled as a function of *β* parameters such that each *β*_{k} is a matrix of the same size as the transition matrix. Therefore, the *β* parameters are a list of matrices. Since the rows of a transition matrix *T* are assumed independent, we also model *β* parameters such that the rows of those matrices are independent. However, to avoid redundant parameters, the rows on each row of *β* are assumed to sum to zero, In practice, we sample the first *K*−1 parameter on a given row and set the *K*th element to equal the negative of the sum of all the other parameters. This constraint ensures that the mapping from *β* to *T* is one-to-one.

To give an interpretation to the parameters *β*, we start with Equation 21 and write the softmax for a single element of a transition matrix *T* (*i, j*):
from which it follows after straightforward algebra:

In our application *K* = 2 and sums of the rows of *β*_{k} are assumed equal to 0, so consequently also the rows in the *U*_{t} sum to zero, which gives us:
where ¬*j* refers to the other element on the row that is not *j*. Therefore, in this 2-dimensional case, all *β*_{0} values correspond to half of the log-odds of the respective transition probability, and similarly, parameters *β*_{k} representing the interactions between the chains correspond to half of the change in the log-odds because of the presence of colonization in the other chain.

### Adaptive Metropolis Within Gibbs Sampling

So far, we have described the design specifications of the CHMM. Since transitions are changing at each time step, it is neither feasible to calculate their sufficient statistics as in Equation 16 nor the full-conditionals required for Gibbs sampling in closed form. To resolve this, we sample the *β* using a Metropolis-Hastings step within the Gibbs sampler because conditional on the other chains, the transition parameters are fully determined by the *β* parameters. To get a draw from *p*(*β* | ** π, x, π**′,

*E*), where

**′ states of the other chains, we use the following steps;**

*π*Make a proposal

*β** ∼*q*(*β*)Transform

*β*samples to transition probabilities according to latent states of the other chains,*π*_{l}Accept the proposal

*β** with probability where

The first quantity on the right-hand side of the Equation 28 can be directly calculated with the forward algorithm described in the Forward-Filtering Backward-Sampling. We used the Gaussian proposal in which the previously sampled *β* is the mean such that *q*(*β** |*β*) = *N* (*β** | *β, c*^{2}Σ). As the variance parameter, we used a fixed diagonal covariance matrix during the warm-up period and then used the covariance matrix estimated from the previous samples. Among metropolis jumping rules, the most efficient *c* ≈ 5.76*/d*, where *d* is the dimension of the sampling [22]. Since this proposal distribution is symmetric (i.e., Normal), *q* (*β* | *β**) and *q* (*β**| *β*) in Equation 27, cancel out each other. Algorithm 2 summarizes the whole algorithm.

### Implementation details

We set the initial covariance matrix for the metropolis proposal as 0.01 × *I*, where *I* is the identity matrix, which corresponds to a step size giving the optimal acceptance rate of ≈ 23% [23]. We set the prior of *β*_{0} as *N* (*β*_{0} | 0, 5), which is almost uninformative so that the estimates are not affected strongly by the prior. For the rest of the *β* parameters, denoted by *β*_{k}, we used sparsity encouraging Horseshoe prior with mean and scale parameters are 0 and 1, respectively (see Suplementary Materials for details). We used a uniform prior on the initial state probabilities *π*_{0}, and very weak Dirichlet priors for the rows of emission probabilities *E* such that we set the value to 0.9 if the latent state emits the same observation and set the value to 0.1*/*(*n*−1), *n* is the number of states, if the latent state emits a different observation. In such a formulation, except for the initialization, the prior has a negligible effect because it is summed with observation counts during inference. We drew 15,000 MCMC samples, and we set the warm-up length as 7,500. Posterior probabilities are calculated using the remaining MCMC samples. HMM and cHMM implicitly assume that time intervals between observations are the same, which is not the case in our data. Therefore, during the model training, we assumed that there are missing observations at 2, 4, and 5 months after enrollment.

## Data Availability

The CLEAR (Changing Lives by Eradicating Antibiotic Resistance) Trial demonstrated that the use of a post-discharge decolonization protocol in MRSA carriers reduces infection and hospitalization rates; ClinicalTrial.gov number NCT01209234, see reference https://doi.org/10.1056/NEJMoa1716771 for informed consent and institutional review board approvals. The data used in the present article is available on request from the corresponding author Susan Huang of the https://doi.org/10.1056/NEJMoa1716771, a co-author of this study.

## Data Availability

The CLEAR (Changing Lives by Eradicating Antibiotic Resistance) Trial demonstrated that the use of a post-discharge decolonization protocol in MRSA carriers reduces infection and hospitalization rates; ClinicalTrial.gov number NCT01209234, see reference [9] for informed consent and institutional review board approvals. The data set used in the present article is available on request from the corresponding author Susan Huang of the [9], a coauthor of this study.

## Code Availability

The R package is available at https://github.com/onurpoyraz/chmmMCMC.

## Author Contributions

YHG and PM perceived the study. OP, PM, and YHG designed the model. OP implemented the model and run the analyses. OP, SSH, YHG, and PM interpreted the results. MRAS, LGM, JAM, and SSH prepared the data. OP, PM, YHG wrote the manuscript with comments from all authors.

## Supplementary Materials

### Model Comparison

We compared three models: the single-HMM as a baseline (i.e., modeling each site independently with an HMM), the Or-cHMM, and the Additive-cHMM model using six different initial parameterizations. First, we considered prior *N* (*β*_{k} | 0, 5), as in *β*_{0}, and prior *N* (*β*_{k} | 0, 1), to see how sensitive the performance is to the strength of the prior. In addition to these Gaussian priors, we also considered the Laplace and Horseshoe distributions, which are sparse priors that encourage values close to zero. Note that *β*_{0} always had the same uninformative Gaussian prior because the *β*_{0} values represent the impact of the history of a chain on its future values, which is expected to be non-sparse. On the other hand, *β*_{k} values correspond to the effects of the other chains, which may be zero in practice.

We compared the models with three different information criteria (IC): 1) The Akaike Information Criteria (AIC), which is the summation of the negative log-likelihood of the model and the number of parameters, 2) the Watanabe-Akaike Information Criteria (WAIC), also known as the Widely applicable information criteria, which is the generalized version of the AIC onto singular statistical models [24] and 3) the leave-one-out cross-validation (LOO-CV), which measures model performance by estimating pointwise out-of-sample prediction accuracy from a fitted Bayesian model using the log-likelihood evaluated at the posterior simulations of the parameter values[15]. WAIC and LOO-CV are more computationally intensive than AIC. Still, they are theoretically preferred over AIC for Bayesian models[15]. For all criteria, a lower score is better.

## Footnotes

↵msater{at}hsph.harvard.edu

lmiller{at}lundquist.org, jmckinnell{at}lundquist.org

sshuang{at}hs.uci.edu

ygrad{at}hsph.harvard.edu

pekka.marttinen{at}aalto.fi