<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/percentile-calculations-in-water-quality-regulations/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://r.prevos.net/percentile-calculations/"> The Devil is in the Data</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><img sizes="(max-width: 300px) 100vw, 300px" data-image-description="" data-orig-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/percentiles.png?fit=1571%2C962" data-orig-size="1571,962" data-image-title="percentile calculation" class="wp-image-363" data-comments-opened="1" data-recalc-dims="1" data-attachment-id="363" title="Percentile Calculations in Water Quality Regulations" data-permalink="http://r.prevos.net/percentile-calculations/percentiles/" alt="Percentile Calculations in Water Quality Regulations" data-medium-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/percentiles.png?fit=300%2C184" data-large-file="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/percentiles.png?fit=584%2C358" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i2.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/percentiles.png?resize=300%2C184">Percentile calculations can be more tricky than at first meets the eye. A <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Percentile">percentile</a> indicates the value below which a percentage of observations fall. Some percentiles have special names, such as the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Quartile">quartile</a> or the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Decile">decile</a>, both of which are <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Quantile">quantiles</a>. This deceivingly simple definition hides the various ways to determine this number. Unfortunately, there is <a target="_blank" rel="nofollow" href="https://www.amherst.edu/media/view/129116/original/Sample+Quantiles.pdf">no standard definition for percentiles</a>, so which method do you use?</p>
<p>The <a target="_blank" rel="nofollow" href="http://stat.ethz.ch/R-manual/R-devel/library/stats/html/quantile.html">quantile</a> function in R generates sample percentiles corresponding to the given probabilities. By default, the quantile function provides the quartiles and the minimum and maximum values. The code snippet below generates semi-random data, plots the histogram and visualises the third quartile.</p>
<pre>
set.seed(1969)
test.data &lt;- rnorm(n=10000, mean=100, sd=15)
library(ggplot2)
ggplot(as.data.frame(test.data), aes(test.data)) + 
 geom_histogram(binwidth=1, aes(y=..density..), fill=&quot;dodgerblue&quot;) + 
 geom_line(stat=&quot;function&quot;, fun=dnorm, args=list(mean=100, sd=15), colour=&quot;red&quot;, size=1) + 
 geom_area(stat=&quot;function&quot;, fun=dnorm, args=list(mean=100, sd=15), 
 colour=&quot;red&quot;, fill=&quot;red&quot;, alpha=0.5, xlim=quantile(test.data, c(0.5, 0.75))) + 
 theme(text=element_text(size=16))
</pre>
<p>The quantile default function and the 95<sup>th</sup> percentile give the following results:</p>
<pre>
> quantile(test.data)
       0%       25%       50%       75%      100% 
 39.91964  89.68041 100.16437 110.01910 153.50195 

> quantile(test.data, probs=0.95)
     95% 
124.7775 
</pre>
<h2>Methods of percentile calculations</h2>
<p>The quantile function in R provides for nine different ways to calculate percentiles. Each of these options uses a different method to interpolate between observed values. I will not discuss the mathematical nuances between these methods. <a target="_blank" rel="nofollow" href="https://www.amherst.edu/media/view/129116/original/Sample+Quantiles.pdf">Hyndman and Fan</a> (1996) provide a detailed overview of these methods.</p>
<p>The differences between the nine available methods only matter in skewed distributions, such as water quality data. For the normal distribution simulated above the outcome for all methods is exactly the same, as illustrated by the following code.</p>
<pre>
> sapply(1:9, function(m) quantile(test.data, 0.95, type=m))

     95%      95%      95%      95%      95%      95%      95%      95%      95% 
124.7775 124.7775 124.7775 124.7775 124.7775 124.7775 124.7775 124.7775 124.7775 
</pre>
<h2>Percentile calculations in water quality</h2>
<p>The <a target="_blank" rel="nofollow" href="https://www.nhmrc.gov.au/_files_nhmrc/file/publications/nhmrc_adwg_6_version_3.3_2.pdf">Australian Drinking Water Quality Guidelines</a> (November 2016) specify that: “based on aesthetic considerations, the turbidity should not exceed 5 NTU at the consumer’s tap”. The Victorian Safe Drinking Water Regulations (2015) relax this requirement and require that:</p>
<blockquote>
<p>“The 95<sup>th</sup> percentile of results for samples in any 12 month period must be less than or equal to 5.0 NTU.”</p>
</blockquote>
<p>The Victorian regulators <a target="_blank" rel="nofollow" href="https://www2.health.vic.gov.au/Api/downloadmedia/%7BA1F6D255-D5C7-4B7E-AAE5-8B7451EDE81A%7D">also specify </a>that the percentile should be calculated with the <em>Weibull Method</em>. This requirement raises two questions: What is the Weibull method? How do you implement this requirement in R?</p>
<p>The term Weibull Method is a bit confusing as this is not a name used by statisticians. In Hyndman &amp; Fan (1996), this method has the less poetic name <i>$latex Q_8 p &amp;s=1</i>. <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Waloddi_Weibull">Waloddi Weibull</a>, a Swedish engineer famous for his distribution, was one of the first to describe this method. Only the regulator in Victoria uses that name, which is based on <a target="_blank" rel="nofollow" href="http://amzn.to/2k8shr8">McBride (2005)</a>. This theoretical interlude aside, how can we practically apply this to water quality data?</p>
<p>In case you are interested in how the Weibull method works, the <em>weibull.quantile</em> function shown below calculates a quantile <em>p</em> for a vector <em>x</em> using this method. This function gives the same result as <em>quantile(x, p, type=6)</em>.</p>
<pre>
weibull.quantile &lt;- function(x, p) {
    # Order Samples from large to small
    x &lt;- x[order(x, decreasing=FALSE)]
    # Determine ranking of percentile according to Weibull (1939)
    r &lt;- p * (length(x) + 1)
    # Linear interpolation
    rfrac &lt;- (r - floor(r))
    return((1 - rfrac) * x[floor(r)] + rfrac * x[floor(r) + 1])
}
</pre>
<h2>Turbidity Data Example</h2>
<p>Turbidity data is not normally distributed as it is always larger than zero. In this example, the turbidity results for the year 2016 for the water system in <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Tarnagulla">Tarnagulla </a>are used to illustrate the percentile calculations. The range of weekly turbidity measurements is between 0.,05 NTU and 0.8 NTU, well below the aesthetic limits.</p>
<div class="wp-caption"><img sizes="(max-width: 584px) 100vw, 584px" data-image-description="" data-orig-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/turbidity.png?fit=2150%2C912" data-orig-size="2150,912" data-image-title="turbidity" data-comments-opened="1" class="wp-image-379" data-recalc-dims="1" data-attachment-id="379" data-permalink="http://r.prevos.net/percentile-calculations/turbidity/" alt="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/turbidity.png?resize=584%2C248" data-medium-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/turbidity.png?fit=300%2C127" data-large-file="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/turbidity.png?fit=584%2C248" data-image-meta="{&quot;aperture&quot;:&quot;0&quot;,&quot;credit&quot;:&quot;&quot;,&quot;camera&quot;:&quot;&quot;,&quot;caption&quot;:&quot;&quot;,&quot;created_timestamp&quot;:&quot;0&quot;,&quot;copyright&quot;:&quot;&quot;,&quot;focal_length&quot;:&quot;0&quot;,&quot;iso&quot;:&quot;0&quot;,&quot;shutter_speed&quot;:&quot;0&quot;,&quot;title&quot;:&quot;&quot;,&quot;orientation&quot;:&quot;0&quot;}" src="https://i1.wp.com/r.prevos.net/wp-content/uploads/sites/11/2017/02/turbidity.png?resize=584%2C248">
<p class="wp-caption-text">Turbidity at customer tap for each zone in the Tarnagulla system in 2016 (n=53).</p>
</div>
<p>When we calculate the percentiles for all nine methods available in the base-R function we see that the so-called Weibull method generally provides the most conservative result.</p>
<table>
<thead>
<tr>
<th>Zone</th>
<th>R1</th>
<th>R2</th>
<th>R3</th>
<th>R4</th>
<th>R5</th>
<th>R6</th>
<th>R7</th>
<th>R8</th>
<th>R9</th>
</tr>
</thead>
<tbody>
<tr>
<td>Bealiba</td>
<td>0.300</td>
<td>0.300</td>
<td>0.200</td>
<td>0.240</td>
<td>0.290</td>
<td>0.300</td>
<td>0.245</td>
<td>0.300</td>
<td>0.300</td>
</tr>
<tr>
<td>Dunolly</td>
<td>0.40000</td>
<td>0.40000</td>
<td>0.30000</td>
<td>0.34000</td>
<td>0.39000</td>
<td>0.43500</td>
<td>0.34500</td>
<td>0.40500</td>
<td>0.40125</td>
</tr>
<tr>
<td>Laanecoorie</td>
<td>0.50000</td>
<td>0.50000</td>
<td>0.40000</td>
<td>0.44000</td>
<td>0.49000</td>
<td>0.53500</td>
<td>0.44500</td>
<td>0.50500</td>
<td>0.50125</td>
</tr>
<tr>
<td>Tarnagulla</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
<td>0.4</td>
</tr>
</tbody>
</table>
<p></p>
<p>The graph and the table were created with the following code snippet:</p>
<pre>
ggplot(turbidity, aes(Result)) + 
 geom_histogram(binwidth=.05, fill=&quot;dodgerblue&quot;, aes(y=..density..)) + 
 facet_wrap(~Zone) + 
 theme(text=element_text(size=16))

tapply(turbidity$Result, turbidity$Zone, 
 function(x) sapply(1:9, function(m) quantile(x, 0.95, type=m)))
</pre>
<p>The post <a target="_blank" rel="nofollow" href="http://r.prevos.net/percentile-calculations/">Percentile Calculations in Water Quality Regulations</a> appeared first on <a target="_blank" rel="nofollow" href="http://r.prevos.net/">The Devil is in the Data</a>.</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://r.prevos.net/percentile-calculations/"> The Devil is in the Data</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/Mt66UgVFwlU"></a></div>