output: github_document
bibliography: paper.bib

performance

knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  out.width = "100%",
  dpi = 150,
  fig.path = "man/figures/",
  comment = "#>"
)

options(
  knitr.kable.NA = "",
  digits = 4,
  width = 100
)

library(performance)

Test if your model is a good model!

A crucial aspect when building regression models is to evaluate the quality of modelfit. It is important to investigate how well models fit to the data and which fit indices to report. Functions to create diagnostic plots or to compute fit measures do exist, however, mostly spread over different packages. There is no unique and consistent approach to assess the model quality for different kind of models.

The primary goal of the performance package is to fill this gap and to provide utilities for computing indices of model quality and goodness of fit. These include measures like r-squared (R2), root mean squared error (RMSE) or intraclass correlation coefficient (ICC) , but also functions to check (mixed) models for overdispersion, zero-inflation, convergence or singularity.

Installation

The performance package is available on CRAN, while its latest development version is available on R-universe (from rOpenSci).

Type	Source	Command
Release	CRAN	`install.packages("performance")`
Development	R-universe	`install.packages("performance", repos = "https://easystats.r-universe.dev")`

Once you have downloaded the package, you can then load it using:

library("performance")

Tip

Instead of library(performance), use library(easystats). This will make all features of the easystats-ecosystem available.

To stay updated, use easystats::install_latest().

Citation

To cite performance in publications use:

citation("performance")

Documentation

There is a nice introduction into the package on youtube.

The performance workflow

knitr::include_graphics("man/figures/figure_workflow.png")

Assessing model quality

R-squared

performance has a generic r2() function, which computes the r-squared for many different models, including mixed effects and Bayesian regression models.

r2() returns a list containing values related to the "most appropriate" r-squared for the given model.

model <- lm(mpg ~ wt + cyl, data = mtcars)
r2(model)

model <- glm(am ~ wt + cyl, data = mtcars, family = binomial)
r2(model)

library(MASS)
data(housing)
model <- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
r2(model)

The different R-squared measures can also be accessed directly via functions like r2_bayes(), r2_coxsnell() or r2_nagelkerke() (see a full list of functions here).

For mixed models, the conditional and marginal R-squared are returned. The marginal R-squared considers only the variance of the fixed effects and indicates how much of the model's variance is explained by the fixed effects part only. The conditional R-squared takes both the fixed and random effects into account and indicates how much of the model's variance is explained by the "complete" model.

For frequentist mixed models, r2() (resp. r2_nakagawa()) computes the mean random effect variances, thus r2() is also appropriate for mixed models with more complex random effects structures, like random slopes or nested random effects [@johnson_extension_2014; @nakagawa_coefficient_2017].

set.seed(123)
library(rstanarm)

model <- stan_glmer(
  Petal.Length ~ Petal.Width + (1 | Species),
  data = iris,
  cores = 4
)

r2(model)

library(lme4)
model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
r2(model)

Intraclass Correlation Coefficient (ICC)

Similar to R-squared, the ICC provides information on the explained variance and can be interpreted as "the proportion of the variance explained by the grouping structure in the population" [@hox_multilevel_2010].

icc() calculates the ICC for various mixed model objects, including stanreg models.

library(lme4)
model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
icc(model)

...and models of class brmsfit.

model <- insight::download_model("brms_mixed_1")

library(brms)
set.seed(123)
model <- brm(mpg ~ wt + (1 | cyl) + (1 + wt | gear), data = mtcars)

icc(model)

Model diagnostics

Check for overdispersion

Overdispersion occurs when the observed variance in the data is higher than the expected variance from the model assumption (for Poisson, variance roughly equals the mean of an outcome). check_overdispersion() checks if a count model (including mixed models) is overdispersed or not.

library(glmmTMB)
data(Salamanders)
model <- glm(count ~ spp + mined, family = poisson, data = Salamanders)
check_overdispersion(model)

Overdispersion can be fixed by either modelling the dispersion parameter (not possible with all packages), or by choosing a different distributional family (like Quasi-Poisson, or negative binomial, see [@gelman_data_2007]).

Check for zero-inflation

Zero-inflation (in (Quasi-)Poisson models) is indicated when the amount of observed zeros is larger than the amount of predicted zeros, so the model is underfitting zeros. In such cases, it is recommended to use negative binomial or zero-inflated models.

Use check_zeroinflation() to check if zero-inflation is present in the fitted model.

model <- glm(count ~ spp + mined, family = poisson, data = Salamanders)
check_zeroinflation(model)

Check for singular model fits

A "singular" model fit means that some dimensions of the variance-covariance matrix have been estimated as exactly zero. This often occurs for mixed models with overly complex random effects structures.

check_singularity() checks mixed models (of class lme, merMod, glmmTMB or MixMod) for singularity, and returns TRUE if the model fit is singular.

library(lme4)
data(sleepstudy)

# prepare data
set.seed(123)
sleepstudy$mygrp <- sample(1:5, size = 180, replace = TRUE)
sleepstudy$mysubgrp <- NA
for (i in 1:5) {
  filter_group <- sleepstudy$mygrp == i
  sleepstudy$mysubgrp[filter_group] <-
    sample(1:30, size = sum(filter_group), replace = TRUE)
}

# fit strange model
model <- lmer(
  Reaction ~ Days + (1 | mygrp / mysubgrp) + (1 | Subject),
  data = sleepstudy
)

check_singularity(model)

Remedies to cure issues with singular fits can be found here.

Check for heteroskedasticity

Linear models assume constant error variance (homoskedasticity).

The check_heteroscedasticity() functions assess if this assumption has been violated:

data(cars)
model <- lm(dist ~ speed, data = cars)

check_heteroscedasticity(model)

Comprehensive visualization of model checks

performance provides many functions to check model assumptions, like check_collinearity(), check_normality() or check_heteroscedasticity(). To get a comprehensive check, use check_model().

# defining a model
model <- lm(mpg ~ wt + am + gear + vs * cyl, data = mtcars)

# checking model assumptions
check_model(model)

Model performance summaries

model_performance() computes indices of model performance for regression models. Depending on the model object, typical indices might be r-squared, AIC, BIC, RMSE, ICC or LOOIC.

Linear model

m1 <- lm(mpg ~ wt + cyl, data = mtcars)
model_performance(m1)

Logistic regression

m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")
model_performance(m2)

Linear mixed model

library(lme4)
m3 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
model_performance(m3)

Models comparison

The compare_performance() function can be used to compare the performance and quality of several models (including models of different types).

counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
m4 <- glm(counts ~ outcome + treatment, family = poisson())

compare_performance(m1, m2, m3, m4, verbose = FALSE)

General index of model performance

One can also easily compute and a composite index of model performance and sort the models from the best one to the worse.

compare_performance(m1, m2, m3, m4, rank = TRUE, verbose = FALSE)

Visualisation of indices of models' performance

Finally, we provide convenient visualisation (the see package must be installed).

plot(compare_performance(m1, m2, m4, rank = TRUE, verbose = FALSE))

Testing models

test_performance() (and test_bf, its Bayesian sister) carries out the most relevant and appropriate tests based on the input (for instance, whether the models are nested or not).

set.seed(123)
data(iris)

lm1 <- lm(Sepal.Length ~ Species, data = iris)
lm2 <- lm(Sepal.Length ~ Species + Petal.Length, data = iris)
lm3 <- lm(Sepal.Length ~ Species * Sepal.Width, data = iris)
lm4 <- lm(Sepal.Length ~ Species * Sepal.Width + Petal.Length + Petal.Width, data = iris)

test_performance(lm1, lm2, lm3, lm4)

test_bf(lm1, lm2, lm3, lm4)

Plotting Functions

Plotting functions are available through the see package.

Code of Conduct

Please note that the performance project is released with a Contributor Code of Conduct. By contributing to this project, you agree to abide by its terms.