Doing Residual Analysis Post Regression in R
In this post, we take a deep dive into the R language by exploring residual analysis and visualizing the results with R.
Join the DZone community and get the full member experience.Join For Free
Residuals are essentially gaps that are left when a given model, in this case, linear regression, does not fit the given observations completely.
A close analogy for residual analysis is found in medical pathology. What remains post-metabolism usually becomes an indicator of what was processed and absorbed versus what was not so.
Residuals are the difference between the original value to be modeled and the estimator of original value that came as a result of your model.
Residual error = Y- Y-hat, where Y is the original value and Y-hat is the computed value.
I would expect this error to be as close to zero as possible and will be interested in finding any outliers or aberrations thereof.
A quick way to find outliers is to use either standardized or studentized residuals. The first is simply scaling the residuals relative to their standard deviation and, hence, anything that lies beyond 3 standard deviations can be considered to be an aberration.
## standardize residuals - scale residuals relative to their standard deviation residueStandard<-rstandard(lmfit) df[residueStandard>3,]
Here is what you get:
days.instant days.atemp days.hum days.windspeed days.casual 442 442 0.505046 0.755833 0.110704 3155 456 456 0.421708 0.738333 0.250617 2301 463 463 0.426129 0.254167 0.274871 3252 470 470 0.487996 0.502917 0.190917 2795 471 471 0.573875 0.507917 0.225129 2846 505 505 0.566908 0.456250 0.083975 3410 512 512 0.642696 0.732500 0.198992 2855 513 513 0.641425 0.697083 0.215171 3283 533 533 0.594708 0.504167 0.166667 2963 624 624 0.585867 0.501667 0.247521 3160 645 645 0.538521 0.664167 0.268025 3031 659 659 0.472842 0.572917 0.117537 2806
Naturally, I would expect my model to be unbiased, at least in intention, and hence any leftovers on either side of the regression line that did not make it on the line are expected to be random, i.e. without any particular pattern.
That is, I expect my residual error distributions to follow a bland, normal distribution.
In R, you can do this elegantly with just two lines of code.
1. Plot a histogram of residuals
2. Add a Quantile-Quantile plot with a line that passes through, namely, the first and third quantiles.
So we know that the plot deviates from normal (represented by the straight line).
But that binary aspect of information is seldom enough. So, we should check the skewness and kurtosis to know the dispersion of the distribution.
First, we will compute the skewness; we will use a simple Galton's skewness formula.
## check on symmetry of distribution: skewness summary<-summary(lmfit$residuals) Q1<-summary[] Q2<-summary[] Q3<-summary[] skewness<-((Q1-Q2)+(Q3-Q2))/(Q3-Q1) skewness  0.2665384
This quantifies which side the distribution has swayed to (the right side or left side, or if it's perfectly symmetrical or even).
Another measuring aspect is Kurtosis, which will tell you whether the distribution is assembled towards the center (+ve value) or away from it. (-ve value). Kurtosis is a metric that quantifies the possibility of outliers. A positive value indicates their existence.
kurtosis(lmfit$residuals)  2.454145
Next, apart from visualizing the distribution of residual errors, I would also like to find whether the residuals are correlated with each other.
For a model to account for all the variation in observations, residual errors must occur at random and be uncorrelated with each other.
Durbin-Watson test allows you to test the independence of residuals from each other.
dwt(lm(df$days.windspeed~df$days.atemp+df$days.hum),method='resample', alternative='two.sided') lag Autocorrelation D-W Statistic p-value 1 0.2616792 1.475681 0 Alternative hypothesis: rho != 0
The statistic value needs to be compared with lower and upper limits established by the critical values table, for a given degree of freedom and number of observations. In this case, the range is
Since the computed D-W statistic falls below the lower value of this range, we reject the null hypothesis that residuals are not correlated. Hence, an alternate hypothesis that there may be a correlation between residuals comes into play.
Visually, this can be learned by studying how a model falls short when attempting to capture increases in the value of the original Y value. As this happens, the residuals appear proportional to Y, as Y increases.
Compare this with plotting the fitted Y-hat value against Y. Where Y-hat values tend to fall behind, residuals appear to run along with Y and hence the past residual value appears to be a better predictor of present values as they appear to continue along a fixed slope value. And correlation shows up in Darwin-Watson tests that try to compare and contrast the sum square of the difference between residual errors and their previous value(s) in relation to the sum of a given residual error over all observations.
Correlation between residual errors is thus like failing on a cardio test. If your model simply can't keep up with the rapid changes in Y, it will start looking more and more like an exasperated runner, with the count of lost steps appearing to be more and more correlated than the actual steps taken by the runner.
Reconstructing the subconscious influence your linear regression model was under.
It is important to remember that when you run a linear regression on any given set of observations, each value of THE computed estimator (symbolically represented as Y-hat) of THE dependent variable (similarly represented as Y), is not only dependant on the current value (say, with the observation), but also dependant on each observation.
After all, your multiple linear regression models computed the coefficients for each contributing factor, only after minimizing the error between the actual value of dependent variable Y and the computed estimator of Y (Y-hat) on all observations.
Mathematically, this can be represented by the simple weighted average model as follows.
Yhat(j)= w1j*Y(1)+w2j*Y(2)+w3j*Y(3)+w4j*Y(4)+.....+wnj*Y(n); where n is number of observations and j is one particular estimated value.
If you want to think of this in matrix notation, matrix W will have dimensions [n,n] where n is the number of observations. Similarly, Y represented in matrix form will have dimensions [n,1], and hence Y-hat[j] is represented as a matrix because it is obtained as a product of two matrices (W and Y) and will have the dimensions [n,1], which are the same as Y.
So, each value of the computed estimator can be expressed as a weighted average of all original Y values, indicating to what extent a given computed estimator is influenced by each original Y value.
To put it more simply, if I see a bunch of observations scattered here and there, and I have to somehow draw a straight line that passes through them with as little bias as possible, that is, if I must leave out some observations, I must do so equally on either side of this regression line.
But while doing so, I (or the regression program of choice) subconsciously estimated to what extent the slope of my regression line was influenced by each observation, irrespective of where it belonged.
All I am doing now is back calculating such influences in retrospective. How much weight did each observation carry in order to sway my regression line to either include them or exclude them?
The statistical term for weights is hat values because they bridge computed estimators and their original counterparts.
This is how it is computed in R.
First, we get the model in a matrix form. Then we compute the hat values.
Alternatively, you can use the following function to get a similar result.
Let's think about the limits that can be imposed on each weight.
Obviously, the lowest possible value for the weight will be equal to the possibility when all original Y values contribute equally (since they must contribute something for the linear regression program to estimate the coefficients by working and optimizing through all observations).
In this case, the lower limit of its range will be 1/n, where n is the total number of observations. In any case, since n will always be +ve, the weights will always be <1.
Now try to sum up your hat values and you will see an interesting result...
sum(hatvalues)  4
They equal the number of contributing factors your linear regression model has considered for computation +1.
For instance, in the example dataset, we have three contributing factors, namely: temperature, humidity, and wind speed.
lmfit$coefficients (Intercept) df$days.atemp df$days.hum df$days.windspeed 490.2956 2312.0163 -875.6955 -992.1411
Next, how do we find the most significant or influencing observations?
An elegant way will be to just:
segment hat values into quantiles .
apply the 95% criterion to filter the most aberrant values.
apply this filtering criterion back on observations.
R allows you to all of this in one step!
We can get specific observations that satisfy this criterion.
days.instant days.atemp days.hum days.windspeed days.casual 9 9 0.1161750 0.434167 0.361950 54 21 21 0.1578330 0.457083 0.353242 75 22 22 0.0790696 0.400000 0.171970 93 26 26 0.2036000 0.862500 0.293850 34 32 32 0.2345300 0.829565 0.053213 47 36 36 0.2430580 0.929167 0.161079 100 45 45 0.3983500 0.375833 0.417908 208 50 50 0.3914040 0.187917 0.507463 532 65 65 0.3662520 0.948261 0.343287 114 69 69 0.3856680 0.000000 0.261877 46 90 90 0.2575750 0.918333 0.217646 179 94 94 0.5429290 0.426250 0.385571 734 95 95 0.3983500 0.642083 0.388067 167 106 106 0.4254920 0.888333 0.340808 121 153 153 0.6439420 0.305000 0.292287 736 239 239 0.6355560 0.850000 0.375617 226 249 249 0.5152000 0.886957 0.343943 204 293 293 0.4665250 0.636250 0.422275 471 302 302 0.2279130 0.882500 0.351371 57 341 341 0.4002460 0.970417 0.266175 50 368 368 0.1262750 0.441250 0.365671 89 383 383 0.2752540 0.443333 0.415429 109 388 388 0.2430580 0.911250 0.110708 145 408 408 0.1016580 0.464583 0.409212 73 421 421 0.2556750 0.395833 0.421642 317 433 433 0.5246040 0.567500 0.441563 486 434 434 0.3970830 0.407083 0.414800 447 463 463 0.4261290 0.254167 0.274871 3252 465 465 0.4766380 0.317500 0.358196 905 478 478 0.3895040 0.835417 0.344546 120 543 543 0.5947040 0.373333 0.347642 1077 627 627 0.5650670 0.872500 0.357587 371 667 667 0.4677710 0.694583 0.398008 998 668 668 0.4394000 0.880000 0.358200 2 694 694 0.2487420 0.404583 0.376871 532 722 722 0.2361130 0.441250 0.407346 205 726 726 0.2203330 0.823333 0.316546 9
Finally, you can visualize them with color coding.
That's it! This is a simple overview of doing residual analysis and why it's important.
The link to the data set and model assumptions are mentioned in the references section.
Opinions expressed by DZone contributors are their own.