<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/auc-meets-the-wilcoxon-mann-whitney-u-statistic/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://blog.revolutionanalytics.com/2017/03/auc-meets-u-stat.html"> Revolutions</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p><em>by Bob Horton, Senior Data Scientist, Microsoft</em></p>
<p>The area under an ROC curve (AUC) is commonly used in machine learning to summarize the performance of a predictive model with a single value. But you might be surprised to learn that the AUC is directly connected to the Mann-Whitney U-Statistic, which is commonly used in a robust, non-parametric alternative to Student’s t-test. Here I’ll use ‘literate analysis’ to demonstrate this connection and illustrate how the two measures are related.</p>
<p>In previous posts on <a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2016/08/roc-curves-in-two-lines-of-code.html">ROC curves</a> and <a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2016/11/calculating-auc.html">AUC</a> I described some simple ways to visualize and calculate these objects. Here is the simple data we used earlier to illustrate AUC:</p>
<pre><code>category &lt;- c(1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0)
prediction &lt;- rev(seq_along(category))
prediction[9:10] &lt;- mean(prediction[9:10])

library('pROC')
(official_auc &lt;- auc(roc(category, prediction)))</code></pre>
<pre><code>## Area under the curve: 0.825</code></pre>
<p>Here <code>category</code> is a vector of Boolean labels marking the true status of a sequence of cases. The other vector, <code>prediction</code>, represents a set of numeric scores as would normally be generated by some type of measurement or classifier algorithm. These scores could represent, for example, the expected probability of an object being a cat. But they don’t need to be probabilities; any value indicating the relative strength of the classifier’s confidence that the object is a cat can work, as long as the scores let us sort the cases into some order. Our fake scores are designed to put the cases in the order they start with, except that the scores of two cases have been replaced with their average; this gives us some instances where the scores are tied, which is a fairly reasonable condition we should be sure to handle. For this dataset the ‘official’ value for AUC is 0.825; when we try various other ways to calculate AUC, this is the number we want to see.</p>
<div>
<h3>Computing AUC from the U Statistic</h3>
<p>From <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test">Wikipedia</a> we learn that</p>
<p><span>\[ {AUC}_1 = {U_1 \over n_1n_2} \]</span></p>
<p>where <span>\(U_1\)</span> is the Mann-Whitney U statistic, also known as the Wilcoxon rank-sum test statistic, or some combination and/or permutation of those names. That seems like a strange claim, but it is easy enough to test:</p>
<pre><code>auc_wmw &lt;- function(labels, scores){
  labels &lt;- as.logical(labels)
  pos &lt;- scores[labels]
  neg &lt;- scores[!labels]
  U &lt;- as.numeric(wilcox.test(pos, neg)$statistic)
  U/(length(pos) * length(neg))
}

auc_wmw(category, prediction)</code></pre>
<pre><code>## Warning in wilcox.test.default(pos, neg): cannot compute exact p-value with
## ties</code></pre>
<pre><code>## [1] 0.825</code></pre>
<p>The <code>wilcox.test</code> function warns us that we “cannot compute exact p-value with ties”, but for the current exercise let’s just savor the fact that it got the AUC exactly right.</p>
</div>
<p>To start to decipher the U statistic, let’s calculate it ourselves instead of using the <code>wilcox.test</code> function as a black box. Going back to Wikipedia we learn that</p>
<p><span>\[ U_1 = R_1 – {n_1(n_1+1) \over 2} \]</span></p>
<p>where <span>\(R_1\)</span> is the sum of the ranks of the positive cases, and <span>\(n_1\)</span> is the number of positive cases. You can calculate a related score (<span>\(U_2\)</span>) with the negative cases (of which there are <span>\(n_2\)</span>); these two scores are complementary in that they always add up to <span>\(n_1 n_2\)</span>, similar to the way that flipping positives and negatives when calculating an ROC curve gives you an AUC that is one minus the AUC of the unflipped curve. Anyway, now we have:</p>
<pre><code>auc_wmw2 &lt;- function(labels, scores){
  labels &lt;- as.logical(labels)
  n1 &lt;- sum(labels)
  n2 &lt;- sum(!labels)
  R1 &lt;- sum(rank(scores)[labels])
  U1 &lt;- R1 - n1 * (n1 + 1)/2
  U1/(n1 * n2)
}

auc_wmw2(category, prediction)</code></pre>
<pre><code>## [1] 0.825</code></pre>
<p>Base R’s <code>rank</code> function assigns the lowest rank value (1, if there are no ties) to the lowest score, and by default it averages the ranks for tied scores, which is exactly what we need. Now I’ll try to convince you that it makes sense that <code>U1/(n1 * n2)</code> is equal to AUC.</p>
<div>
<h3>Graphical interpretation of the U Statistic</h3>
<p>First we plot the ranks of all the scores as a stack of horizontal bars, and color them by the labels.</p>
<pre><code># simplify the syntax for drawing rectangles
rectangle &lt;- function(x, y, width, height, density=12, angle=-45, ...) 
  polygon(c(x,x,x+width,x+width), c(y,y+height,y+height,y), 
          density=density, angle=angle, ...)

U_illustration_part1 &lt;- function(labels, scores){
  # put cases in order by score
  sort_order &lt;- order(scores)
  labels &lt;- labels[sort_order]
  scores &lt;- scores[sort_order]
  
  # count the cases
  n &lt;- length(labels)
  
  # find overall rank for each case by score
  ranks &lt;- rank(scores)
  
  # start with an empty plot
  plot(c(0, n), c(0, n), type='n',
       xlab=&quot;rank&quot;, ylab=&quot;case&quot;, asp=1)
  
  # draw a grid in the background
  abline(h=0:n, col=&quot;lightblue&quot;)
  abline(v=0:n, col=&quot;lightblue&quot;)
  
  # plot bars representing ranks of all cases
  mapply(rectangle, x=0, y=(n - 1):0,  # starting from the top 
         width=ranks, height=1, 
         density=NA, border=&quot;black&quot;, lwd=2, col=c(&quot;red&quot;, &quot;green&quot;)[1 + labels])
  
  legend(&quot;topright&quot;, legend=c(&quot;negative case&quot;, &quot;positive case&quot;), 
         text.col=c(&quot;red&quot;, &quot;green&quot;), bty='o', box.lwd=1, inset=0.1)
}

U_illustration_part1(labels=category, scores=prediction)</code></pre>
<p><a target="_blank" rel="nofollow" href="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01b8d26ab202970c-pi"><img title="U_illustration_part1-1" alt="U_illustration_part1-1" src="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01b8d26ab202970c-800wi" border="0"></a></p>
<p>Now consider only the green bars, representing the ranks of the positive cases. We’ll stack them on top of one another, and slide them horizontally as needed to get a nice even stairstep on the right edge:</p>
<pre><code>U_illustration_part2 &lt;- function(labels, scores){
  # sort the cases
  sort_order &lt;- order(scores)
  labels &lt;- labels[sort_order]
  scores &lt;- scores[sort_order]

  # count positive and negative cases
  n1 &lt;- sum(labels)  # number of positive cases
  n2 &lt;- sum(!labels) # number of negative cases
  
  # find the overall ranks for the positive cases
  ranks &lt;- rank(scores)
  pos_ranks &lt;- ranks[as.logical(labels)]
  
  # how far to slide each bar to make stairsteps on the right hand edge
  x_offset &lt;- n2 + (1:n1) - pos_ranks

  # start with an empty plot  
  plot(c(0, n2 + n1), c(0, n1), type='n', asp=1,
       xlab=&quot;n2 + n1 divisions&quot;, ylab=&quot;n1 divisions&quot;)
  
  # plot bars for ranks of positive cases
  mapply(rectangle, x=x_offset, y=(n1 - 1):0, 
         width=pos_ranks, height=1,
         density=NA, border=&quot;darkgreen&quot;, lwd=2, col=&quot;green&quot;)
  
  # draw the grid
  abline(h=0:n1, col=&quot;lightblue&quot;)
  abline(v=0:(n1 + n2), col=&quot;lightblue&quot;)
  
  # mark the area we remove, and the area we keep
  rectangle(n2, 0, n1, n1, density=10, col=&quot;red&quot;, lwd=1)
  rectangle(0, 0, n2, n1, density=0, col=&quot;black&quot;, lty=2, lwd=3)
  
  # draw a scaled version of the &quot;official&quot; ROC curve on top
  roc_obj &lt;- roc(labels, scores)
  roc_df &lt;- with(roc_obj, data.frame(FPR=rev(1 - specificities), <br>    TPR=rev(sensitivities)))
  with(roc_df, lines(n2*FPR, n1*TPR, type='l', lwd=4, col=&quot;blue&quot;))
}

U_illustration_part2(labels=category, scores=prediction)</code></pre>
<p><a target="_blank" rel="nofollow" href="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb0983835b970d-pi"><img title="U_illustration_part2-1" alt="U_illustration_part2-1" src="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb0983835b970d-800wi" border="0"></a></p>
<p>The total area of all the green bars equals the sum of the ranks of the positive cases (<span>\(R_1\)</span> in the equation). The red hatched box on the right hand side represents the part we’ll subtract. Note that the total area of green inside the red hatched box is <span>\({n_1(n_1+1) \over 2}\)</span>, which you may recognize as the sum of the integers <span>\(1\)</span> through <span>\(n_1\)</span>. The dashed rectangle on the left side shows the part we keep; the green area inside this dashed rectangle represents the value of <span>\(U_1\)</span>.</p>
<p>Here <code>n1</code>, the number of steps on the y axis, is obviously just the number of green bars. It may take a bit more contemplation to see that <code>n2</code>, the number of steps on the x axis inside the dashed rectangle, reflects the number of red bars we removed (the negative cases). This plot is basically a transliteration of the formula for <span>\(U_1\)</span> into a figure.</p>
<p>The whole grid of the dashed rectangle thus has an area of <code>n1 * n2</code> steps squared, and <code>U1</code> is the area under the curve (in “square steps”). <code>U1/(n1 * n2)</code> is the area under the curve as a fraction of the total area of the rectangle. In a traditional ROC curve the number of steps along each axis is normalized to 1, so that AUC is a fraction of a 1 by 1 area; once we normalize <span>\(U_1\)</span> to the area of the <span>\(n_1\)</span> by <span>\(n_2\)</span> rectangle, it equals AUC.</p>
<p>The blue line is the ‘official’ ROC curve scaled by <code>n2</code> on the x (FPR) axis and by <code>n1</code> on the y (TPR) axis, and it matches the left edge of our stacked bars except where scores are tied. The half-unit area (where two scores are tied) is split vertically rather than diagonally, but it is still a half unit, so you can see that the area of stacked green bars within the dashed rectangle corresponds exactly to the area under the ROC curve.</p>
</div>
<div>
<h3>Visualizing the U Statistic on a larger dataset</h3>
<p>Let’s generate similar figures using the example from the <a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2016/08/roc-curves-in-two-lines-of-code.html">earlier post</a>, where the <code>test_set</code> dataframe has a Boolean column called <code>bad_widget</code> that labels cases, and a vector called <code>glm_response_scores</code> containing scores from a logistic regression. We’ll use the same plotting code as above.</p>
<pre><code>U_illustration_part1(test_set$bad_widget, glm_response_scores)</code></pre>
<p><a target="_blank" rel="nofollow" href="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb0983838f970d-pi"><img title="Big_example_U_fig1-1" alt="Big_example_U_fig1-1" src="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01bb0983838f970d-800wi" border="0"></a></p>
<p>Now we remove the red bars, stack the green ones together, and slide them to make even stair-steps on the right side. The left edge of these stacked bars forms an ROC curve.</p>
<pre><code>U_illustration_part2(test_set$bad_widget, glm_response_scores)</code></pre>
<p><a target="_blank" rel="nofollow" href="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01b8d26ab258970c-pi"><img title="Big_example_U_fig2-1" alt="Big_example_U_fig2-1" src="http://revolution-computing.typepad.com/.a/6a010534b1db25970b01b8d26ab258970c-800wi" border="0"></a></p>
<p>Again, the blue line is the ‘official’ ROC curve calculated with the <code>pROC</code> package, then scaled by <code>n1</code> and <code>n2</code> on the y and x axes, respectively. In this dataset there are no tied scores so the ROC curve has no diagonal segments, and it matches the jagged contour of the left edge of our stacked rank bars exactly.</p>
</div>
<div>
<h3>Check results with tied scores</h3>
<p>Although this graphical approach does not draw the diagonal segments of an ROC curve correctly, calculation based on the U-statistic does get the AUC right even in the presence of ties. To demonstrate this more dramatically, we’ll round off the original scores to one decimal place so that many of the scores round to the same value, producing groups of tied scores (see the ROC plots for the rounded values in the earlier <a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2016/11/calculating-auc.html">AUC</a> post).</p>
<pre><code>rounded_scores &lt;- round(glm_response_scores, digits=1)</code></pre>
<p>Now we can use our U statistic – based functions to compute the AUC on both the original and rounded versions:</p>
<pre><code>data.frame(
  auc_pROC = auc(roc(test_set$bad_widget, glm_response_scores)),
  auc_wmw = auc_wmw(test_set$bad_widget == &quot;TRUE&quot;, glm_response_scores),
  auc_wmw2 = auc_wmw2(test_set$bad_widget == &quot;TRUE&quot;, glm_response_scores),
  auc_pROC_ties = auc(roc(test_set$bad_widget, rounded_scores)),
  auc_wmw_ties = auc_wmw(test_set$bad_widget, rounded_scores),
  auc_wmw2_ties = auc_wmw2(test_set$bad_widget, rounded_scores)
)</code></pre>
<pre><code>##   auc_pROC  auc_wmw auc_wmw2 auc_pROC_ties auc_wmw_ties auc_wmw2_ties
## 1 0.903698 0.903698 0.903698     0.8972779    0.8972779     0.8972779</code></pre>
<p>Note that the “official” AUC is slightly smaller when the scores are rounded, and that both versions of our Wilcoxon-Mann-Whitney U-Statistic based functions get the AUC exactly right, with or without ties.</p>
<p>This post has focused on the value and interpretation of the <strong>U statistic</strong> itself, but we have not discussed the <strong>U test</strong>, which uses this statistic (along with sample sizes) to evaluate statistical significance. According to <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Mann%E2%80%93Whitney_U_test">Wikipedia</a>, the U test evauates the null hypothesis that it is “equally likely that a randomly selected value from one sample will be less than or greater than a randomly selected value from a second sample.” As <a target="_blank" rel="nofollow" href="http://blog.revolutionanalytics.com/2016/11/calculating-auc.html">illustrated earlier</a>, AUC reflects a similar probability, that a randomly chosen positive case will receive a higher score from your model than a randomly chosen negative case. Although exploring how this test calculates p-values is beyond the scope of this post, it should not be surprising that metrics embodying these sorts of probabilities can form the basis of such tests.</p>
</div>
</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="http://blog.revolutionanalytics.com/2017/03/auc-meets-u-stat.html"> Revolutions</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 title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/pDOi-86OEiM"></a></div>