Interpretable classifiers using rules and Bayesian analysis: Building a better stroke prediction model

Introduction

Our goal is to build predictive models that are highly accurate, yet are highly interpretable. These predictive models will be in the form of sparse decision lists, which consist of a series of if…then… statements where the if statements define a partition of a set of features and the then statements correspond to the predicted outcome of interest. Because of this form, a decision list model naturally provides a reason for each prediction that it makes. Figure 1 presents an example decision list that we created using the Titanic data set available in R. This data set provides details about each passenger on the Titanic, including whether the passenger was an adult or child, male or female, and their class (1st, 2nd, 3rd or crew). The goal is to predict whether the passenger survived based on his or her features. The list provides an explanation for each prediction that is made. For example, we predict that a passenger is less likely to survive than not because he or she was in the 3rd class. The list in Figure 1 is one accurate and interpretable decision list for predicting survival on the Titanic, possibly one of many such lists. Our goal is to learn these lists from data.

Our model, called Bayesian Rule Lists (BRL), produces a posterior distribution over permutations of if…then… rules, starting from a large, pre-mined set of possible rules. The decision lists with high posterior probability tend to be both accurate and interpretable, where the interpretability comes from a hierarchical prior over permutations of rules. The prior favors concise decision lists that have a small number of total rules, where the rules have few terms in the left-hand side.

BRL provides a new type of balance between accuracy, interpretability and computation. Consider the challenge of constructing a predictive model that discretizes the input space in the same way as decision trees [Breiman et al. (1984); Quinlan (1993)], decision lists [Rivest (1987)] or associative classifiers [Liu, Hsu and Ma (1998)]. Greedy construction methods like classification and regression trees (CART) or C5.0 are not particularly computationally demanding, but, in practice, the greediness heavily affects the quality of the solution, both in terms of accuracy and interpretability. At the same time, optimizing a decision tree over the full space of all possible splits is not a tractable problem. BRL strikes a balance between these extremes, in that its solutions are not constructed in a greedy way involving splitting and pruning, yet it can solve problems at the scale required to have an impact in real problems in science or society, including modern healthcare.

A major source of BRL’s practical feasibility is the fact that it uses pre-mined rules, which reduces the model space to that of permutations of rules as opposed to all possible sets of splits. The complexity of the problem then depends on the number of pre-mined rules rather than on the full space of feature combinations; in a sense, this algorithm scales with the sparsity of the data set rather than the number of features. As long as the pre-mined set of rules is sufficiently expressive, an accurate decision list can be found and, in fact, the smaller model space might improve generalization [through the lens of statistical learning theory, Vapnik (1995)]. An additional advantage to using pre-mined rules is that each rule is independently both interpretable and informative about the data.

BRL’s prior structure encourages decision lists that are sparse. Sparse decision lists serve the purpose of not only producing a more interpretable model, but also reducing computation, as most of the sampling iterations take place within a small set of permutations corresponding to the sparse decision lists. In practice, BRL is able to compute predictive models with accuracy comparable to state-of-the-art machine learning methods, yet maintain the same level of interpretability as medical scoring systems.

The motivation for our work lies in developing interpretable patient-level predictive models using massive observational medical data. To this end, we use BRL to construct an alternative to the CHADS2 score of Gage et al. (2001). CHADS2 is widely used in medical practice to predict stroke in patients with atrial fibrillation. A patient’s CHADS2 score is computed by assigning one “point” each for the presence of congestive heart failure (C), hypertension (H), age 75 years or older (A) and diabetes mellitus (D), and by assigning 2 points for history of stroke, transient ischemic attack or thromoembolism (S2). The CHADS2 score considers only 5 factors, whereas the updated CHA2DS2-VASc score [Lip et al. (2010b)] includes three additional risk factors: vascular disease (V), age 65 to 74 years old (A) and female gender (Sc). Higher scores correspond to increased risk. In the study defining the CHADS2 score [Gage et al. (2001)], the score was calibrated with stroke risks using a database of 1733 Medicare beneficiaries followed for, on average, about a year.

Bayesian rule lists

In Sections 2.1 and 2.2 we provide the association rule concepts and notation upon which the method is built. Section 2.3 introduces BRL by outlining the generative model. Sections 2.4 and 2.5 provide detailed descriptions of the prior and likelihood, and then Sections 2.6 and 2.7 describe sampling and posterior predictive distributions.

An association rule $a\rightarrow b$ is an implication with an antecedent $a$ and a consequent $b$. For the purposes of classification, the antecedent is an assertion about the feature vector $x_{i}$ that is either true or false, for example, “$x_{i,1}=1\mbox{ and }x_{i,2}=0$.” This antecedent contains two conditions, which we call the cardinality of the antecedent. The consequent $b$ would typically be a predicted label $y$. A Bayesian association rule has a multinomial distribution over labels as its consequent rather than a single label:

The multinomial probability is then given a prior, leading to a prior consequent distribution:

Given observations $(\mathbf{x},\mathbf{y})$ classified by this rule, we let $N_{\cdot,l}$ be the number of observations with label $y_{i}=l$, and $N=(N_{\cdot,1},\ldots,N_{\cdot,L})$. We then obtain a posterior consequent distribution:

The core of a Bayesian decision list is an ordered antecedent list $d=(a_{1},\ldots,a_{m})$. Let $N_{j,l}$ be the number of observations $x_{i}$ that satisfy $a_{j}$ but not any of $a_{1},\ldots,a_{j-1}$, and that have label $y_{i}=l$. This is the number of observations to be classified by antecedent $a_{j}$ that have label $l$. Let $N_{0,l}$ be the number of observations that do not satisfy any of $a_{1},\ldots,a_{m}$ and that have label $l$. Let $\mathbf{N}_{j}=(N_{j,1},\ldots,N_{j,L})$ and $\mathbf{N}=(\mathbf{N}_{0},\ldots,\mathbf{N}_{m})$.

A Bayesian decision list $D=(d,\bolds{\alpha},\mathbf{N})$ is an ordered list of antecedents together with their posterior consequent distributions. The posterior consequent distributions are obtained by excluding data that have satisfied an earlier antecedent in the list. A Bayesian decision list then takes the form:

if $a_{1}$ then $y\sim\operatorname{Multinomial}(\bolds{\theta}_{1})$, $\bolds{\theta}_{1}\sim\operatorname{Dirichlet}(\bolds{\alpha}+\mathbf{N}_{1})$

else if $a_{2}$ then $y\sim\operatorname{Multinomial}(\bolds{\theta}_{2})$, $\bolds{\theta}_{2}\sim\operatorname{Dirichlet}(\bolds{\alpha}+\mathbf{N}_{2})$

else if $a_{m}$ then $y\sim\operatorname{Multinomial}(\bolds{\theta}_{m})$, $\bolds{\theta}_{m}\sim\operatorname{Dirichlet}(\bolds{\alpha}+\mathbf{N}_{m})$

else $y\sim\operatorname{Multinomial}(\bolds{\theta}_{0})$, $\bolds{\theta}_{0}\sim\operatorname{Dirichlet}(\bolds{\alpha}+\mathbf{N}_{0})$.

Any observations that do not satisfy any of the antecedents in $d$ are classified using the parameter $\theta_{0}$, which we call the default rule parameter.

2 Antecedent mining

We are interested in forming Bayesian decision lists whose antecedents are a subset of a preselected collection of antecedents. For data with binary or categorical features this can be done using frequent itemset mining, where itemsets are used as antecedents. In our experiments, the features were binary and we used the FP-Growth algorithm [Borgelt (2005)] for antecedent mining, which finds all itemsets that satisfy constraints on minimum support and maximum cardinality. This means each antecedent applies to a sufficiently large amount of data and does not have too many conditions. For binary or categorical features the particular choice of the itemset mining algorithm is unimportant, as the output is an exhaustive list of all itemsets satisfying the constraints. Other algorithms, such as Apriori or Eclat [Agrawal and Srikant (1994); Zaki (2000)], would return an identical set of antecedents as FP-Growth if given the same minimum support and maximum cardinality constraints. Because the goal is to obtain decision lists with few rules and few conditions per rule, we need not include any itemsets that apply only to a small number of observations or have a large number of conditions. Thus, frequent itemset mining allows us to significantly reduce the size of the feature space, compared to considering all possible combinations of features.

The frequent itemset mining that we do in our experiments produces only antecedents with sets of features, such as “diabetes and heart disease.” Other techniques could be used for mining antecedents with negation, such as “not diabetes” [Wu, Zhang and Zhang (2004)]. For data with continuous features, a variety of procedures exist for antecedent mining [Fayyad and Irani (1993); Dougherty, Kohavi and Sahami (1995); Srikant and Agrawal (1996)]. Alternatively, one can create categorical features using interpretable thresholds (e.g., ages 40–49, 50–59, etc.) or interpretable quantiles (e.g., quartiles)—we took this approach in our experiments.

We let $\mathcal{A}$ represent the complete, pre-mined collection of antecedents, and suppose that $\mathcal{A}$ contains $|\mathcal{A}|$ antecedents with up to $C$ conditions in each antecedent.

3 Generative model

We now sketch the generative model for the labels $\mathbf{y}$ from the observations $\mathbf{x}$ and antecedents $\mathcal{A}$. Define $a_{<j}$ as the antecedents before $j$ in the rule list if there are any, for example, $a_{<3}=\{a_{1},a_{2}\}$. Similarly, let $c_{j}$ be the cardinality of antecedent $a_{j}$, and $c_{<j}$ the cardinalities of the antecedents before $j$ in the rule list. The generative model is then:

Sample a decision list length $m\sim p(m|\lambda)$.

Sample the default rule parameter $\theta_{0}\sim\operatorname{Dirichlet}(\alpha)$.

Sample the cardinality of antecedent $a_{j}$ in $d$ as $c_{j}\sim p(c_{j}|c_{<j},\mathcal{A},\eta)$.

Sample $a_{j}$ of cardinality $c_{j}$ from $p(a_{j}|a_{<j},c_{j},\mathcal{A})$.

Sample rule consequent parameter $\theta_{j}\sim\operatorname{Dirichlet}(\alpha)$.

Find the antecedent $a_{j}$ in $d$ that is the first that applies to $x_{i}$.

If no antecedents in $d$ apply, set $j=0$.

Sample $y_{i}\sim\operatorname{Multinomial}(\theta_{j})$.

Our goal is to sample from the posterior distribution over antecedent lists:

Given $d$, we can compute the posterior consequent distributions required to construct a Bayesian decision list as in Section 2.1. Three prior hyperparameters must be specified by the user: $\bolds{\alpha}$, $\lambda$ and $\eta$. We will see in Sections 2.4 and 2.5 that these hyperparameters have natural interpretations that suggest the values to which they should be set.

4 The hierarchical prior for antecedent lists

Suppose the list of antecedents $d$ has length $m$ and antecedent cardinalities $c_{1},\ldots,c_{m}$. The prior probability of $d$ is defined hierarchically as

We take the distributions for list length $m$ and antecedent cardinality $c_{j}$ to be Poisson with parameters $\lambda$ and $\eta$, respectively, with proper truncation to account for the finite number of antecedents in $\mathcal{A}$. Specifically, the distribution of $m$ is Poisson truncated at the total number of preselected antecedents:

The distribution of $c_{j}$ must be truncated at zero and at the maximum antecedent cardinality $C$. Additionally, any cardinalities that have been exhausted by point $j$ in the decision list sampling must be excluded. Let $R_{j}(c_{1},\ldots,c_{j},\mathcal{A})$ be the set of antecedent cardinalities that are available after drawing antecedent $j$. For example, if $\mathcal{A}$ contains antecedents of size $1$, $2$ and $4$, then we begin with $R_{0}(\mathcal{A})=\{1,2,4\}$. If $\mathcal{A}$ contains only $2$ rules of size $4$ and $c_{1}=c_{2}=4$, then $R_{2}(c_{1},c_{2},\mathcal{A})=\{1,2\}$ as antecedents of size $4$ have been exhausted. We now take $p(c_{j}|c_{<j},\mathcal{A},\eta)$ as Poisson truncated to remove values for which no rules are available with that cardinality:

If the number of rules of different sizes is large compared to $\lambda$, and $\eta$ is small compared to $C$, the prior expected average antecedent cardinality is close to $\eta$. Thus, $\eta$ can be set to the prior belief of the antecedent cardinality required to model the data.

Once the antecedent cardinality $c_{j}$ has been selected, the antecedent $a_{j}$ must be sampled from all available antecedents in $\mathcal{A}$ of size $c_{j}$. Here, we use a uniform distribution over antecedents in $\mathcal{A}$ of size $c_{j}$, excluding those in $a_{<j}$:

It is straightforward to sample an ordered antecedent list $d$ from the prior by following the generative model, using the provided distributions.

5 The likelihood function

The likelihood function follows directly from the generative model. Let $\bolds{\theta}=(\theta_{0},\theta_{1},\ldots,\theta_{m})$ be the consequent parameters corresponding to each antecedent in $d$, together with the default rule parameter $\theta_{0}$. Then, the likelihood is the product of the multinomial probability mass functions for the observed label counts at each rule:

We can marginalize over $\theta_{j}$ in each multinomial distribution in the above product, obtaining, through the standard derivation of the Dirichlet-multinomial distribution,

The prior hyperparameter $\bolds{\alpha}$ has the usual Bayesian interpretation of pseudocounts. In our experiments, we set $\alpha_{l}=1$ for all $l$, producing a uniform prior. Other approaches for setting prior hyperparameters such as empirical Bayes are also applicable.

6 Markov chain Monte Carlo sampling

We do Metropolis–Hastings sampling of $d$, generating the proposed $d^{*}$ from the current $d^{t}$ using one of three options: (1) Move an antecedent in $d^{t}$ to a different position in the list. (2) Add an antecedent from $\mathcal{A}$ that is not currently in $d^{t}$ into the list. (3) Remove an antecedent from $d^{t}$. Which antecedents to adjust and their new positions are chosen uniformly at random at each step. The option to move, add or remove is also chosen uniformly. The probabilities for the proposal distribution $Q(d^{*}|d^{t})$ depend on the size of the antecedent list, the number of pre-mined antecedents, and whether the proposal is a move, addition or removal. For the uniform distribution that we used, the proposal probabilities for a $d^{*}$ produced by one of the three proposal types is

To explain these probabilities, if there is a move proposal, we consider the number of possible antecedents to move and the number of possible positions for it; if there is an add proposal, we consider the number of possible antecedents to add to the list and the number of positions to place a new antecedent; for remove proposals we consider the number of possible antecedents to remove. This sampling algorithm is related to those used for Bayesian Decision Tree models [Chipman, George and McCulloch (1998; 2002), Wu, Tjelmeland and West (2007)] and to methods for exploring tree spaces [Madigan, Mittal and Roberts (2011)].

For every MCMC run, we ran 3 chains, each initialized independently from a random sample from the prior. We discarded the first half of simulations as burn-in, and then assessed chain convergence using the Gelman–Rubin convergence diagnostic applied to the log posterior density [Gelman and Rubin (1992)]. We considered chains to have converged when the diagnostic $\hat{R}<1.05$.

7 The posterior predictive distribution and point estimates

Additionally, (3) allows for the estimation of 95% credible intervals using the Dirichlet distribution function.

The posterior mean is often a good choice for a point estimate, but the interpretation of “mean” here is not clear since the posterior is a distribution over antecedent lists. We thus look for an antecedent list whose statistics are similar to the posterior mean statistics. Specifically, we are interested in finding a point estimate $\hat{d}$ whose length $m$ and whose average antecedent cardinality $\bar{c}=\frac{1}{m}\sum_{j=1}^{m}c_{j}$ are close to the posterior mean list length and average cardinality. Let $\bar{m}$ be the posterior mean decision list length and $\bar{\bar{c}}$ the posterior mean average antecedent cardinality, as estimated from the MCMC samples. Then, we choose our point estimate $\hat{d}$ as the list with the highest posterior probability among all samples with $m\in\{\lfloor{\bar{m}}\rfloor,\lceil{\bar{m}}\rceil\}$ and $\bar{c}\in[\lfloor{\bar{\bar{c}}}\rfloor,\lceil{\bar{\bar{c}}}\rceil]$. We call this point estimate BRL-point.

Another possible point estimate is the decision list with the highest posterior probability—the maximum a posteriori estimate. Given two list lengths, there are many more possible lists of the longer length than of the shorter length, so prior probabilities in (1) are generally higher for shorter lists. The maximum a posteriori estimate might yield a list that is much shorter than the posterior mean decision list length, so we prefer the BRL-point.

In addition to point estimates, we can use the entire posterior $p(d|\mathbf{x},\mathbf{y},\mathcal{A},\break\alpha,\lambda,\eta)$ to estimate $y$. The posterior predictive distribution for $y$ is

where $\mathbf{D}$ is the set of all ordered subsets of $\mathcal{A}$. The posterior samples obtained by MCMC simulation, after burn-in, can be used to approximate this sum. We call the classifier that uses the full collection of posterior samples BRL-post. Using the entire posterior distribution to make a prediction means the classifier is no longer interpretable. One could, however, use the posterior predictive distribution to classify, and then provide several point estimates from the posterior to the user as example explanations for the prediction.

Simulation studies

We use simulation studies and a deterministic data set to show that when data are generated by a decision list model, the BRL (Bayesian Rule Lists; see Section 1) method is able to recover the true decision list.

Given observations with arbitrary features and a collection of rules on those features, we can construct a binary matrix where the rows represent observations and the columns represent rules, and the entry is $1$ if the rule applies to that observation and otherwise. We need only simulate this binary matrix to represent the observations without losing generality. For our simulations, we generated independent binary rule sets with $100$ rules by setting each feature value to $1$ independently with probability $1/2$.

We generated a random decision list of size $5$ by selecting 5 rules at random, and adding the default rule. Each rule in the decision list was assigned a consequent distribution over labels using a random draw from the $\operatorname{Beta}(1/2,1/2)$ distribution, which ensures that the rules are informative about labels. Labels were then assigned to each observation using the decision list: For each observation, the label was taken as a draw from the label distribution corresponding to the first rule that applied to that observation.

For each number of observations $N\in\{100,250,500,1000,2500,5000\}$, we generated $100$ independent data sets $(\mathbf{x},\mathbf{y})$, for a total of $600$ simulated data sets. We did MCMC sampling with three chains as described in Section 2 for each data set. For all data sets, 20,000 samples were sufficient for the chains to converge.

To appropriately visualize the posterior distribution, we binned the posterior antecedent lists according to their distance from the true antecedent list, using the Levenshtein string edit distance [Levenshtein (1965)] to measure the distance between two antecedent lists. This metric measures the minimum number of antecedent substitutions, additions or removals to transform one decision list into the other. The results of the simulations are given in Figure 2.

Figure 2(a) shows that as the number of observations increases, the posterior mass concentrates on the true decision list. Figure 2(b) illustrates this concentration with two choices of the distribution of posterior distances to the true decision list, for $n$ small and for $n$ large.

2 A deterministic problem

Stroke prediction

We used Bayesian Rule Lists to derive a stroke prediction model using the MarketScan Medicaid Multi-State Database (MDCD). MDCD contains administrative claims data for 11.1 million Medicaid enrollees from multiple states. This database forms part of the suite of databases from the Innovation in Medical Evidence Development and Surveillance (IMEDS, http://imeds.reaganudall.org/) program that have beenmapped to a common data model [Stang et al. (2010)].

We extracted every patient in the MDCD database with a diagnosis of atrial fibrillation, one year of observation time prior to the diagnosis and one year of observation time following the diagnosis ($n=12\mbox{,}586$). Of these, 1786 (14%) had a stroke within a year of the atrial fibrillation diagnosis.

As candidate predictors, we considered all drugs and all conditions. Specifically, for every drug and condition, we created a binary predictor variable indicating the presence or absence of the drug or condition in the full longitudinal record prior to the atrial fibrillation diagnosis. These totaled 4146 unique medications and conditions. We included features for age and gender. Specifically, we used the natural values of 50, 60, 70 and 80 years of age as split points, and for each split point introduced a pair of binary variables indicating if age was less than or greater than the split point. Considering both patients and features, here we apply our method to a data set that is over 6000 times larger than that originally used to develop the CHADS2 score (which had $n=1733$ and considered 5 features).

We did five folds of cross-validation. For each fold, we pre-mined the collection of possible antecedents using frequent itemset mining with a minimum support threshold of $10\%$ and a maximum cardinality of $2$. The total number of antecedents used ranged from $2162$ to $2240$ across the folds. We set the antecedent list prior hyperparameters $\lambda$ and $\eta$ to $3$ and $1$, respectively, to obtain a Bayesian decision list of similar complexity to the CHADS2 score. For each fold, we evaluated the performance of the BRL point estimate by constructing a receiver operating characteristic (ROC) curve and measuring area under the curve (AUC) for each fold.

In Figure 3 we show the BRL point estimate recovered from one of the folds. The list indicates that past history of stroke reveals a lot about the vulnerability toward future stroke. In particular, the first half of the decision list focuses on a history of stroke, in order of severity. Hemiplegia, the paralysis of an entire side of the body, is often a result of a severe stroke or brain injury. Cerebrovascular disorder indicates a prior stroke, and transient ischaemic attacks are generally referred to as “mini-strokes.” The second half of the decision list includes age factors and vascular disease, which are known risk factors and are included in the CHA2DS2-VASc score. The BRL-point lists that we obtained in the 5 folds of cross-validation were all of length 7, a similar complexity to the CHADS2 and CHA2DS2-VASc scores which use 5 and 8 features, respectively.

The point estimate lists for all five of the folds of cross-validation are given in the supplemental material [Letham et al. (2015)]. There is significant overlap in the antecedents in the point estimates across the folds. This suggests that the model may be more stable in practice than decision trees, which are notorious for producing entirely different models after small changes to the training set [Breiman (1996a; 1996b)].

In Figure 4 we give ROC curves for all 5 folds for BRL-point, CHADS2 and CHA2DS2-VASc, and in Table 2 we report mean AUC across the folds. These results show that with complexity and interpretability similar to CHADS2, the BRL point estimate decision lists performed significantly better at stroke prediction than both CHADS2 and CHA2DS2-VASc. Interestingly, we also found that CHADS2 outperformed CHA2DS2-VASc despite CHA2DS2-VASc being an extension of CHADS2. This is likely because the model for the CHA2DS2-VASc score, in which risk factors are added linearly, is a poor model of actual stroke risk. For instance, the stroke risks estimated by CHA2DS2-VASc are not a monotonic function of score. Within the original CHA2DS2-VASc calibration study, Lip et al. (2010a) estimate a stroke risk of 9.6% with a CHA2DS2-VASc score of 7, and a 6.7% risk with a score of 8. The indication that more stroke risk factors can correspond to a lower stroke risk suggests that the CHA2DS2-VASc model may be misspecified, and highlights the difficulty in constructing these interpretable models manually.

All of the methods were applied to the data on the same, single Amazon Web Services virtual core with a processor speed of approximately 2.5 GHz and 4 GB of memory. Bayesian CART was unable to fit the data since it ran out of memory, and so it is not included in Table 2.

The BRL MCMC chains were simulated until convergence, which required 50,000 iterations for 4 of the 5 folds, and 100,000 for the fifth. The three chains for each fold were simulated in serial, and the total CPU time required per fold is given in Table 2, together with the CPU times required for training the comparison algorithms on the same processor. Table 2 shows that the BRL MCMC simulation was more than ten times faster than training SVM, and more than thirty times faster than training random forests, using standard implementations of these methods as described in the \hyperref[app]Appendix.

We further investigated the properties and performance of the BRL by applying it to two subsets of the data, female patients only and male patients only. The female data set contained 8368 observations, and the number of pre-mined antecedents in each of 5 folds ranged from 1982 to 2197. The male data set contained 4218 observations, and the number of pre-mined antecedents in each of 5 folds ranged from 1629 to 1709. BRL MCMC simulations and comparison algorithm training were done on the same processor as the full experiment. The AUC and training time across five folds for each of the data sets is given in Table 3.

The BRL point estimate again outperformed the other interpretable models (CHADS2, CHA2DS2-VASc, CART and C5.0), and the BRL-post performance matched that of random forests for the best performing method. As before, BRL MCMC simulation required significantly less time than SVM or random forests training. Point estimate lists for these additional experiments are given in the supplemental materials [Letham et al. (2015)].

Related work and discussion

Most widely used medical scoring systems are designed to be interpretable, but are not necessarily optimized for accuracy, and generally are derived from few factors. The Thrombolysis In Myocardial Infarction (TIMI) Score [Antman et al. (2000)], Apache II score for infant mortality in the ICU [Knaus et al. (1985)], the CURB-65 score for predicting mortality in community-acquired pneumonia [Lim et al. (2003)] and the CHADS2 score [Gage et al. (2001)] are examples of interpretable predictive models that are very widely used. Each of these scoring systems involves very few calculations and could be computed by hand during a doctor’s visit. In the construction of each of these models, heuristics were used to design the features and coefficients for the model; none of these models was fully learned from data.

In contrast with these hand-designed interpretable medical scoring systems, recent advances in the collection and storing of medical data present unprecedented opportunities to develop powerful models that can predict a wide variety of outcomes [Shmueli (2010)]. The front-end user interface of medical risk assessment tools are increasingly available online (e.g., http://www.r-calc.com). At the end of the assessment, a patient may be told he or she has a high risk for a particular outcome but without understanding why the predicted risk is high, particularly if many pieces of information were used to make the prediction.

In general, humans can handle only a handful of cognitive entities at once [Miller (1956); Jennings, Amabile and Ross (1982)]. It has long since been hypothesized that simple models predict well, both in the machine learning literature [Holte (1993)] and in the psychology literature [Dawes (1979)]. The related concepts of explanation and comprehensibility in statistical modeling have been explored in many past works [Bratko (1997); Madigan, Mosurski and Almond (1997); Giraud-Carrier (1998); Rüping (2006); Huysmans et al. (2011); Vellido, Martín-Guerrero and Lisboa (2012); Freitas (2014), e.g.].

Decision lists have the same form as models used in the expert systems literature from the 1970s and 1980s [Leondes (2002)], which were among the first successful types of artificial intelligence. The knowledge base of an expert system is composed of natural language statements that are if…then… rules. Decision lists are a type of associative classifier, meaning that the list is formed from association rules. In the past, associative classifiers have been constructed from heuristic greedy sorting mechanisms [Rivest (1987); Liu, Hsu and Ma (1998); Marchand and Sokolova (2005); Rudin, Letham and Madigan (2013)]. Some of these sorting mechanisms work provably well in special cases, for instance, when the decision problem is easy and the classes are easy to separate, but are not optimized to handle more general problems. Sometimes associative classifiers are formed by averaging several rules together, or having the rules each vote on the label and then combining the votes, but the resulting classifier is not generally interpretable [Li, Han and Pei (2001); Yin and Han (2003); Friedman and Popescu (2008); Meinshausen (2010)].

In a previous paper we proved that the VC dimension of decision list classifiers equals $|\mathcal{A}|$, the number of antecedents used to learn the model [Theorem 3, Rudin, Letham and Madigan (2013)]. This result leads to a uniform generalization bound for decision lists [Corollary 4, Rudin, Letham and Madigan (2013)]. This is the same as the VC dimension obtained by using the antecedents as features in a linear model, thus we have the same prediction guarantees. We then expect similar generalization behavior for decision lists and weighted linear combination models.

BRL interacts with the feature space only through the collection of antecedents $\mathcal{A}$. The computational effort scales with the number of antecedents, not the number of features, meaning there will generally be less computation when the data are sparse. This means that BRL tends to scale with the sparsity of the data rather than the number of features.

Decision trees are closely related to decision lists, and are in some sense equivalent: any decision tree can be expressed as a decision list, and any decision list is a one-sided decision tree. Decision trees are almost always constructed greedily from the top down, and then pruned heuristically upward and cross-validated to ensure accuracy. Because the trees are not fully optimized, if the top of the decision tree happened to have been chosen badly at the start of the procedure, it could cause problems with both accuracy and interpretability. Bayesian decision trees [Chipman, George and McCulloch (1998; 2002), Denison, Mallick and Smith (1998)] use Markov chain Monte Carlo (MCMC) to sample from a posterior distribution over trees. Since they were first proposed, several improvements and extensions have been made in both sampling methods and model structure [Wu, Tjelmeland and West (2007); Chipman, George and McCulloch (2010); Taddy, Gramacy and Polson (2011)]. The space of decision lists using pre-mined rules is significantly smaller than the space of decision trees, making it substantially easier to obtain MCMC convergence and to avoid the pitfalls of local optima. Moreover, rule mining allows for the rules to be individually powerful. Constructing a single decision tree is extremely fast, but sampling over the space of decision trees is extremely difficult (unless one is satisfied with local maxima). To contrast this with our approach, the rule mining step is extremely fast, yet sampling over the space of decision lists is very practical.

There is a subfield of artificial intelligence, Inductive Logic Programming [Muggleton and De Raedt (1994)], whose goal is to mine individual conjunctive rules. It is possible to replace the frequent itemset miner with an inductive logic programming technique, but this generally leads to losses in predictive accuracy; ideally, we would use a large number of diverse rules as antecedents, rather than a few (highly overlapping) complex rules as would be produced by an ILP algorithm. In our experiments to a follow-up work [Wang and Rudin (2015)], the use of an ILP algorithm resulted in a substantial loss in performance.

Interpretable models are generally not unique (stable), in the sense that there may be many equally good models, and it is not clear in advance which one will be returned by the algorithm. For most problems, the space of high quality predictive models is fairly large [called the “Rashomon Effect” Breiman (2001b)], so we cannot expect uniqueness. In practice, as we showed, the rule lists across test folds were very similar, but if one desires stability to small perturbations in the data generally, we recommend using the full posterior rather than a point estimate. The fact that many high performing rule lists exist can be helpful, since it means the user has many choices of which model to use.

This work is related to the Hierarchical Association Rule Model (HARM), a Bayesian model that uses rules [McCormick, Rudin and Madigan (2012)]. HARM estimates the conditional probabilities of each rule jointly in a conservative way. Each rule acts as a separate predictive model, so HARM does not explicitly aim to learn an ordering of rules.

There are related works on learning decision lists from an optimization perspective. In particular, the work of Rudin and Ertekin (2015) uses mixed-integer programming to build a rule list out of association rules, which has guarantees on optimality of the solution. Similarly to that work, Goh and Rudin (2014) fully learn sparse disjunctions of conjunctions using optimization methods.

There have been several follow-up works that directly extend and apply Bayesian Rule Lists. The work of Wang and Rudin (2015) on Falling Rule Lists provides a nontrivial extension to BRL whereby the probabilities for the rules are monotonically decreasing down the list. Wang et al. (2015) build disjunctions of conjunctive rules using a Bayesian framework similar to the one in this work. Zhang et al. (2015) have taken an interesting approach to constructing optimal treatment regimes using a BRL-like method, where, in addition to the criteria of accuracy, the rule list has a decision cost for evaluating it. It is possible to use BRL itself for that purpose as well, as one could give preference to particular antecedents that cost less. This sort of preference could be expressed in the antecedent prior distribution in (2). King, Lam and Roberts (2014) have taken a Bayesian Rule List approach to handle a challenging problem in text analysis, which is to build a keyword-based classifier that is easier to understand in order to solicit high quality human input. Souillard-Mandar et al. (2015) applied Bayesian Rule Lists and Falling Rule Lists to the problem of screening for cognitive disorders such as Alzheimer’s disease based on the digitized pen strokes of patients during the Clock Drawing test.

Shorter preliminary versions of this work are those of Letham et al. (2013; 2014). Letham et al. (2013) used a different prior and called the algorithm the Bayesian List Machine.

Conclusion

We are working under the hypothesis that many real data sets permit predictive models that can be surprisingly small. This was hypothesized over two decades decade ago [Holte (1993)]; however, we now are starting to have the computational tools to truly test this hypothesis. The BRL method introduced in this work aims to hit the “sweet spot” between predictive accuracy, interpretability and tractability.

Interpretable models have the benefits of being both concise and convincing. A small set of trustworthy rules can be the key to communicating with domain experts and to allowing machine learning algorithms to be more widely implemented and trusted. In practice, a preliminary interpretable model can help domain experts to troubleshoot the inner workings of a complex model, in order to make it more accurate and tailored to the domain. We demonstrated that interpretable models lend themselves to the domain of predictive medicine, and there is a much wider variety of domains in science, engineering and industry, where these models would be a natural choice.

Appendix

Acknowledgments

The authors thank Zachary Shahn and the OMOP team for help with the data.

Computer code \slink[doi]10.1214/15-AOAS848SUPPA \sdatatype.zip \sfilenameaoas848_suppa.zip \sdescriptionOur Python code used to fit decision lists to data, along with an example data set.

BRL point estimates \slink[doi]10.1214/15-AOAS848SUPPB \sdatatype.pdf \sfilenameaoas848_suppb.pdf \sdescriptionThe BRL point estimates for all of the cross-validation folds for the stroke prediction experiment, and BRL-point estimates for the female-only and male-only experiments.

References