<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/why-no-exact-permutation-tests-at-scale/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.win-vector.com/blog/2018/02/why-no-exact-permutation-tests-at-scale/"> R – Win-Vector Blog</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Here at <a rel="nofollow" href="http://www.win-vector.com/" target="_blank">Win-Vector LLC</a> we like permutation tests. Our team has written on them (for example: <a rel="nofollow" href="http://www.win-vector.com/blog/2015/08/how-do-you-know-if-your-data-has-signal/" target="_blank">How Do You Know if Your Data Has Signal?</a>) and they are used to estimate significances in our <a rel="nofollow" href="https://cran.r-project.org/package=sigr" target="_blank">sigr</a> and <a rel="nofollow" href="https://cran.r-project.org/package=WVPlots" target="_blank">WVPlots</a> <a rel="nofollow" href="https://cran.r-project.org/" target="_blank">R</a> packages. For example permutation methods are used to estimate the significance reported in the following <a rel="nofollow" href="https://winvector.github.io/WVPlots/reference/ROCPlot.html" target="_blank">ROC plot</a>.</p>
<p><a rel="nofollow" href="https://github.com/WinVector/Examples/blob/master/WVPlots/WVPlots_concept.md" target="_blank"><img border="0" src="https://i2.wp.com/www.win-vector.com/blog/wp-content/uploads/2017/04/NewImage-12.png?w=660" alt="NewImage" title="NewImage.png" data-recalc-dims="1"></a></p>
<p>Permutation tests have their own literature and issues (examples: <a rel="nofollow" href="http://www.springer.com/us/book/9780387202792" target="_blank"><em>Permutation, Parametric and Bootstrap Tests of Hypotheses,</em> Springer-Verlag, NY, 1994 (3rd edition, 2005)</a>, <a rel="nofollow" href="http://faculty.wcas.northwestern.edu/~iac879/wp/RDDPermutations.pdf" target="_blank">2</a>, <a rel="nofollow" href="http://thomasleeper.com/Rcourse/Tutorials/permutationtests.html" target="_blank">3</a>, and <a rel="nofollow" href="http://polmeth.rice.edu/files/2016/07/CaugheyDafoeMiratrix-uaw652.pdf" target="_blank">4</a>).</p>
<p>In our <a rel="nofollow" href="https://www.r-project.org/" target="_blank">R</a> packages the permutation tests are estimated by a sampling procedure, and not computed exactly (or deterministically). It turns out this is likely a necessary concession; a complete <em>exact</em> permutation test procedure at scale would be big news. Please read on for my comments on this issue.</p>
<p><span></span></p>
<p>
<h2>Introduction</h2>
<p>In a basic permutation test we check if a relation is important by deliberately breaking the relation and seeing if outcome quality changes. The idea is: we build a deliberately damaged data set by permuting the values in one or more columns independent of other columns (breaking any per-instance relations between the permuted and non-permuted columns). If damage is unnoticeable, then what was damaged may not be critical. If damage is large, then what was damaged may be critical.</p>
<p>Permutation tests are related to other empirical re-sampling methods (such as the <a rel="nofollow" href="https://projecteuclid.org/euclid.aos/1176344552" target="_blank">bootstrap</a>) but have their own flavor. One of the most famous uses is <a rel="nofollow" href="https://www.statistik.uni-dortmund.de/useR-2008/slides/Strobl+Zeileis.pdf" target="_blank">random forest variable importance testing</a>.</p>
<h2>Definitions</h2>
<p>Let’s define our terms using a simple example:</p>
<p>Suppose we have a sequence of examples <code>(x_i, v_i, y_i)</code> where <code>i</code> takes integer values from 1 to n and is our “instance number”. <code>x_i</code> is one set of variable values (variables we are using, but not estimating the importance of) associated with the <code>i</code>-th example and <code>v_i</code> is another set of observations associated (a variable or group or variables we are using and estimating the importance of), and <code>y_i</code> is the known outcome we are trying to match.</p>
<p>Define:</p>
<ul>
<li><strong><code>predict(x_i,v_i)</code></strong>: a function or model (probably trained on earlier data) that takes two arguments and returns a prediction (a prediction being a number in the case of regression, or a probability or class label in the case of classification). For n-vectors <code>x</code> and <code>v</code> we will say <code>predict(x,v)</code> is the n-vector whose <code>i</code>-th entry is <code>predict(x_i,v_i)</code>.</li>
<li><strong><code>loss(y, yhat)</code></strong>: a function that takes the n-vector <code>y</code> of known outcomes and an n-vector <code>yhat</code> of matching predictions and returns a criticism or score. Often <code>loss()</code> is a simple additive function (such as <a rel="nofollow" href="https://en.wikipedia.org/wiki/Residual_sum_of_squares" target="_blank">SSE</a> = <code>sum_{i:1 ...n} (y_i - yhat_i)^2</code>), or it can be more involved such as <a rel="nofollow" href="http://www.win-vector.com/blog/2016/07/a-budget-of-classifier-evaluation-measures/" target="_blank">1-AUC</a>.</li>
<li>Let <code>Sn</code> be the set of all <a rel="nofollow" href="https://en.wikipedia.org/wiki/Permutation" target="_blank">permutations</a> of integers 1 through n. For a permutation <code>p</code> in <code>Sn</code> and an n-vector <code>v</code> we write <code>p(v)</code> to represent the n-vector whose <code>i</code>-th entry is <code>v_{p(i)}</code>.</li>
</ul>
<p>For each of these functions we will say</p>
<p>In this notation a permutation test of the effect of the <code>v</code> variables is constructed as follows:</p>
<ul>
<li>Compute: lV = <code>loss(y, predict(x,v))</code>, the observed loss value.</li>
<li>For an arbitrary real number T, compute or estimate: pV(T) = <code>(1/|Sn|) sum_{p in Sn} (loss(y, predict(x, p(v))) ≥ T)</code> (where <code>loss(y, predict(x, p(v))) ≥ T</code> is valued at 1 if true and 0 if false).</li>
</ul>
<p>The idea is: lV is the loss we saw, and pV(lV) is how often a random relation achieves a loss at least as large. If 1-pV(lV) is small we can say the difference is significant in a frequentist sense.</p>
<p>The problem is that <code>|Sn| = n!</code> which is typically enormous. So we don’t want to actually compute the above sum in practice. We can <em>estimate</em> pV(T) by sampling terms uniformly at random (which works especially well as the terms all have bounded size), but this has the (minor) undesirable side-effect that our estimator isn’t deterministic (two different honest experimenters will typically get different results for the same data).</p>
<h2>Why no exact solution at scale?</h2>
<p>The question is: can we directly and exactly calculate pV(T) (instead of estimating it)?</p>
<p>The answer is: in the general case to do so would be a very big new result. There is definitely not a publicly known trusted completely general fast method to do so.</p>
<p>Let’s go over why that is.</p>
<p>Take the simple case of “additive loss”, where: <code>loss(y, yhat)</code> is defined as <code>sum_{i = 1...n} loss(y_i, yhat_i)</code>. Losses such as sum of square errors, and number of mis-classifications are of this simple form. For such a loss we can define a matrix <code>l</code> such that:</p>
<ul>
<li><code>l(i,j) = loss(y_i, predict(x_i, v_j))</code> where <code>i,j</code> are both integers chosen from 1 to n.</li>
</ul>
<p>We are going to further restrict to the case where <code>loss(y_i, predict(x_i, v_j))</code> is always 0 or 1 (such as in counting mis-classifications). The issue is: it is easy to confirm that every n by n matrix of 0/1 entries can be generated by picking the appropriate <code>predict()</code> example. For example: encode a given 0/1 matrix <code>M</code> as: <code>y = 0</code>, <code>x = v = (1, 2, ..., n)</code>, <code>predict(x_i, v_j) = M(x_i, x_j)</code>, and <code>loss(y, p)</code> = 0 if <code>y==p</code> and <code>1</code> otherwise.</p>
<p>If we could perfectly calculate arbitrary pV(T) we would be able to compute <code>P[loss(y, predict(x, p(v))) ≥ n] = pV(n)</code> where <code>P[]</code> is the probability computed over the uniform distribution of <code>p in Sn</code>. If we could always exactly compute that efficiently we could already solve a deep open problem: computing the 0/1 <a rel="nofollow" href="https://en.wikipedia.org/wiki/Permanent" target="_blank">permanent</a> quickly. So such a method would be big news.</p>
<p>Here is our derivation:</p>
<pre>
 P[loss(y, predict(x, p(v))) ≥ n] 
  = (1/|Sn|) sum_{p in Sn} (loss(y, predict(x, p(v))) ≥ n)
  = (1/|Sn|) sum_{p in Sn} ((sum_{i = 1...n} l(i,p(i))) ≥ n)
  = (1/|Sn|) sum_{p in Sn} product_{i = 1...n} l(i,p(i))
  = (1/|Sn|) permanent(l)
  = (1/(n!)) permanent(l)
</pre>
<p>We are using the fact that while sum doesn’t equal product in general, for a 0/1 n-term loss we have <code>loss(y, predict(x, p(v)) ≥ n</code> if and only if <code>product_{i = 1...n} l(i,p(i)) == 1</code>. Now <code>n!</code> is easy to compute; so if <code>P[loss(y, predict(x, p(v)) ≥ n]</code> were easy to compute, then <code>permanent(l)</code> would also be easy to compute.</p>
<p>And we now have the catch: <a rel="nofollow" href="https://en.wikipedia.org/wiki/Sharp-P-completeness_of_01-permanent" target="_blank">it would be big news if somebody could efficiently exactly compute arbitrary 0/1-permanents</a>. As we said: our restricted loss function is already rich enough to every 0/1 matrix.</p>
<p>To be able to exactly compute pV(n) for <em>all</em> examples (our desired goal) would mean we are also able to compute the permanent for <em>all</em> 0/1 matrices (something thought to be hard). It is the expressive power that hangs us: being able to compute the loss for only <em>some</em> examples would not be an issue, as that would not imply the power exactly compute <em>all</em> 0/1 permanents.</p>
	<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://www.win-vector.com/blog/2018/02/why-no-exact-permutation-tests-at-scale/"> R – Win-Vector Blog</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/ZhZdnjKKGNs" height="1"></a></div>