Boxplot

Boxplots are great graphs to visualize outliers. Most of them identify outliers when those observations lie outside \(1.5 * IQR\). The IQR (Inter Quartile Range) is the value difference between 75th and 25th quartiles. In a boxplot, we will look at the points outside the whiskers.

data(algae)
outliers<-boxplot.stats(algae$Cl)$out
boxplot(algae$Cl, main="Pressure Height", boxwex=0.1)
mtext(paste("Outliers: ", paste(outliers, collapse=", ")), cex=0.5,at = c(0.8,300),outer = F)

calculate_outlier <- function(x) {
  return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}

data("algae")
outlier.data<-  dplyr::mutate(algae,obs=seq_along(algae$season))%>%
                gather(key = "variables",value = "values",4:18,na.rm = T)%>%
                dplyr::group_by(variables)%>%
                dplyr::mutate(outlier=ifelse(calculate_outlier(values),obs,NA))

ggplot()+
  geom_boxplot(outlier.data,mapping=aes(x=variables,y=values))+
  geom_text_repel(outlier.data,mapping=aes(x=variables,y=values,label=outlier))+
  facet_wrap(variables~.,scales = "free")+
  theme_bw()

geom_text_repel does not work as expected, so we use geom_text and move the text to the right to have some space between the points.

p=ggplot(outlier.data)+
  geom_boxplot(mapping=aes(x=variables,y=values))+
  geom_text(mapping=aes(x=variables,y=values,label=outlier),nudge_x = 0.1)+
  facet_wrap(variables~.,scales = "free")+
  theme_bw()
plotly::ggplotly(p)

Cook’s Distance

Cook’s distance is a mathematical measure of the impact of deleting a case on the estimated regression coefficients not a statistical test. Each case or sample has one Cook’s distance \(D\) value, and the plot orders them in terms of their influence on the fitted \(hat_y\).

It is generally helpful to study cases with \(D_i>0.5\) and essential to review cases with \(D_i>1\).

fit1<-lm(NO3~mxPH+Cl,data=algae) #
cooksd <- cooks.distance(fit1) #store the D of each point in a named numeric form
plot(cooksd, pch="*", cex=2, main="Influential Obs by Cooks distance")  # plot cook's distance
abline(h = 4*mean(cooksd, na.rm=T), col="red")  # add cutoff line
text(x=1:length(cooksd)+1, y=cooksd, col="red",adj = c(0.15,0.1), 
     labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""))  # add labels

Get the information about these extreme observations.

influential=cooksd[cooksd>(4*mean(cooksd, na.rm=T))]  # get the influential row numbers
filtered_data<-algae[-which(row.names(algae) %in% names(influential)),] #remove the rows inside influential
head(algae[names(influential), ])  # influential observations.
##     season   size  speed mxPH mnO2      Cl    NO3   NH4 oPO4 PO4 Chla  a1  a2
## 134 spring medium medium  7.9  8.3 391.500  6.045   380  173 317  5.5 2.4 1.7
## 153 autumn medium   high  7.3 11.8  44.205 45.650 24064   44  34 53.1 2.2 0.0
##      a3  a4  a5   a6  a7
## 134 4.2 8.3 1.7  0.0 2.4
## 153 0.0 1.2 5.9 77.6 0.0

We do not have a variable unique to each observation, so we will need to use row.names() as well as the names(), which takes the row number for each value inside influential.

We can use a different package for this style of plot.

a=olsrr::ols_plot_cooksd_chart(model = fit1,print_plot = T)+
    theme_bw()

In this case we got values that are different to the rows in the main data set so we can have a look at the plot data to get the actual rows.

values=row.names(a$data[(a$data$fct_color=="outlier"),])
values
## [1] "134" "153"
fit2<-lm(NO3~mxPH+Cl,data=algae[-as.integer(values),]) #our second regression

Test for outliers

The function outlierTest from the car package gives the most extreme observation based on the given model using the Bonferroni Outlier Test.

car::outlierTest(fit1,n.max=Inf) #the largest influence cause by this value 
##     rstudent unadjusted p-value Bonferroni p
## 153 19.80916         1.3401e-47   2.5328e-45
fit3<-lm(NO3~mxPH+Cl,data=algae[-153,]) #remove that observation
car::outlierTest(fit3) #the largest influence in the new model
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 5 3.490858         0.00060268       0.1133

With the package Outliers, we calculate the extreme observations. First, let’s look at the extremest value with outlier.

outlier(algae$Cl)
## [1] 391.5

Use the scores function to calculates standard, t, chi-squared, IQR, and MAD scores of given data.

head(scores(na.omit(algae$Cl), type="chisq"))
## [1] 0.134322963 0.090826050 0.005962813 0.518682218 0.062562936 0.222972355

We can test which observations are beyond a given probability (percentile).

head(scores(na.omit(algae$Cl), type="chisq",prob = 0.95),40)
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [37] FALSE FALSE FALSE FALSE