<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/outlier-detection-and-treatment-with-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://datascienceplus.com/outlier-detection-and-treatment-with-r/"> DataScience+</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Outliers in data can distort predictions and affect the accuracy, if you don’t detect and handle them appropriately especially in regression models.</p>
<h2>Why outliers treatment is important?</h2>
<p>Because, it can drastically bias/change the fit estimates and predictions. Let me illustrate this using the cars dataset.</p>
<p>To better understand the implications of outliers better, I am going to compare the fit of a simple linear regression model on cars dataset with and without outliers. In order to distinguish the effect clearly, I manually introduce extreme values to the original cars dataset. Then, I predict on both the datasets.</p>
<pre>
# Inject outliers into data.
cars1 &lt;- cars[1:30, ]  # original data
cars_outliers &lt;- data.frame(speed=c(19,19,20,20,20), dist=c(190, 186, 210, 220, 218))  # introduce outliers.
cars2 &lt;- rbind(cars1, cars_outliers)  # data with outliers.

# Plot of data with outliers.
par(mfrow=c(1, 2))
plot(cars2$speed, cars2$dist, xlim=c(0, 28), ylim=c(0, 230), main=&quot;With Outliers&quot;, xlab=&quot;speed&quot;, ylab=&quot;dist&quot;, pch=&quot;*&quot;, col=&quot;red&quot;, cex=2)
abline(lm(dist ~ speed, data=cars2), col=&quot;blue&quot;, lwd=3, lty=2)

# Plot of original data without outliers. Note the change in slope (angle) of best fit line.
plot(cars1$speed, cars1$dist, xlim=c(0, 28), ylim=c(0, 230), main=&quot;Outliers removed \n A much better fit!&quot;, xlab=&quot;speed&quot;, ylab=&quot;dist&quot;, pch=&quot;*&quot;, col=&quot;red&quot;, cex=2)
abline(lm(dist ~ speed, data=cars1), col=&quot;blue&quot;, lwd=3, lty=2)
</pre>
<p>Plot of original data without outliers:<br>
<a title="" data-rel="lightbox-0" target="_blank" rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2016/12/outliers_effect.png"><img sizes="(max-width: 490px) 100vw, 490px" alt="outliers effect on model fit" class="wp-image-4941" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2016/12/outliers_effect-490x245.png?resize=456%2C228" data-recalc-dims="1"></a></p>
<p>Notice the change in slope of the best fit line after removing the outliers. Had we used the outliers to train the model(left chart), our predictions would be exagerated (high error) for larger values of speed because of the larger slope.</p>
<h2>Detect Outliers </h2>
<p><strong>Univariate approach</strong><br>
For a given continuous variable, outliers are those observations that lie outside 1.5 * IQR, where IQR, the ‘Inter Quartile Range’ is the difference between 75th and 25th quartiles. Look at the points outside the whiskers in below box plot.</p>
<pre>outlier_values &lt;- boxplot.stats(inputData$pressure_height)$out  # outlier values.
boxplot(inputData$pressure_height, main=&quot;Pressure Height&quot;, boxwex=0.1)
mtext(paste(&quot;Outliers: &quot;, paste(outlier_values, collapse=&quot;, &quot;)), cex=0.6)</pre>
<h3>Bivariate approach</h3>
<p>Visualize in box-plot of the X and Y, for categorical X’s</p>
<pre>url &lt;- &quot;http://rstatistics.net/wp-content/uploads/2015/09/ozone.csv&quot;
ozone &lt;- read.csv(url)

# For categorical variable
boxplot(ozone_reading ~ Month, data=ozone, main=&quot;Ozone reading across months&quot;)  # clear pattern is noticeable.
boxplot(ozone_reading ~ Day_of_week, data=ozone, main=&quot;Ozone reading for days of week&quot;)  # this may not be significant, as day of week variable is a subset of the month var.</pre>
<p>Box-plot:<br>
<a title="" data-rel="lightbox-1" target="_blank" rel="nofollow" href="https://i0.wp.com/datascienceplus.com/wp-content/uploads/2016/12/bivariate_boxplot.png"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-4942" src="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2016/12/bivariate_boxplot-490x245.png?resize=456%2C228" data-recalc-dims="1"></a></p>
<p>What is the inference? The change in the level of boxes suggests that Month seem to have an impact in ozone_reading while Day_of_week does not. Any outliers in respective categorical level show up as dots outside the whiskers of the boxplot.</p>
<pre># For continuous variable (convert to categorical if needed.)
boxplot(ozone_reading ~ pressure_height, data=ozone, main=&quot;Boxplot for Pressure height (continuos var) vs Ozone&quot;)
boxplot(ozone_reading ~ cut(pressure_height, pretty(inputData$pressure_height)), data=ozone, main=&quot;Boxplot for Pressure height (categorial) vs Ozone&quot;, cex.axis=0.5)</pre>
<p><a title="" data-rel="lightbox-2" target="_blank" rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2016/12/boxplot_continuous.png"><img sizes="(max-width: 490px) 100vw, 490px" alt="Boxplot Continuous" class="wp-image-4943" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2016/12/boxplot_continuous-490x245.png?resize=456%2C228" data-recalc-dims="1"></a></p>
<p>You can see few outliers in the box plot and how the <code>ozone_reading</code> increases with <code>pressure_height</code>. Thats clear.</p>
<h3>Multivariate Model Approach</h3>
<p><strong>Cook’s Distance</strong><br>
Cook’s distance is a measure computed with respect to a given regression model and therefore is impacted only by the X variables included in the model. But, what does cook’s distance mean? It computes the influence exerted by each data point (row) on the predicted outcome.</p>
<p>The cook’s distance for each observation i measures the change in Ŷ Y^ (fitted Y) for all observations with and without the presence of observation i, so we know how much the observation i impacted the fitted values.</p>
<pre>mod &lt;- lm(ozone_reading ~ ., data=ozone)
cooksd &lt;- cooks.distance(mod)</pre>
<p><strong>Influence measures</strong><br>
In general use, those observations that have a cook’s distance greater than 4 times the mean may be classified as influential. This is not a hard boundary.</p>
<pre>plot(cooksd, pch=&quot;*&quot;, cex=2, main=&quot;Influential Obs by Cooks distance&quot;)  # plot cook's distance
abline(h = 4*mean(cooksd, na.rm=T), col=&quot;red&quot;)  # add cutoff line
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),&quot;&quot;), col=&quot;red&quot;)  # add labels</pre>
<p>Influential Obs by Cooks distance:<br>
<a title="" data-rel="lightbox-3" target="_blank" rel="nofollow" href="https://i2.wp.com/datascienceplus.com/wp-content/uploads/2016/12/cooks_distance_plot.png"><img sizes="(max-width: 490px) 100vw, 490px" alt="" class="wp-image-4963" src="https://i1.wp.com/datascienceplus.com/wp-content/uploads/2016/12/cooks_distance_plot-490x245.png?resize=456%2C228" data-recalc-dims="1"></a></p>
<p>Now lets find out the influential rows from the original data. If you extract and examine each influential row 1-by-1 (from below output), you will be able to reason out why that row turned out influential. It is likely that one of the X variables included in the model had extreme values.</p>
<pre>
influential  4*mean(cooksd, na.rm=T))])  # influential row numbers
head(ozone[influential, ])  # influential observations.<em>
#> Month Day_of_month Day_of_week ozone_reading pressure_height Wind_speed Humidity
#> 19 1 19 1 4.07 5680 5 73
#> 23 1 23 5 4.90 5700 5 59
#> 58 2 27 5 22.89 5740 3 47
#> 133 5 12 3 33.04 5880 3 80
#> 135 5 14 5 31.15 5850 4 76
#> 149 5 28 5 4.82 5750 3 76
#> Temperature_Sandburg Temperature_ElMonte Inversion_base_height Pressure_gradient
#> 19 52 56.48 393 -68
#> 23 69 51.08 3044 18
#> 58 53 58.82 885 -4
#> 133 80 73.04 436 0
#> 135 78 71.24 1181 50
#> 149 65 51.08 3644 86
#> Inversion_temperature Visibility
#> 19 69.80 10
#> 23 52.88 150
#> 58 67.10 80
#> 133 86.36 40
#> 135 79.88 17
#> 149 59.36 70</em></pre>
<p>Lets examine the first 6 rows from above output to find out why these rows could be tagged as influential observations.</p>
<ul>
<li>Row 58, 133, 135 have very high <code>ozone_reading</code>.</li>
<li>Rows 23, 135 and 149 have very high <code>Inversion_base_height</code>.</li>
<li>Row 19 has very low <code>Pressure_gradient</code>.</li>
</ul>
<h2>Outliers Test</h2>
<p>The function outlierTest from car package gives the most extreme observation based on the given model.</p>
<pre>car::outlierTest(mod)
#> No Studentized residuals with Bonferonni p  Largest |rstudent|:
#>     rstudent unadjusted p-value Bonferonni p
#> 243 3.045756          0.0026525      0.53845</pre>
<p>This output suggests that observation in row 243 is most extreme.</p>
<h2>Outliers package</h2>
<p>The outliers package provides a number of useful functions to systematically extract outliers. Some of these are convenient and come handy, especially the <code>outlier()</code> and <code>scores()</code> functions.</p>
<p><strong>Outliers</strong><br>
<code>outliers</code> gets the extreme most observation from the mean. If you set the argument opposite=TRUE, it fetches from the other side.</p>
<pre>set.seed(1234)
y=rnorm(100)
outlier(y)<em>
#> [1] 2.548991</em>
outlier(y,opposite=TRUE)<em>
#> [1] -2.345698</em>
dim(y) &lt;- c(20,5)  # convert it to a matrix
outlier(y)<em>
#> [1] 2.415835 1.102298 1.647817 2.548991 2.121117</em>
outlier(y,opposite=TRUE)<em>
#> [1] -2.345698 -2.180040 -1.806031 -1.390701 -1.372302</em></pre>
<p><strong>Scores</strong><br>
There are two aspects the the <code>scores()</code> function.<br>
Compute the normalised scores based on “z”, “t”, “chisq” etc<br>
Find out observations that lie beyond a given percentile based on a given score.</p>
<pre>set.seed(1234)
x = rnorm(10)
scores(x)  # z-scores => (x-mean)/sd
scores(x, type=&quot;chisq&quot;)  # chi-sq scores => (x - mean(x))^2/var(x)<em>
#> [1] 0.68458034 0.44007451 2.17210689 3.88421971 0.66539631 . . .</em>
scores(x, type=&quot;t&quot;)  # t scores
scores(x, type=&quot;chisq&quot;, prob=0.9)  # beyond 90th %ile based on chi-sq<em>
#> [1] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE</em>
scores(x, type=&quot;chisq&quot;, prob=0.95)  # beyond 95th %ile
scores(x, type=&quot;z&quot;, prob=0.95)  # beyond 95th %ile based on z-scores
scores(x, type=&quot;t&quot;, prob=0.95)  # beyond 95th %ile based on t-scores</pre>
<h2>Treatment</h2>
<p>Once the outliers are identified, you may rectify it by using one of the following approaches.</p>
<p><strong>Imputation</strong><br>
Imputation with mean / median / mode. This method has been dealt with in detail in the discussion about <a target="_blank" rel="nofollow" href="http://datascienceplus.com/missing-value-treatment/">treating missing values</a>. Another robust method which we covered at DataScience+ is <a target="_blank" rel="nofollow" href="http://datascienceplus.com/handling-missing-data-with-mice-package-a-simple-approach/">multivariate imputation by chained equations</a>.</p>
<p><strong>Capping</strong><br>
For missing values that lie outside the 1.5 * IQR limits, we could cap it by replacing those observations outside the lower limit with the value of 5th %ile and those that lie above the upper limit, with the value of 95th %ile. Below is a sample code that achieves this.</p>
<pre>x &lt;- ozone$pressure_height
qnt &lt;- quantile(x, probs=c(.25, .75), na.rm = T)
caps &lt;- quantile(x, probs=c(.05, .95), na.rm = T)
H &lt;- 1.5 * IQR(x, na.rm = T)
x[x &lt; (qnt[1] - H)]  (qnt[2] + H)] &lt;- caps[2]</pre>
<p><strong>Prediction</strong><br>
In yet another approach, the outliers can be replaced with missing values <code>NA</code> and then can be predicted by considering them as a response variable. We already discussed how to <a target="_blank" rel="nofollow" href="http://datascienceplus.com/missing-value-treatment/">predict missing values</a>.</p>
<blockquote><p>If you liked this post, you might find my video courses <a target="_blank" rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/introduction-r-programming-video">Introduction to R Programming</a> and <a target="_blank" rel="nofollow" href="https://www.packtpub.com/big-data-and-business-intelligence/mastering-r-programming-video">Mastering R Programming</a> or to visit <a target="_blank" rel="nofollow" href="http://r-statistics.co/Outlier-Treatment-With-R.html">My Blog</a>.</p>
</blockquote>
<p>Happy Learning!</p>
<ul></ul>
<p><strong>
<p>Related Post</p>
<p>
<ol>
<li><a title="Permanent Link toR for Publication: Lesson 6, Part 2 – Linear Mixed Effects Models" target="_blank" rel="nofollow" href="http://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-6-part-2-linear-mixed-effects-models-lmem/">R for Publication: Lesson 6, Part 2 – Linear Mixed Effects Models</a></li>
<li><a title="Permanent Link toR for Publication: Lesson 6, Part 1 – Linear Mixed Effects Models" target="_blank" rel="nofollow" href="http://datascienceplus.com/r-for-publication-by-page-piccinini-lesson-6-part-1-linear-mixed-effects-models-lmem/">R for Publication: Lesson 6, Part 1 – Linear Mixed Effects Models</a></li>
<li><a title="Permanent Link toCross-Validation: Estimating Prediction Error" target="_blank" rel="nofollow" href="http://datascienceplus.com/cross-validation-estimating-prediction-error/">Cross-Validation: Estimating Prediction Error</a></li>
<li><a title="Permanent Link toInteractive Performance Evaluation of Binary Classifiers" target="_blank" rel="nofollow" href="http://datascienceplus.com/interactive-performance-evaluation-of-binary-classifiers/">Interactive Performance Evaluation of Binary Classifiers</a></li>
<li><a title="Permanent Link toPredicting wine quality using Random Forests" target="_blank" rel="nofollow" href="http://datascienceplus.com/predicting-wine-quality-using-random-forests/">Predicting wine quality using Random Forests</a></li>
</ol></strong>
<div>
<a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=FBlus9jLZPQ:jdEHqHGW_VE:D7DqB2pKExk"><img src="http://feeds.feedburner.com/~ff/datascienceplus?i=FBlus9jLZPQ:jdEHqHGW_VE:D7DqB2pKExk" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=FBlus9jLZPQ:jdEHqHGW_VE:qj6IDK7rITs"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=qj6IDK7rITs" border="0"></a> <a target="_blank" rel="nofollow" href="http://feeds.feedburner.com/~ff/datascienceplus?a=FBlus9jLZPQ:jdEHqHGW_VE:yIl2AUoC8zA"><img src="http://feeds.feedburner.com/~ff/datascienceplus?d=yIl2AUoC8zA" border="0"></a>
</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://datascienceplus.com/outlier-detection-and-treatment-with-r/"> DataScience+</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/ZEvHgmMVKMg"></a></div>