This vignette shows how to use the *performance*
package to check the fit of a model, how to detect misspecification
and how to improve your model. The basic workflow of the
*performance* package can be summarized as follows:

- fit a regression model
- check the model fit and assess model fit indices
- if necessary, fit another model that could potentially improve the fit
- compare the model fit indices and perform statistical tests to determine which model is the best fit

In the following, we will demonstrate this workflow using a model
with a count response variable. We will fit a Poisson regression model
to the Salamanders dataset from the *glmmTMB* package. The
dataset contains counts of salamanders in different sites, along with
information on the number of mines and the species of salamanders. We
will check the model fit and assess the model fit indices.

Problems that may arise with count response variables are *zero
inflation* and *overdispersion*. Zero inflation occurs when
there are more zeros in the data than expected under the Poisson
distribution. Overdispersion occurs when the variance of the data is
greater than the mean, which violates the assumption of equidispersion
in the Poisson distribution.

We will check for these problems and suggest ways to improve the model fit, i.e.Â if necessary, we will fit another model that could potentially improve the fit. Finally, we will compare the model fit indices and perform statistical tests to determine which model is the best fit.

We start with a generalized mixed effects model, using a Poisson distribution.

```
library(performance)
model1 <- glmmTMB::glmmTMB(
count ~ mined + spp + (1 | site),
family = poisson,
data = glmmTMB::Salamanders
)
```

First, let us look at the summary of the model.

```
library(parameters)
model_parameters(model1)
#> # Fixed Effects
#>
#> Parameter | Log-Mean | SE | 95% CI | z | p
#> ---------------------------------------------------------------
#> (Intercept) | -1.62 | 0.24 | [-2.10, -1.15] | -6.76 | < .001
#> mined [no] | 2.26 | 0.28 | [ 1.72, 2.81] | 8.08 | < .001
#> spp [PR] | -1.39 | 0.22 | [-1.81, -0.96] | -6.44 | < .001
#> spp [DM] | 0.23 | 0.13 | [-0.02, 0.48] | 1.79 | 0.074
#> spp [EC-A] | -0.77 | 0.17 | [-1.11, -0.43] | -4.50 | < .001
#> spp [EC-L] | 0.62 | 0.12 | [ 0.39, 0.86] | 5.21 | < .001
#> spp [DES-L] | 0.68 | 0.12 | [ 0.45, 0.91] | 5.75 | < .001
#> spp [DF] | 0.08 | 0.13 | [-0.18, 0.34] | 0.60 | 0.549
#>
#> # Random Effects
#>
#> Parameter | Coefficient | 95% CI
#> -------------------------------------------------
#> SD (Intercept: site) | 0.58 | [0.38, 0.87]
#>
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#> using a Wald z-distribution approximation.
#>
#> The model has a log- or logit-link. Consider using `exponentiate =
#> TRUE` to interpret coefficients as ratios.
```

We see a lot of statistically significant estimates here. No matter, which philosophy you follow in terms of interpreting statistical test results, our conclusions we draw from our regression models will be inaccurate if our modeling assumptions are a poor fit for the situation. Hence, checking model fit is essential.

In *performance*, we can conduct a comprehensive visual
inspection of our model fit using `check_model()`

. We wonâ€™t
go into details of all the plots here, but you can find more information
on all created diagnostic plots in the dedicated
vignette.

For now, we want to focus on the *posterior predictive
checks*, *dispersion and zero-inflation* as well as the Q-Q
plot (*uniformity of residuals*).

Note that unlike `plot()`

, which is a base R function to
create diagnostic plots, `check_model()`

relies on
*simulated residuals* for the Q-Q plot, which is more accurate
for non-Gaussian models. See this
vignette and the documentation of `simulate_residuals()`

for further details.

The above plot suggests that we may have issues with overdispersion
and/or zero-inflation. We can check for these problems using
`check_overdispersion()`

and
`check_zeroinflation()`

, which will perform statistical tests
(based on simulated residuals). These tests can additionally be used
beyond the visual inspection.

```
check_overdispersion(model1)
#> # Overdispersion test
#>
#> dispersion ratio = 2.324
#> Pearson's Chi-Squared = 1475.875
#> p-value = < 0.001
#> Overdispersion detected.
check_zeroinflation(model1)
#> # Check for zero-inflation
#>
#> Observed zeros: 387
#> Predicted zeros: 311
#> Ratio: 0.80
#> Model is underfitting zeros (probable zero-inflation).
```

As we can see, our model seems to suffer both from overdispersion and zero-inflation.

We can try to improve the model fit by fitting a model with zero-inflation component:

```
model2 <- glmmTMB::glmmTMB(
count ~ mined + spp + (1 | site),
ziformula = ~ mined + spp,
family = poisson,
data = glmmTMB::Salamanders
)
check_model(model2)
#> `check_outliers()` does not yet support models of class `glmmTMB`.
```

Looking at the above plots, the zero-inflation seems to be addressed
properly (see especially *posterior predictive checks* and
*uniformity of residuals*, the Q-Q plot). However, the
overdispersion still could be present. We can check for these problems
using `check_overdispersion()`

and
`check_zeroinflation()`

again.

```
check_overdispersion(model2)
#> # Overdispersion test
#>
#> dispersion ratio = 1.679
#> p-value = 0.008
#> Overdispersion detected.
check_zeroinflation(model2)
#> # Check for zero-inflation
#>
#> Observed zeros: 387
#> Predicted zeros: 387
#> Ratio: 1.00
#> Model seems ok, ratio of observed and predicted zeros is within the
#> tolerance range (p > .999).
```

Indeed, the overdispersion is still present.

We can try to address this issue by fitting a negative binomial model instead of using a Poisson distribution.