README.Rmd

---
output: 
  github_document:
    toc: false
    fig_width: 10.08
    fig_height: 6
tags: [r, reports]
vignette: >
  \usepackage[utf8]{inputenc}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
bibliography: paper/paper.bib
---

# performance <img src='man/figures/logo.png' align="right" height="139" />

```{r, echo = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  tidy.opts = list(width.cutoff = 60),
  tidy = TRUE,
  fig.path = "README-",
  comment = "#>"
)
options(knitr.kable.NA = '',
        digits = 4,
        width = 60)
```


[![CRAN](http://www.r-pkg.org/badges/version/performance)](https://cran.r-project.org/package=performance) [![Documentation](https://img.shields.io/badge/documentation-performance-orange.svg?colorB=E91E63)](https://easystats.github.io/performance/) [![downloads](http://cranlogs.r-pkg.org/badges/performance)](https://cran.r-project.org/package=performance) [![total](http://cranlogs.r-pkg.org/badges/grand-total/performance)](http://cranlogs.r-pkg.org/)

***Test if your model is a good model!***

The primary goal of the **performance** package is to provide utilities for computing **indices of model quality** and **goodness of fit**. This includes 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

Run the following:

```{r eval=FALSE, message=FALSE, warning=FALSE}
install.packages("devtools")
devtools::install_github("easystats/performance")
```
```{r message=FALSE, warning=FALSE}
library("performance")
```


# Examples

[![Features](https://img.shields.io/badge/features-performance-orange.svg?colorB=2196F3)](https://easystats.github.io/performance/reference/index.html) [![Blog](https://img.shields.io/badge/blog-easystats-orange.svg?colorB=FF9800)](https://easystats.github.io/blog/posts/)

## 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.

```{r, warning=FALSE, message=FALSE}
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](https://easystats.github.io/performance/reference/index.html)).

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].

```{r, warning=FALSE, message=FALSE}
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.

```{r, warning=FALSE, message=FALSE}
library(lme4)
model <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
icc(model)
```

For models of class `brmsfit`, an ICC based on variance decomposition is returned (for details, see the [documentation](https://easystats.github.io/performance/reference/icc.html)).
 
```{r, warning=FALSE, message=FALSE, echo=FALSE}
model <- insight::download_model("brms_mixed_1")
```
```{r, warning=FALSE, message=FALSE, eval=FALSE}
library(brms)
set.seed(123)
model <- brm(mpg ~ wt + (1 | cyl) + (1 + wt | gear), data = mtcars)
```
```{r, warning=FALSE, message=FALSE}
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.

```{r, warning=FALSE, message=TRUE}
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.

```{r, warning=FALSE, message=TRUE}
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.

```{r, warning=FALSE, message=FALSE}
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](https://easystats.github.io/performance/reference/check_singularity.html). 

## 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

```{r, warning=FALSE, message=FALSE, eval=FALSE}
m1 <- lm(mpg ~ wt + cyl, data = mtcars)
model_performance(m1)
```
```{r, warning=FALSE, message=FALSE, echo=FALSE}
m1 <- lm(mpg ~ wt + cyl, data = mtcars)
knitr::kable(model_performance(m1), digits = 2)
```

### Logistic regression

```{r, warning=FALSE, message=FALSE, eval=FALSE}
m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")
model_performance(m2)
```
```{r, warning=FALSE, message=FALSE, echo=FALSE}
m2 <- glm(vs ~ wt + mpg, data = mtcars, family = "binomial")
knitr::kable(model_performance(m2), digits = 2)
```

### Linear mixed model

```{r, warning=FALSE, message=FALSE, eval=FALSE}
library(lme4)
m3 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
model_performance(m3)
```
```{r, warning=FALSE, message=FALSE, echo=FALSE}
library(lme4)
m3 <- lmer(Reaction ~ Days + (1 + Days | Subject), data = sleepstudy)
knitr::kable(model_performance(m3), digits = 2)
```

### Comparing different models

```{r, warning=FALSE, message=FALSE, eval=FALSE}
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)
```
```{r, warning=FALSE, message=FALSE, echo=FALSE}
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())

knitr::kable(compare_performance(m1, m2, m3, m4), digits = 2)
```

# References