Quickly test linear regression assumptions

Getting Started

Load the rempsyc package:

library(rempsyc)

Note: If you haven’t installed this package yet, you will need to install it via the following command: install.packages("rempsyc"). Furthermore, you may be asked to install the following packages if you haven’t installed them already (you may decide to install them all now to avoid interrupting your workflow if you wish to follow this tutorial from beginning to end):

pkgs <- c(
  "flextable", "performance", "see", "lmtest", "ggplot2",
  "qqplotr", "ggrepel", "patchwork", "boot"
)
install_if_not_installed(pkgs)

The raw output doesn’t look so nice in the console because the column names are so long, but it does look good in the viewer or when exported to a word processing software because the column names get wrapped. You can try it in RStudio yourself with the following command:

View(nice_assumptions(model))

nice_table(nice_assumptions(model), col.format.p = 2:4)

Interpretation: (p) values < .05 imply assumptions are not respected. Diagnostic is how many assumptions are not respected for a given model or variable.

This function is particularly useful when testing several models simultaneously as it makes for a nice table of results. Let’s make a short demonstration of this.

# Define our dependent variables
DV <- names(mtcars[-1])
# Make list of all formulas
formulas <- paste(DV, "~ mpg")
# Make list of all models
models.list <- lapply(X = formulas, FUN = lm, data = mtcars)
# Make diagnostic table
assumptions.table <- nice_assumptions(models.list)

Use the Viewer for better results

View(assumptions.table)

Or alternatively with the nice_table() function

nice_table(assumptions.table, col.format.p = 2:4)

Make the basic plot

nice_qq(
  data = iris,
  variable = "Sepal.Length",
  group = "Species"
)

Generally speaking, as long as the dots lie within the confidence band, all is good and the data distribute normally.

Further customization

Change (or reorder) colours, x-axis title or y-axis title, names of groups, no grid lines, or add the p-value from the Shapiro-Wilk (if you wish so).

nice_qq(
  data = iris,
  variable = "Sepal.Length",
  group = "Species",
  colours = c("#00BA38", "#619CFF", "#F8766D"),
  groups.labels = c("(a) Setosa", "(b) Versicolor", "(c) Virginica"),
  grid = FALSE,
  shapiro = TRUE,
  title = NULL
)

A p > .05 suggests that your data is normally distributed (according to the Shapiro-Wilk test). Inversely, a p < .05 suggests your data is not normally distributed. Keep in mind however that there are several known problems with objective tests, so visual assessment is recommended. Nonetheless, for beginners, it can sometimes be useful (or rather, reassuring) to be able to pair the visual assessment to the tests values.

Make the basic plot

nice_density(
  data = iris,
  variable = "Sepal.Length",
  group = "Species"
)

Pro tip: Save density plots with a narrower width for better-looking results.

Further customization

Change (or reorder) colours, x-axis title or y-axis title, names of groups, no grid lines, or add the p-value from the Shapiro-Wilk (if you wish so).

nice_density(
  data = iris,
  variable = "Sepal.Length",
  group = "Species",
  colours = c("#00BA38", "#619CFF", "#F8766D"),
  xtitle = "Sepal Length",
  ytitle = "Density (vs. Normal Distribution)",
  groups.labels = c("(a) Setosa", "(b) Versicolor", "(c) Virginica"),
  grid = FALSE,
  shapiro = TRUE,
  histogram = TRUE,
  title = "Density (Sepal Length)"
)

Make the basic plot

plot_outliers(
  airquality,
  group = "Month",
  response = "Ozone"
)

Be aware that checking outliers by group and not by group will yield different results, e.g.,:

plot_outliers(
  airquality,
  response = "Ozone"
)

Further customization

Change (or reorder) colours, y-axis title, names of groups, or detection method.

plot_outliers(
  airquality,
  group = "Month",
  response = "Ozone",
  method = "sd",
  criteria = 3.29,
  colours = c("white", "black", "purple", "grey", "pink"),
  ytitle = "Ozone",
  xtitle = "Month of the Year"
)

Flagging outliers

Outliers can also be identifed with find_mad and multiple variables:

find_mad(airquality, names(airquality), criteria = 3)

## 8 outlier(s) based on 3 median absolute deviations for variable(s): 
##  Ozone, Solar.R, Wind, Temp, Month, Day 
## 
## Outliers per variable: 
## 
## $Ozone
##   Row Ozone_mad
## 1  30  3.218284
## 2  62  3.989131
## 3  99  3.488081
## 4 101  3.025573
## 5 117  5.261028
## 6 121  3.333911
## 
## $Wind
##   Row Wind_mad
## 1   9 3.049871
## 2  48 3.225825

Winsorizing outliers

Outliers can also be winsorized with winsorize_mad, meaning that outlier values outside 3 MADs will be brought back to 3 MADs:

winsorize_mad(airquality$Ozone, criteria = 3) |>
  head(30)

##  [1]  41.0000  36.0000  12.0000  18.0000       NA  28.0000  23.0000  19.0000
##  [9]   8.0000       NA   7.0000  16.0000  11.0000  14.0000  18.0000  14.0000
## [17]  34.0000   6.0000  30.0000  11.0000   1.0000  11.0000   4.0000  32.0000
## [25]       NA       NA       NA  23.0000  45.0000 109.3365

Multivariate outliers

For multivariate outliers, it is recommended to use the Minimum Covariance Determinant (MCD), as available in the performance package from easystats mentioned in the beginning of this post.

check_outliers(na.omit(airquality), method = "mcd")

## 7 outliers detected: cases 6, 7, 11, 14, 23, 34, 77.
## - Based on the following method and threshold: mcd (22.458).
## - For variables: Ozone, Solar.R, Wind, Temp, Month, Day.

We note also that ideally the handling of outliers should be preregistered a priori. For more up-to-date recommendations on handling outliers, please see:

Leys, C., Delacre, M., Mora, Y. L., Lakens, D., & Ley, C. (2019). How to classify, detect, and manage univariate and multivariate outliers, with emphasis on pre-registration. International Review of Social Psychology, 32(1). https://doi.org/10.5334/irsp.289

Make the basic table

View(nice_var(
  data = iris,
  variable = "Sepal.Length",
  group = "Species"
))

Let’s now try it for many variables to see how handy it can be.

# Define our dependent variables
DV <- names(iris[1:4])

# Make diagnostic table
var.table <- nice_var(
  data = iris,
  variable = DV,
  group = "Species"
)

Use the Viewer for better results (or export it to Word using my nice_table() function)

View(var.table)

Visualizing Variance

It can also be interesting to look at variance visually for learning purposes.

Make the basic plot

nice_varplot(
  data = iris,
  variable = "Sepal.Length",
  group = "Species"
)

Further customization

Change (or reorder) colours, y-axis title, or names of groups.

nice_varplot(
  data = iris,
  variable = "Sepal.Length",
  group = "Species",
  colours = c("#00BA38", "#619CFF", "#F8766D"),
  ytitle = "Sepal Length",
  groups.labels = c("(a) Setosa", "(b) Versicolor", "(c) Virginica")
)

Thanks for checking in

Make sure to check out this page again if you use the code after a time or if you encounter errors, as I periodically update or improve the code. Feel free to contact me for comments, questions, or requests to improve this function at https://github.com/rempsyc/rempsyc/issues. See all tutorials here: https://remi-theriault.com/tutorials.