<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/beware-the-argument-the-flint-water-crisis-and-quantiles/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://staff.math.su.se/hoehle/blog/2017/06/18/quantiles.html"> Theory meets practice...</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<h2>Abstract</h2>
<p>If your tap water suddenly becomes brown while authorities claim everything is okay, you start to worry. <span data-cites="langkjaerbain2017">Langkjær-Bain (2017)</span> tells the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Flint_water_crisis" target="_blank">Flint Water Crisis</a> story from a statistical viewpoint: essentially the interest is in whether the 90th percentile in a sample of lead concentration measurements is above a certain threshold or not. We illustrate how to perform the necessary calculations with R’s quantile function and show that the type-argument of the function matters.</p>
<p><center><br>
<img src="https://i1.wp.com/staff.math.su.se/hoehle/blog/figure/source/2017-06-18-quantiles/SAWTOOTH-1.png?w=456" data-recalc-dims="1"><br>
</center></p>
<p> <a rel="nofollow" href="http://creativecommons.org/licenses/by-sa/4.0/" target="_blank"><img alt="Creative Commons License" src="https://i1.wp.com/i.creativecommons.org/l/by-sa/4.0/88x31.png?w=456&ssl=1" data-recalc-dims="1"></a> This work is licensed under a <a rel="nofollow" href="http://creativecommons.org/licenses/by-sa/4.0/" target="_blank">Creative Commons Attribution-ShareAlike 4.0 International License</a>. The markdown+Rknitr source code of this blog is available under a <a rel="nofollow" href="https://www.gnu.org/licenses/gpl-3.0.html" target="_blank">GNU General Public License (GPL v3)</a> license from <a rel="nofollow" href="https://github.com/hoehleatsu/hoehleatsu.github.io" target="_blank"><img alt="github" src="https://i2.wp.com/staff.math.su.se/hoehle/blog/images/GitHub-Mark-32px.png?w=456" data-recalc-dims="1"></a>.</p>
<h2>Introduction</h2>
<p>In a recent <a rel="nofollow" href="https://www.significancemagazine.com/" target="_blank">Significance</a> article, <span data-cites="langkjaerbain2017">Langkjær-Bain (2017)</span> tells the story about the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Flint_water_crisis" target="_blank">Flint water crisis</a>. In 2014 the city of Flint, Michigan, USA, decided to change its water supply to Flint River. Due to insufficient treatment of the water with corrosion inhibitors the lead concentration in the drinking water increased, because lead in the aging pipes leached into the water supply. This created serious health problems – as explained in this <a rel="nofollow" href="https://youtu.be/QnDQFivtCd0" target="_blank">short video summary</a>. In this blog post we investigate further the computation of the 90th percentile of the tap water lead concentration samples described in <span data-cites="langkjaerbain2017">Langkjær-Bain (2017)</span>. Quantile estimation in this context has already been discussed in a recent blog entry entitled <a rel="nofollow" href="https://blogs.sas.com/content/iml/2017/05/17/quantiles-flint-water-crisis.html" target="_blank">Quantiles and the Flint water crisis</a> by <a rel="nofollow" href="https://blogs.sas.com/content/author/rickwicklin/" target="_blank">Rick Wicklin</a>.</p>
<p>The monitoring of drinking water in the US is regulated by the <a rel="nofollow" href="https://www.epa.gov/dwreginfo/lead-and-copper-rule" target="_blank">Lead and Copper Rule</a> of the United States Environmental Protection Agency. The entire text of the rule is available as electronic code of federal regulation (e-CFR) Title 40: Protection of Environment, <a rel="nofollow" href="https://www.ecfr.gov/cgi-bin/text-idx?SID=531617f923c3de2cbf5d12ae4663f56d&mc=true&node=sp40.23.141.i&rgn=div6" target="_blank">Part 141 – National Primary Drinking Water Regulations</a>. In particular the regulation defines a sampling plan for collecting tap water samples. The size of the sample depends on the number of people the water system serves. In case this number is bigger than 100,000 a sample of 100 sites is needed. If there are 10,001-100,000 people served, then a sample from 60 sites is needed. For systems serving fewer than 10,000 sizes of 40, 20, 10 and 5 are defined – see <a rel="nofollow" href="https://www.ecfr.gov/cgi-bin/text-idx?SID=531617f923c3de2cbf5d12ae4663f56d&mc=true&node=sp40.23.141.i&rgn=div6#se40.25.141_186" target="_blank">§141.86(c)</a> of the rule for details. Of interest for this blog post is that action needs to be taken, if too many of the samples are above a given threshold of 15 part per billion (ppb):</p>
<p><center><br>
<img src="https://i1.wp.com/staff.math.su.se/hoehle/blog/figure/source/2017-06-18-quantiles/ecfr141.png?w=456" data-recalc-dims="1"><br>
</center><br>
<font color="bbbbbb">Figure: Excerpt of <a rel="nofollow" href="https://www.ecfr.gov/cgi-bin/text-idx?SID=531617f923c3de2cbf5d12ae4663f56d&mc=true&node=sp40.23.141.i&rgn=div6#se40.25.141_180" target="_blank">CFR 40 §141.80</a> with the quantile calculation highlighted in yellow.</font></p>
<p><p>We note the explicit duality in the CFR between the <em>in more than 10%</em> and the 90% quantile in the text. However, we also note that it is not clear how to proceed, if the number calculated in (c)(3)(ii) is <strong>not an integer</strong>. This is not a problem per se, because the CFR itself only operates with samples sizes 100, 60, 40, 20, 10 and 5, where 0.9 times the sample size always gives an integer number. But if one for some reason does not obtain exactly the required number this quickly can become an issue as we shall see below.</p>
<h2>The Flint Data</h2>
<p>The data of the spring 2015 Flint water supply monitoring conducted by the Michigan Department of Environmental Quality are presented in the figure on p. 20 of <span data-cites="langkjaerbain2017">Langkjær-Bain (2017)</span>. Alternatively, they can also be taken directly from the Blog entry <a rel="nofollow" href="https://blogs.sas.com/content/iml/2017/05/17/quantiles-flint-water-crisis.html" target="_blank">Quantiles and the Flint water crisis</a> by <a rel="nofollow" href="https://blogs.sas.com/content/author/rickwicklin/" target="_blank">Rick Wicklin</a>.</p>
<div>
<pre><code>##Read flint water monitoring data (the column lead is the measurement)
flint &lt;-<span> </span><span>read.csv</span>(<span>file=</span><span>file.path</span>(filePath,<span>&quot;flint.csv&quot;</span>))
##Sort the measured lead values in ascending order
lead  &lt;-<span> </span><span>sort</span>(flint$lead)
##Number of observations
n &lt;-<span> </span><span>length</span>(lead)</code></pre>
</div>
<p><img src="https://i1.wp.com/staff.math.su.se/hoehle/blog/figure/source/2017-06-18-quantiles/HISTOGRAM-1.png?w=456" data-recalc-dims="1"></p>
<p>The proportion of observations above the 15 ppb threshold is</p>
<div>
<pre><code><span>mean</span>(lead ><span> </span><span>15</span>)</code></pre>
</div>
<pre><code>## [1] 0.1126761</code></pre>
<p>In other words, the proportion 11.3% is above the 10% action threshold and, hence, something needs to be done. However, as <span data-cites="langkjaerbain2017">Langkjær-Bain (2017)</span> describes, the story is a little more complicated, because two of the values above 15 were removed with the argument that they originated from sites which are not covered by the sampling plan. Only private households at high risk, i.e. with lead pipelines, are supposed to be sampled. As one can read in the article the removal is highly controversial, in particular, because the proportion of critical observations falls below the 10% action threshold when these two values are removed. For this blog entry, we will, however, work with the full <span>\(n=71\)</span> sample and focus on the <strong>quantile aspect</strong> of the calculation.</p>
<h2>Percentages and Quantiles</h2>
<p>Let <span>\(n\)</span> denote the size of the selected sample, which in our case is <span>\(n=71\)</span>. If more than 10% of the sample is to be above 15 ppb, this means that <span>\(\lfloor 0.1\cdot n + 1\rfloor\)</span> of the samples need to be above 15 ppb, where <span>\(\lfloor y \rfloor\)</span> denotes the largest integer less than or equal to <span>\(y\)</span>. We shall denote this the <strong>number of critical samples</strong>. If we divide this number by <span>\(n\)</span> we get the actual proportion of critical samples needed before action. It is worthwhile noting the difference between this critical proportion and the 10% threshold illustrated by the sawtooth curve in the figure below. The explanation for these sawtooth step-spikes is the discreteness of the decision problem (i.e. <span>\(x\)</span> out of <span>\(n\)</span>).</p>
<p><img src="https://i1.wp.com/staff.math.su.se/hoehle/blog/figure/source/2017-06-18-quantiles/SAWTOOTH-1.png?w=456" data-recalc-dims="1"></p>
<p>Turning to the equivalent <em>90% quantile is above 15 ppm</em> decision criterion of the CFR, one will need to determine the 90% quantile from the finite sample of lead concentration measurements. How to estimate quantiles with statistical software is discussed in the excellent survey of <span data-cites="hyndman_fan1996">Hyndman and Fan (1996)</span>. In their work they describe the nine different quantile estimators implemented in R’s <code>quantile</code> function. The estimators are all based on the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Order_statistic" target="_blank">order statistic</a> of the sample, i.e. let <span>\(x_{(1)} \leq x_{(2)} &lt; \cdots \leq x_{(n)}\)</span> denote the ordered lead concentration measurements. Then the simplest estimator for the <span>\(p\cdot 100\%\)</span> quantile is</p>
<p><span>\[<br>
\hat{x}_p = \min_{k} \left\{\hat{F}(x_{(k)}) \geq p\right\} =<br>
x_{(\lceil n \cdot p\rceil)},<br>
\]</span></p>
<p>where <span>\(\hat{F}\)</span> is the <a rel="nofollow" href="https://en.wikipedia.org/wiki/Empirical_distribution_function" target="_blank">empirical cumulative distribution</a> function of the sample and <span>\(\lceil y \rceil\)</span> denotes the smallest integer greater than or equal to <span>\(y\)</span>. This method corresponds to R’s <code>type=1</code>. For our application <span>\(p\)</span> would be 0.9 and the 90th percentil for the Flint data is thus</p>
<div>
<pre><code><span>c</span>(<span>&quot;manual.90%&quot;</span>=lead[<span>ceiling</span>(n*p)], <span>R=</span><span>quantile</span>(lead, <span>probs=</span>p, <span>type=</span><span>1</span>))</code></pre>
</div>
<pre><code>## manual.90%      R.90% 
##         18         18</code></pre>
<p>which is above the action threshold of 15 ppb. It is also important to note that when <span>\(x_{(\lceil n \cdot 0.9\rceil)} > 15\> \text{ppb}\)</span> then a total <span>\(n – \lceil n \cdot 0.9\rceil + 1\)</span> samples are above 15 ppm. In other words the proportion of samples above 15 ppm in this situation is <span>\[ \frac{n – \lceil n \cdot 0.9\rceil + 1}{n}.\]</span> To show the duality between the <em>more than 10% critical samples</em> and the 90% quantile being above 15 ppm we thus need to show that <span>\((n – \lceil n \cdot 0.9\rceil + 1)/n = (\lfloor 0.1\cdot n + 1\rfloor)/n\)</span>. This is possible since the following <a rel="nofollow" href="https://en.wikipedia.org/wiki/Floor_and_ceiling_functions" target="_blank">relations hold for the floor and ceiling functions</a>: <span>\[ \begin{align*}<br>
– \lceil y \rceil &amp;= \lfloor -y \rfloor \quad \text{and} \\<br>
n + \lfloor y \rfloor &amp;= \lfloor n + y \rfloor,<br>
\end{align*}<br>
\]</span> with <span>\(n\)</span> being an integer and <span>\(y\)</span> any real number. Thus, <span>\[<br>
(n+1) – \lceil n \cdot 0.9\rceil = (n+1) + \lfloor -n \cdot 0.9\rfloor =<br>
\lfloor (n+1) – n \cdot 0.9\rfloor = \lfloor 0.1n+1 \rfloor.<br>
\]</span></p>
<p><strong>Important conclusion</strong>: We have thus shown that with the <code>type=1</code> quantile method we have the duality between having more than 10% critical samples and the 90th percentile of the measurements being above 15 ppm.</p>
<h3>Other Quantile Types</h3>
<p>Since <span>\(\hat{F}\)</span> has jumps of size <span>\(1/n\)</span> the actual value of <span>\(\hat{F}(\hat{x}_{p})\)</span> can end up being somewhat larger than the desired <span>\(p\)</span>. Therefore, <span data-cites="hyndman_fan1996">Hyndman and Fan (1996)</span> prefer estimators interpolating between two adjacent order statistics. Also because such estimators have a lower mean squared error in most cases <span data-cites="dielman_etal1994">(Dielman, Lowry, and Pfaffenberger 1994)</span>. As an example of such an extended estimator, the <code>type=5</code> quantile estimator is defined by letting <span>\(h=p n + 1/2\)</span> and then computing: <span>\[<br>
\hat{x}_p = x_{\lfloor h \rfloor} + (h – \lfloor h \rfloor) (x_{\lfloor h<br>
\rfloor + 1} – x_{\lfloor h \rfloor}).<br>
\]</span></p>
<p>Doing this either manually or using the <code>quantile</code> function one obtains:</p>
<div>
<pre><code>## Small function to manually compute the type=5 p*100th quantile
## for the sample x
manual_quantile_5 &lt;-<span> </span>function(x, p) {
  h &lt;-<span> </span><span>length</span>(x)*p +<span> </span><span>0.5</span>
  x &lt;-<span> </span><span>sort</span>(x)
  <span>return</span>(x[<span>floor</span>(h)] +<span> </span>(h -<span> </span><span>floor</span>(h))*<span> </span>(x[<span>floor</span>(h)+<span>1</span>] -<span> </span>x[<span>floor</span>(h)]))
}

<span>c</span>(<span>&quot;manual.90%&quot;</span>=<span>manual_quantile_5</span>(lead, <span>p=</span><span>0.9</span>), <span>R=</span><span>quantile</span>(lead, <span>probs=</span><span>0.9</span>,<span>type=</span><span>5</span>))</code></pre>
</div>
<pre><code>## manual.90%      R.90% 
##       18.8       18.8</code></pre>
<p>Instead of reading the above or using the R code one can also instead watch a more <a rel="nofollow" href="https://youtu.be/9pql00zr700" target="_blank">didactic whiteboard explanation</a> by Professor <a rel="nofollow" href="http://www.emich.edu/math/faculty/cgardiner.php" target="_blank">Christopher Gardiner</a> for <a rel="nofollow" href="http://michiganradio.org/post/video-how-dropping-two-flints-lead-test-numbers-changed-things-state" target="_blank">Michigan Radio</a> on how to calculate the 90% quantile using a <code>type=5</code> argument for the Flint sample. However, the important point of the above calculations is that this quantile type is <strong>of limited interest</strong>, because the Lead and Copper Rule implicitly defines that one has to use the <code>type=1</code> quantile. To make this point even more explicit, we use sampling with replacement from the Flint data to construct a dataset, where the 90% <code>type=5</code>-quantile is above 15 ppm, but the percentage of samples above the 15 ppm threshold is less than 10%.</p>
<div>
<pre><code>##Function to compute the proportion critical as well as the 90% quantile
##using type (type)-quantiles. Returns the quantile as well as the proportion
##above the 15 ppm threshold
prop_critical_and_q90 &lt;-<span> </span>function(x, <span>type=</span><span>5</span>) {
  q90 &lt;-<span> </span><span>as.numeric</span>(<span>quantile</span>(x, <span>type=</span>type,<span>probs=</span>p))
  prop &lt;-<span> </span><span>mean</span>(x><span>15</span>)
  <span>c</span>(<span>q90=</span>q90,<span>prop=</span> prop)
}

##Make 100 datasets by sampling with replacement
r &lt;-<span> </span><span>data.frame</span>(<span>seed=</span><span>seq_len</span>(<span>100</span>)) %>%<span> </span>rowwise %>%<span> </span><span>do</span>({
  <span>set.seed</span>(.$seed)
  newdata &lt;-<span> </span><span>sample</span>(lead, <span>replace=</span><span>TRUE</span>)
  <span>as.data.frame</span>(<span>seed=</span>.$seed, <span>t</span>(<span>prop_critical_and_q90</span>(newdata)))
})

##Check which datasets violate the duality between quantile and
##percentage above threshold assumption
r %&lt;>%<span> </span><span>mutate</span>(<span>violates_duality =</span>  q90 ><span> </span><span>15</span> &amp;<span> </span>prop <span> </span><span>0.1</span>)

##Do the stats for this dataset
(five &lt;-<span> </span>r %>%<span> </span><span>filter</span>(violates_duality) %>%<span> </span><span>slice</span>(<span>1</span>:<span>5</span>))</code></pre>
</div>
<pre><code>## # A tibble: 5 × 3
##     q90       prop violates_duality
##   </code></pre>