<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/exploring-models-with-lime/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://nielsenmark.us/2018/11/09/exploring-models-with-lime/"> Programming on nielsenmark.us</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Recently at work I’ve been asked to help some clinicians understand why my risk model classifies specific patients as high risk. Just prior to this work I stumbled across the work of some data scientists at the University of Washington called <code>lime</code>. <a rel="nofollow" href="https://github.com/marcotcr/lime" target="_blank">LIME</a> stands for “Local Interpretable Model-Agnostic Explanations”. The idea is that I can answer those questions I’m getting from clinicians for a specific patient by locally fitting a linear (aka “interpretable”) model in the parameter space just around my data point. I decided to pursue <code>lime</code> as a solution and the last few months I’ve been focusing on implementing this explainer for my risk model. Happily, I also discovered an <a rel="nofollow" href="https://github.com/thomasp85/lime" target="_blank">R package</a> that implements this solution that originated in python.</p>
<div>
<h3>Sample Data</h3>
<p>So the first step to this blog was to find some public data for illustration. I remembered an example used in an <a rel="nofollow" href="http://www-bcf.usc.edu/~gareth/ISL/index.html" target="_blank">Introduction to Statistical Learning by James, Witten, Hastie and Tibshirani</a>. I will use the <code>Heart.csv</code> data which can be downloaded using the link in the code below:</p>
<pre><code>library(readr)
library(ranger)
library(tidyverse)
library(lime)

dat &lt;- read_csv(&quot;http://www-bcf.usc.edu/~gareth/ISL/Heart.csv&quot;)
dat$X1 &lt;- NULL</code></pre>
<p>Now let’s take a quick look at the data:</p>
<pre><code>Hmisc::describe(dat)</code></pre>
<pre><code>## dat 
## 
##  14  Variables      303  Observations
## ---------------------------------------------------------------------------
## Age 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      303        0       41    0.999    54.44     10.3       40       42 
##      .25      .50      .75      .90      .95 
##       48       56       61       66       68 
## 
## lowest : 29 34 35 37 38, highest: 70 71 74 76 77
## ---------------------------------------------------------------------------
## Sex 
##        n  missing distinct     Info      Sum     Mean      Gmd 
##      303        0        2    0.653      206   0.6799   0.4367 
## 
## ---------------------------------------------------------------------------
## ChestPain 
##        n  missing distinct 
##      303        0        4 
##                                                               
## Value      asymptomatic   nonanginal   nontypical      typical
## Frequency           144           86           50           23
## Proportion        0.475        0.284        0.165        0.076
## ---------------------------------------------------------------------------
## RestBP 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      303        0       50    0.995    131.7    19.41      108      110 
##      .25      .50      .75      .90      .95 
##      120      130      140      152      160 
## 
## lowest :  94 100 101 102 104, highest: 174 178 180 192 200
## ---------------------------------------------------------------------------
## Chol 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      303        0      152        1    246.7    55.91    175.1    188.8 
##      .25      .50      .75      .90      .95 
##    211.0    241.0    275.0    308.8    326.9 
## 
## lowest : 126 131 141 149 157, highest: 394 407 409 417 564
## ---------------------------------------------------------------------------
## Fbs 
##        n  missing distinct     Info      Sum     Mean      Gmd 
##      303        0        2    0.379       45   0.1485   0.2538 
## 
## ---------------------------------------------------------------------------
## RestECG 
##        n  missing distinct     Info     Mean      Gmd 
##      303        0        3     0.76   0.9901    1.003 
##                             
## Value          0     1     2
## Frequency    151     4   148
## Proportion 0.498 0.013 0.488
## ---------------------------------------------------------------------------
## MaxHR 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      303        0       91        1    149.6    25.73    108.1    116.0 
##      .25      .50      .75      .90      .95 
##    133.5    153.0    166.0    176.6    181.9 
## 
## lowest :  71  88  90  95  96, highest: 190 192 194 195 202
## ---------------------------------------------------------------------------
## ExAng 
##        n  missing distinct     Info      Sum     Mean      Gmd 
##      303        0        2     0.66       99   0.3267   0.4414 
## 
## ---------------------------------------------------------------------------
## Oldpeak 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      303        0       40    0.964     1.04    1.225      0.0      0.0 
##      .25      .50      .75      .90      .95 
##      0.0      0.8      1.6      2.8      3.4 
## 
## lowest : 0.0 0.1 0.2 0.3 0.4, highest: 4.0 4.2 4.4 5.6 6.2
## ---------------------------------------------------------------------------
## Slope 
##        n  missing distinct     Info     Mean      Gmd 
##      303        0        3    0.798    1.601   0.6291 
##                             
## Value          1     2     3
## Frequency    142   140    21
## Proportion 0.469 0.462 0.069
## ---------------------------------------------------------------------------
## Ca 
##        n  missing distinct     Info     Mean      Gmd 
##      299        4        4    0.783   0.6722   0.9249 
##                                   
## Value          0     1     2     3
## Frequency    176    65    38    20
## Proportion 0.589 0.217 0.127 0.067
## ---------------------------------------------------------------------------
## Thal 
##        n  missing distinct 
##      301        2        3 
##                                            
## Value           fixed     normal reversable
## Frequency          18        166        117
## Proportion      0.060      0.551      0.389
## ---------------------------------------------------------------------------
## AHD 
##        n  missing distinct 
##      303        0        2 
##                       
## Value         No   Yes
## Frequency    164   139
## Proportion 0.541 0.459
## ---------------------------------------------------------------------------</code></pre>
<p>Our target variable in this data is <code>AHD</code>. This flag identifies whether or not a patient has <a rel="nofollow" href="https://g.co/kgs/hT5ibs" target="_blank">Coronary Artery Disease</a>. If we can predict this accurately, clinicians could probably better treat these patients and hopefully help them avoid the symptoms of AHD like chest pain or worse, heart attacks.</p>
</div>
<div>
<h3>Data Wrangling</h3>
<p>For a predictive model I’ve opted to use a random forest model using the <code>ranger</code> implmentation which parallelizes the random forests algorithm in R. But first, some data cleaning is necessary. After replacing missing values, I’m going to split the data into test and training dataframes.</p>
<pre><code># Replace missing values
dat$Ca[is.na(dat$Ca)] &lt;- -1
dat$Thal[is.na(dat$Thal)] &lt;- &quot;missing&quot;

## 75% of the sample size
smp_size &lt;- floor(0.75 * nrow(dat))

## set the seed to make your partition reproducible
set.seed(123)
train_ind &lt;- sample(seq_len(nrow(dat)), size = smp_size)

train &lt;- dat[train_ind, ]
test &lt;- dat[-train_ind, ]

mod &lt;- ranger(AHD~., data=train, probability = TRUE, importance = &quot;permutation&quot;)

mod$prediction.error</code></pre>
<pre><code>## [1] 0.1326235</code></pre>
<p>Our quick and dirty check of the OOB prediction error tells us that our model appears to be doing okay at predicting <code>AHR</code>. Now the trick is to describe to our physicians and nurses why we believe someone is high risk for <code>AHR</code>. Before I learned of <code>lime</code>, I would have probably done something similar to the code below by first looking at which variables were most important in my trees.</p>
<pre><code>plot_importance &lt;- function(mod){
  tmp &lt;- mod$variable.importance
  dat &lt;- data.frame(variable=names(tmp),importance=tmp)
  ggplot(dat, aes(x=reorder(variable,importance), y=importance))+ 
    geom_bar(stat=&quot;identity&quot;, position=&quot;dodge&quot;)+ coord_flip()+
    ylab(&quot;Variable Importance&quot;)+
    xlab(&quot;&quot;)
}

# Plot the variable importance
plot_importance(mod)</code></pre>
<p><img src="https://i1.wp.com/nielsenmark.us/post/2018-11-09-exploring-models-with-lime_files/figure-html/plot-importance-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>After this, I probably would have taken a look at some partial dependence plots to get an idea of how those important variables are changing over the range of that variable. However, often the weakness of this approach is that I need to hold all other variables constant. And if I truly believe there are interactions between my variables, the partial dependence plot could change dramatically when other variables are changed.</p>
</div>
<div>
<h3>Explain the model with LIME</h3>
<p>Enter <code>lime</code>. As discussed above, the entire purpose of <code>lime</code> is to provide a local interpretable model to help us understand how our prediction would change if we tweak the other variables slightly in a lot of permutations. The first step to using <code>lime</code> in this specific case is to add some functions so that the <code>lime</code> package knows how to deal with the output of the <code>ranger</code> package. Once I have these I can use the combination of the <code>lime()</code> and <code>explain()</code> functions to get what I need. As in all multivariate linear models, we still have an issue… correlated explanatory varaibles. And depending on the number of variables in our original model, we may need to pair down our models to only look at the most “influential” or “important” variables. By default lime is going to use either forward-selection or pick the variables with the larges coefficients after correcting for multicollinearity using ridge regression or L2 penalization. As seen below, you can also select variables for the explanation using Lasso (aka L1 penalization) or use <code>xgboost</code> most important variables using the <code>&quot;tree&quot;</code> method.</p>
<pre><code># Train LIME Explainer
expln &lt;- lime(train, model = mod)


preds &lt;- predict(mod,train,type = &quot;response&quot;)
# Add ranger to LIME
predict_model.ranger &lt;- function(x, newdata, type, ...) {
  res &lt;- predict(x, data = newdata, ...)
  switch(
    type,
    raw = data.frame(Response = ifelse(res$predictions[,&quot;Yes&quot;] >= 0.5,&quot;Yes&quot;,&quot;No&quot;), stringsAsFactors = FALSE),
    prob = as.data.frame(res$predictions[,&quot;Yes&quot;], check.names = FALSE)
  )
}

model_type.ranger &lt;- function(x, ...) 'classification'


reasons.forward &lt;- explain(x=test[,names(test)!=&quot;AHD&quot;], explainer=expln, n_labels = 1, n_features = 4)
reasons.ridge &lt;- explain(x=test[,names(test)!=&quot;AHD&quot;], explainer=expln, n_labels = 1, n_features = 4, feature_select = &quot;highest_weights&quot;)
reasons.lasso &lt;- explain(x=test[,names(test)!=&quot;AHD&quot;], explainer=expln, n_labels = 1, n_features = 4, feature_select = &quot;lasso_path&quot;)
reasons.tree &lt;- explain(x=test[,names(test)!=&quot;AHD&quot;], explainer=expln, n_labels = 1, n_features = 4, feature_select = &quot;tree&quot;)</code></pre>
<p>Note: Using the current version of <code>lime</code> you may have issues with the <code>feature_select = &quot;lasso_path&quot;</code> option. To get the above code to run above you can install my tweaked version of <code>lime</code> <a rel="nofollow" href="https://github.com/nielsenmarkus11/lime" target="_blank">here</a>.</p>
</div>
<div>
<h3>Plotting explanations</h3>
<p>Now that we have all the explanations, one of my favorite features in the <code>lime</code> package is the <code>plot_explanations()</code> function. You can easily show the most important variables for each of our selection methods above and we can see that they are all very consistent in the choice of the top 4 most influential variables in predicting <code>AHD</code>.</p>
<pre><code>plot_explanations(reasons.forward)</code></pre>
<p><img src="https://i1.wp.com/nielsenmark.us/post/2018-11-09-exploring-models-with-lime_files/figure-html/unnamed-chunk-1-1.png?w=450&ssl=1" data-recalc-dims="1"></p>
<pre><code>plot_explanations(reasons.ridge)</code></pre>
<p><img src="https://i1.wp.com/nielsenmark.us/post/2018-11-09-exploring-models-with-lime_files/figure-html/unnamed-chunk-1-2.png?w=450&ssl=1" data-recalc-dims="1"></p>
<pre><code>plot_explanations(reasons.lasso)</code></pre>
<p><img src="https://i0.wp.com/nielsenmark.us/post/2018-11-09-exploring-models-with-lime_files/figure-html/unnamed-chunk-1-3.png?w=450&ssl=1" data-recalc-dims="1"></p>
<pre><code>plot_explanations(reasons.tree)</code></pre>
<p><img src="https://i0.wp.com/nielsenmark.us/post/2018-11-09-exploring-models-with-lime_files/figure-html/unnamed-chunk-1-4.png?w=450&ssl=1" data-recalc-dims="1"></p>
<p>Thanks for reading this quick tutorial on <code>lime</code>. There is much more of this package that I want to explore. Particulary its use for image and text classifications. Then the only real question left is… How do I get one of those cool hex stickers for <code>lime</code>? <img alt="😉" src="https://s.w.org/images/core/emoji/2.4/72x72/1f609.png" class="wp-smiley"></p>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://nielsenmark.us/2018/11/09/exploring-models-with-lime/"> Programming on nielsenmark.us</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/-geXQH1VO6k" height="1"></a></div>