Rule of Thumb

Observations that stray outside of the 95% confidence envelope are statistically significant outliers.

For illustration, R sample data (state.x77) is used in the plot example. This data set give us information such as population size, average income, illiteracy rate, temperature, and murder rate for 50 states in USA.

The qqPlot function from {car} is used to detect outliers. It draws theoretical quantile-comparison plots for variables and for studentized residuals from a linear model. A comparison line is drawn on the plot either through the quartiles of the two distributions, or by robust regression.

library(car)

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])

# fit the data into the linear regression model
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
# summary(fit)

# From our regression model example, we can start
# investigating outliers observation by using Q-Q plot.
car::qqPlot(fit1,labels=row.names(states), id.method="identify",
       simulate=TRUE, main="Q-Q Plot")

##       Nevada Rhode Island 
##           28           39

As you can observed the plot above, Nevada (28th observation) and Rhode Island (39th observation) are states that detected as potential outliers.

Rule of Thumb

Rousseeuw and Leroy (1987) suggest using \(h_{i}>2p/n\) as benchmarks for leverages.

Some researchers believe Hat-values exceeding about twice the average hat-value should be considered noteworthy.

Hat-value can be easily calculated using the hatvalues() function from {stats} package.

Similar to QQ-Plot example, R sample data (state.x77) is also used in the hat-value example.

# data states 
states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])

# fit the data into the linear regression model
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

# fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
# You can compute the high leverage observation by looking at the
# ratio of number of parameters estimated in model and sample size.
# If an observation has a ratio greater than 2-3 times the average ratio,
# then the observation considers as high-leverage points.
# high leverage function

p <- length(coefficients(fit1)) # number of IVs + Intercept
n <- length(fitted(fit1))       # number of observations (N)

# 2-3 times average ratio as cutoff
ratio <- p/n

# display hatvalues
hatvalues(fit1)

##        Alabama         Alaska        Arizona       Arkansas     California 
##     0.09819211     0.43247319     0.12264914     0.08722029     0.34087628 
##       Colorado    Connecticut       Delaware        Florida        Georgia 
##     0.05469964     0.09474658     0.05859637     0.12940379     0.05803142 
##         Hawaii          Idaho       Illinois        Indiana           Iowa 
##     0.22443594     0.06905425     0.09683616     0.03871740     0.04627447 
##         Kansas       Kentucky      Louisiana          Maine       Maryland 
##     0.05245106     0.05131253     0.17022388     0.09966458     0.06942663 
##  Massachusetts       Michigan      Minnesota    Mississippi       Missouri 
##     0.02687166     0.05740920     0.04638897     0.14748152     0.04199088 
##        Montana       Nebraska         Nevada  New Hampshire     New Jersey 
##     0.05253317     0.04421328     0.09508977     0.06387335     0.06616617 
##     New Mexico       New York North Carolina   North Dakota           Ohio 
##     0.15998871     0.23298637     0.05218913     0.10529028     0.08652446 
##       Oklahoma         Oregon   Pennsylvania   Rhode Island South Carolina 
##     0.05117056     0.19057634     0.11185698     0.04562377     0.10445404 
##   South Dakota      Tennessee          Texas           Utah        Vermont 
##     0.07553161     0.04580225     0.13392638     0.06970746     0.08788063 
##       Virginia     Washington  West Virginia      Wisconsin        Wyoming 
##     0.03230467     0.21662901     0.05839687     0.04173326     0.06012359

# High leverage function. 
highleverage <- function(fit) {

  p <- length(coefficients(fit)) # number of IVs + Intercept
  n <- length(fitted(fit))       # number of observations (N)
  ratio <- p/n
  # leverage ratio
  plot(hatvalues(fit), main="Index Plot of Ratio")

  # 2-3 times the average ratio as cutoff lines.
  abline(h=c(2,3)*ratio, col="red", lty=2)

  # identify high leverage points from graph.
  graphics::identify(1:n, hatvalues(fit), names(hatvalues(fit)))
}

# run function 
# Since the graphics::identify function is included, we will not actually run this code in the rmarkdown.
# highleverage(fit1)

The graph above indicates that Alaska and California have the highest ratio. It means these two states are quite unusual if we compared their predictor values compared to other states. Here is why:

• For Alaska, this state has a much higher income compared to other states, while having smaller population and lower temperature.

• For California, this state has much a bigger population than other states, while having higher income and higher temperature.

Both of these states are indeed quite particular compared to the other 48 states. One thing to remember that not all high-leverage observations are influential observations. It will depend if they are outliers or not.

Unusual observations typically have large residuals but not necessarily so-high leverage observations can have small residuals because they pull the line towards them.

The ordinary or simple residuals (observed - predicted values) are the most commonly used measures for detecting outliers.

Rule of Thumb

A large standardized residual (\(d_{i}>3\)) indicates the existence of outliers (Montgomery, D.C., et al.).

However, The numerator and denominator are not independent in the standardized residual. Thus, \(E_{i}^{'}\) does not follow a t-distribution: If \(|E_{i}|\) is large, the standard error is also large:

\[S_{E}=\sqrt{\sum{E_{i}^2}/(n-k-1)}\]

Rule of Thumb

• \(r_{i}>3\) indicates that \(e_{i}\) is an outlier, \(i =1,2,…,n\) (Montgomery, D. C., et al.).

• Also, some researchers suggested that attention should be paid to studentized residuals that exceed +2 or -2 and get even more concerned about residuals that exceed \(|2|\). Observations that have a studentized residual outside the \(\pm2\) range are considered statistically significant at the 95% a level and even yet more concerned about residuals that exceed \(|3|\).

Now lets look at the leverage’s to identify observations that will have potential great influence on regression coefficient estimates. Generally, a point with leverage greater than \(\frac{2k+2}{n}\) should be carefully examined, where \(k\) is the number of predictors (IV) and \(n\) is the number of observations.

Studentized Residuals-the Bonferroni adjustment

Since the furthest outlier are going to be selected, it is not legitimate to use a simple t-test. We would expect that \(5%\) of the studentized residuals would be beyond \(t_{.025}\pm2\) by chance alone. In order to fix this, we can make a Bonferroni adjustment to the p-value. The Bonferroni p-value for the largest outlier is: \(p=2np\)
where \(p\)
is the unadjusted p-value from a t-test with n-k-2 degrees of freedom. A special t-table is needed if you do this calculation manually, but this is done for you automatically, using the outlierTest function in the {car} package.

The Bonferroni-adjusted outlier test in {car} tests the largest absolute studentized residual.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])

# fit the data into the linear regression model
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

# Applying outlierTest function is helping us to
# confirm if potential outliers are indeed outliers.
car::outlierTest(fit1)

##        rstudent unadjusted p-value Bonferroni p
## Nevada 3.542929         0.00095088     0.047544

As you can see the studentized residual and its Bonferroni’s P-value, state-Nevada could be the potential outliers.

In {car} package, there isn’t a single-alone function to plot studentized residual, but through function-influencePlot and influenceIndexPlot, you can incorporate all information from outlier, high-leverage, and influential observations into a single informative plot.

# id.n indicates how many the most unusual observations are identified.
car::influenceIndexPlot(fit1, id.n=3)

In this diagnostic plots, you can find 4 different plots:

• Hat-value
• Bonferroni’s P-value
• Studentdized residual
• Cook’s Distance

The id.n keyword indicates how many the most unusual observations are identified.

Rule of thumb

A standard cut-off for DFBetas is: \(D_{ij}^{*} = 2\times n^{-.5}\).

The R function dfbetas from {stats} package is used to calculate DFBetas value.

# data states 
states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])

# fit the data into the linear regression model
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)


## an influential observation is one that
#  combines discrepancy with leverage.
fit1.dfbetas <- dfbetas(fit1)


# c(2,3) specifies the coefficients of interests (Population and illiteracy)
plot(fit1.dfbetas[,c(2,3)],
     xlim = c(-2,2),ylim=c(-1,1),
     main="DFBetas for the Population and Illiteracy coefficients")

# adds the rule of thumb cut-off line
abline(h=2/sqrt(length(states)), lty=2)

Then, you can identify the influential point above the cut-off line.

Rule of Thumb

Rousseeuw and Leroy (1987) suggest using \(MD_{i}^2>\chi_{p-1;0.95}^2\) as benchmarks for Mahalanobis distances.

Rule of Thumb

The lowest value that Cook’s D can assume is zero (0). There is no significance test for \(D_{i}\) (i.e., the F value here measures only distance) but a cut-off rule of thumb is: \[D_{i}>\frac{4}{n-k-1}\] Generally, when the statistics , \(CD_{i}\), \(h_{i}\) and \(MD_{i}\) are large, case \(i\) may be an outlier or influential case.

Cook’s Distance can be easily calculated by cooks.distance function from {stats} package. You can also directly through plot function, as long as you specify the cook.levels.

Please find the Cook’s D calculation example below.

states <- as.data.frame(state.x77[,c("Murder", "Population",
                                     "Illiteracy", "Income", "Frost")])

# fit the data into the linear regression model
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
##  Influence = Leverage x Discrepancy
# 1. Cook's distance calculation.
fit1_cookd <- cooks.distance(fit1)

cutoff <- 4/(nrow(states)-length(fit1$coefficients)-2)

# plot Cook's Distance
plot(fit1, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")

Cook’s distance, leverages, and Mahalanobis distance can be effective for finding influential cases when a single outlier exist, but can fail if there are two or more outliers.

Nevertheless, these numerical diagnostics combined with plots such as residuals versus fitted values and fitted values versus the response are probably the most effective techniques for detecting cases that affect the fitted values when the multiple linear regression model is a good approximation for the bulk of the data.

We are going to use the same data example and linear regression model (fit1) for the Cook’s distance plot.

The function influencePlot() from {car} package, is a very handy function to represent these unusual observation issues. This function creates a bubble plot of Studentized residuals versus hat values, with the areas of the circles representing the observations proportional to the value Cook’s distance.

• Vertical reference lines are drawn at twice and three times
• the average hat value, horizontal reference lines at -2, 0, and 2 on the Studentized-residual scale.
• x-axis: Hat-value (with cutoffs)
• Y-axis: studentized residual
• size of bubble (Cook’s D)

# id.n indicates how many most unusual observations are identified.
car::influencePlot(fit1, id.n=3)

##                 StudRes        Hat       CookD
## Alaska        1.7536917 0.43247319 0.448050997
## California   -0.2761492 0.34087628 0.008052956
## Nevada        3.5429286 0.09508977 0.209915743
## Rhode Island -2.0001631 0.04562377 0.035858963

Rule of Thumb

Some analysts suggest investigating observations for which \(|DFITS_{i}|\) is greater than \(2\sqrt{(p+1)/(n−p−1)}\). Cook’s D and DFITS give similar answers.