diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index afc762333..22db08405 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 0.10.9 -Date: 2024-02-17 07:56:08 UTC -SHA: 051016febd197937ad083266b630c871fa9e1623 +Version: 0.11.0 +Date: 2024-03-22 21:30:58 UTC +SHA: 051b9bb2b7721c632ce145f85c55aa55c8eebf90 diff --git a/DESCRIPTION b/DESCRIPTION index 0f23a585c..68ce7de42 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: performance Title: Assessment of Regression Models Performance -Version: 0.10.9.10 +Version: 0.11.0 Authors@R: c(person(given = "Daniel", family = "Lüdecke", @@ -154,4 +154,3 @@ Config/Needs/website: r-lib/pkgdown, easystats/easystatstemplate Config/rcmdcheck/ignore-inconsequential-notes: true -Remotes: easystats/see diff --git a/NAMESPACE b/NAMESPACE index 24daed512..2ecdbf751 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -72,6 +72,7 @@ S3method(check_normality,htest) S3method(check_normality,lmerModLmerTest) S3method(check_normality,merMod) S3method(check_normality,numeric) +S3method(check_normality,performance_simres) S3method(check_outliers,BFBayesFactor) S3method(check_outliers,DHARMa) S3method(check_outliers,character) @@ -517,6 +518,7 @@ S3method(r2_tjur,nestedLogit) S3method(residuals,BFBayesFactor) S3method(residuals,check_normality_numeric) S3method(residuals,iv_robust) +S3method(residuals,performance_simres) S3method(rstudent,check_normality_numeric) S3method(test_bf,ListModels) S3method(test_bf,default) diff --git a/NEWS.md b/NEWS.md index 6fb76ca53..b12a34e96 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# performance 0.10.10 +# performance 0.11.0 ## New supported models diff --git a/R/check_model.R b/R/check_model.R index f75a103a8..cf33f3eef 100644 --- a/R/check_model.R +++ b/R/check_model.R @@ -144,12 +144,13 @@ #' inside the error bounds. See [`binned_residuals()`] for further details. #' #' @section Residuals for (Generalized) Linear Models: -#' Plots that check the normality of residuals (Q-Q plot) or the homogeneity of -#' variance use standardized Pearson's residuals for generalized linear models, -#' and standardized residuals for linear models. The plots for the normality of -#' residuals (with overlayed normal curve) and for the linearity assumption use -#' the default residuals for `lm` and `glm` (which are deviance residuals for -#' `glm`). +#' Plots that check the homogeneity of variance use standardized Pearson's +#' residuals for generalized linear models, and standardized residuals for +#' linear models. The plots for the normality of residuals (with overlayed +#' normal curve) and for the linearity assumption use the default residuals +#' for `lm` and `glm` (which are deviance residuals for `glm`). The Q-Q plots +#' use simulated residuals (see [`simulate_residuals()`]) for non-Gaussian +#' models and standardized residuals for linear models. #' #' @section Troubleshooting: #' For models with many observations, or for more complex models in general, diff --git a/R/check_model_diagnostics.R b/R/check_model_diagnostics.R index 39a3c843d..431a2bc0f 100644 --- a/R/check_model_diagnostics.R +++ b/R/check_model_diagnostics.R @@ -301,7 +301,7 @@ d <- data.frame(Predicted = predicted) # residuals based on simulated residuals - but we want normally distributed residuals - d$Residuals <- stats::residuals(simres, quantileFunction = stats::qnorm, ...) + d$Residuals <- stats::residuals(simres, quantile_function = stats::qnorm, ...) d$Res2 <- d$Residuals^2 d$StdRes <- insight::get_residuals(model, type = "pearson") diff --git a/R/check_normality.R b/R/check_normality.R index 2b9f071d5..3c2902915 100644 --- a/R/check_normality.R +++ b/R/check_normality.R @@ -92,6 +92,23 @@ check_normality.glm <- function(x, ...) { invisible(out) } +# simulated residuals ---------- + +#' @export +check_normality.performance_simres <- function(x, ...) { + # check for normality of residuals + res <- stats::residuals(x, quantile_function = stats::qnorm) + p.val <- .check_normality(res[!is.infinite(res) & !is.na(res)], x) + + attr(p.val, "data") <- x + attr(p.val, "object_name") <- insight::safe_deparse_symbol(substitute(x)) + attr(p.val, "effects") <- "fixed" + class(p.val) <- unique(c("check_normality", "see_check_normality", class(p.val))) + + p.val +} + + # numeric ------------------- #' @export diff --git a/R/check_overdispersion.R b/R/check_overdispersion.R index 0a9cfb595..fa46edfef 100644 --- a/R/check_overdispersion.R +++ b/R/check_overdispersion.R @@ -62,19 +62,10 @@ #' multilevel/hierarchical models. Cambridge; New York: Cambridge University #' Press. #' -#' @examplesIf getRversion() >= "4.0.0" && require("glmmTMB", quietly = TRUE) -#' -#' library(glmmTMB) -#' data(Salamanders) +#' @examplesIf getRversion() >= "4.0.0" && require("glmmTMB") +#' data(Salamanders, package = "glmmTMB") #' m <- glm(count ~ spp + mined, family = poisson, data = Salamanders) #' check_overdispersion(m) -#' -#' m <- glmmTMB( -#' count ~ mined + spp + (1 | site), -#' family = poisson, -#' data = Salamanders -#' ) -#' check_overdispersion(m) #' @export check_overdispersion <- function(x, ...) { UseMethod("check_overdispersion") diff --git a/R/check_residuals.R b/R/check_residuals.R index 2bb284c88..d8840e165 100644 --- a/R/check_residuals.R +++ b/R/check_residuals.R @@ -19,7 +19,8 @@ #' #' @inheritSection simulate_residuals Tests based on simulated residuals #' -#' @seealso [`simulate_residuals()`] and [`check_predictions()`]. +#' @seealso [`simulate_residuals()`], [`check_zeroinflation()`], +#' [`check_overdispersion()`] and [`check_predictions()`]. #' #' @return The p-value of the test statistics. #' diff --git a/R/check_zeroinflation.R b/R/check_zeroinflation.R index aaa5e7641..5f87941f3 100644 --- a/R/check_zeroinflation.R +++ b/R/check_zeroinflation.R @@ -30,7 +30,7 @@ #' #' @section Tests based on simulated residuals: #' For certain models, resp. model from certain families, tests are based on -#' simulated residuals (see [`simulated_residual()`]). These are usually more +#' simulated residuals (see [`simulate_residuals()`]). These are usually more #' accurate for testing such models than the traditionally used Pearson residuals. #' However, when simulating from more complex models, such as mixed models or #' models with zero-inflation, there are several important considerations. diff --git a/R/simulate_residuals.R b/R/simulate_residuals.R index de9ee0f4d..41fccea55 100644 --- a/R/simulate_residuals.R +++ b/R/simulate_residuals.R @@ -8,17 +8,25 @@ #' @param x A model object. #' @param iterations Number of simulations to run. #' @param ... Arguments passed on to [`DHARMa::simulateResiduals()`]. +#' @param object A `performance_simres` object, as returned by `simulate_residuals()`. +#' @param quantile_function A function to apply to the residuals. If `NULL`, the +#' residuals are returned as is. If not `NULL`, the residuals are passed to this +#' function. This is useful for returning normally distributed residuals, for +#' example: `residuals(x, quantile_function = qnorm)`. +#' @param outlier_values A vector of length 2, specifying the values to replace +#' `-Inf` and `Inf` with, respectively. #' #' @return Simulated residuals, which can be further processed with #' [`check_residuals()`]. The returned object is of class `DHARMa` and #' `performance_simres`. #' -#' @seealso [`check_residuals()`] and [`check_predictions()`]. +#' @seealso [`check_residuals()`], [`check_zeroinflation()`], +#' [`check_overdispersion()`] and [`check_predictions()`]. #' #' @details This function is a small wrapper around [`DHARMa::simulateResiduals()`]. #' It basically only sets `plot = FALSE` and adds an additional class attribute #' (`"performance_sim_res"`), which allows using the DHARMa object in own plotting -#' functions in the **see** package. See also `vignette("DHARMa")`. There is a +#' functions from the **see** package. See also `vignette("DHARMa")`. There is a #' `plot()` method to visualize the distribution of the residuals. #' #' @section Tests based on simulated residuals: @@ -50,6 +58,9 @@ #' m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) #' simulate_residuals(m) #' +#' # extract residuals +#' head(residuals(simulate_residuals(m))) +#' #' @export simulate_residuals <- function(x, iterations = 250, ...) { insight::check_if_installed("DHARMa") @@ -78,9 +89,10 @@ print.performance_simres <- function(x, ...) { # DHARMa's method. msg <- paste0( "Simulated residuals from a model of class `", class(x$fittedModel)[1], - "` based on ", x$nSim, " simulations. Use `check_residuals()` to check ", - "uniformity of residuals. It is recommended to refer to `?DHARMa::simulateResiudals`", - " and `vignette(\"DHARMa\")` for more information about different settings", + "` based on ", x$nSim, " simulations. Use `check_residuals()` to check", + " uniformity of residuals or `residuals()` to extract simulated residuals.", + " It is recommended to refer to `?DHARMa::simulateResiudals` and", + " `vignette(\"DHARMa\")` for more information about different settings", " in particular situations or for particular models.\n" ) cat(insight::format_message(msg)) @@ -93,6 +105,37 @@ plot.performance_simres <- function(x, ...) { } +# methods -------------------------- + +#' @rdname simulate_residuals +#' @export +residuals.performance_simres <- function(object, quantile_function = NULL, outlier_values = NULL, ...) { + # check for DHARMa argument names + dots <- list(...) + if (!is.null(dots$quantileFunction)) { + quantile_function <- dots$quantileFunction + } + if (!is.null(dots$outlierValues)) { + outlier_values <- dots$outlierValues + } + + if (is.null(quantile_function)) { + res <- object$scaledResiduals + } else { + res <- quantile_function(object$scaledResiduals) + if (!is.null(outlier_values)) { + # check for correct length of outlier_values + if (length(outlier_values) != 2) { + insight::format_error("`outlier_values` must be a vector of length 2.") + } + res[res == -Inf] <- outlier_values[1] + res[res == Inf] <- outlier_values[2] + } + } + res +} + + # helper functions --------------------- .simres_statistics <- function(x, statistic_fun, alternative = "two.sided") { diff --git a/WIP/check_model_logistic.Rmd b/WIP/check_model_logistic.Rmd new file mode 100644 index 000000000..e6b7cafe8 --- /dev/null +++ b/WIP/check_model_logistic.Rmd @@ -0,0 +1,204 @@ +--- +title: "Checking model assumption - logistic regression models" +output: + rmarkdown::html_vignette: + toc: true + fig_width: 10.08 + fig_height: 6 +tags: [r, performance, r2] +vignette: > + \usepackage[utf8]{inputenc} + %\VignetteIndexEntry{Checking model assumption - logistic regression models} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +```{r , include=FALSE} +library(knitr) +library(performance) +options(knitr.kable.NA = "") +knitr::opts_chunk$set( + comment = ">", + message = FALSE, + warning = FALSE, + out.width = "100%", + dpi = 450 +) +options(digits = 2) + +pkgs <- c("see", "ggplot2", "datawizard", "parameters") +successfully_loaded <- vapply(pkgs, requireNamespace, FUN.VALUE = logical(1L), quietly = TRUE) +can_evaluate <- all(successfully_loaded) + +if (can_evaluate) { + knitr::opts_chunk$set(eval = TRUE) + vapply(pkgs, require, FUN.VALUE = logical(1L), quietly = TRUE, character.only = TRUE) +} else { + knitr::opts_chunk$set(eval = FALSE) +} +``` + +# Make sure your model inference is accurate! + +Model diagnostics is crucial, because parameter estimation, p-values and confidence interval depend on correct model assumptions as well as on the data. If model assumptions are violated, estimates can be statistically significant "even if the effect under study is null" (_Gelman/Greenland 2019_). + +There are several problems associated with model diagnostics. Different types of models require different checks. For instance, normally distributed residuals are assumed to apply for linear regression, but is no appropriate assumption for logistic regression. Furthermore, it is recommended to carry out visual inspections, i.e. to generate and inspect so called diagnostic plots of model assumptions - formal statistical tests are often too strict and warn of violation of the model assumptions, although everything is fine within a certain tolerance range. But how should such diagnostic plots be interpreted? And if violations have been detected, how to fix them? + +This vignette introduces the `check_model()` function of the **performance** package, shows how to use this function for logistic regression models and how the resulting diagnostic plots should be interpreted. Furthermore, recommendations are given how to address possible violations of model assumptions. + +Most plots seen here can also be generated by their dedicated functions, e.g.: + +- Posterior predictive checks: `check_predictions()` +- Homogeneity of variance: `check_heteroskedasticity()` +- Normality of residuals: `check_normality()` +- Multicollinearity: `check_collinearity()` +- Influential observations: `check_outliers()` +- Binned residuals: `binned_residuals()` +- Check for overdispersion: `check_overdispersion()` + +# Logistic regression models: Are all assumptions met? + +We start with a simple example for a logistic regression model. + +```{r} +data(Titanic) +d <- as.data.frame(Titanic) +d <- tidyr::uncount(d, Freq) +m1 <- glm(Survived ~ Class + Sex + Age, data = d, family = binomial()) +``` + +Before we go into details of the diagnostic plots, let's first look at the summary table. + +```{r eval=successfully_loaded["parameters"]} +library(parameters) +model_parameters(m1) +``` + +There is nothing suspicious so far. Now let's start with model diagnostics. We use the `check_model()` function, which provides an overview with the most important and appropriate diagnostic plots for the model under investigation. + +```{r eval=all(successfully_loaded[c("see", "ggplot2")]), fig.height=11} +library(performance) +check_model(m1) +``` + +Now let's take a closer look for each plot. To do so, we ask `check_model()` to return a single plot for each check, instead of arranging them in a grid. We can do so using the `panel` argument. This returns a list of *ggplot* plots. + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# return a list of single plots +diagnostic_plots <- plot(check_model(m1, panel = FALSE)) +``` + +## Posterior predictive checks + +The first plot is based on `check_predictions()`. Posterior predictive checks can be used to "look for systematic discrepancies between real and simulated data" (_Gelman et al. 2014, p. 169_). It helps to see whether the type of model (distributional family) fits well to the data (_Gelman and Hill, 2007, p. 158_). + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# posterior predicive checks +diagnostic_plots[[1]] +``` + +In case of logistic regression our count models, the plot shows by default _dots_ for the observed and simulated data, not _lines_ (as for linear models). The blue dots are simulated data based on the model, if the model were true and distributional assumptions met. The green dots represents the actual observed data of the response variable. + +This plot looks good, because the green dots are inside the range of the blue error bars, and thus we would not assume any violations of model assumptions here. + +### How to fix this? + +The best way, if there are serious concerns that the model does not fit well to the data, is to use a different type (family) of regression models. + +## Binned residuals + + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# linearity +diagnostic_plots[[2]] +``` + + +### How to fix this? + + +## Influential observations - outliers + +Outliers can be defined as particularly influential observations, and this plot helps detecting those outliers. Cook's distance (_Cook 1977_, _Cook & Weisberg 1982_) is used to define outliers, i.e. any point in this plot that falls outside of Cook's distance (the dashed lines) is considered an influential observation. + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# influential observations - outliers +diagnostic_plots[[4]] +``` + +In our example, everything looks well. + +### How to fix this? + +Dealing with outliers is not straightforward, as it is not recommended to automatically discard any observation that has been marked as "an outlier". Rather, your _domain knowledge_ must be involved in the decision whether to keep or omit influential observation. A helpful heuristic is to distinguish between error outliers, interesting outliers, and random outliers (_Leys et al. 2019_). _Error outliers_ are likely due to human error and should be corrected before data analysis. _Interesting outliers_ are not due to technical error and may be of theoretical interest; it might thus be relevant to investigate them further even though they should be removed from the current analysis of interest. _Random outliers_ are assumed to be due to chance alone and to belong to the correct distribution and, therefore, should be retained. + +## Multicollinearity + +This plot checks for potential collinearity among predictors. In a nutshell multicollinearity means that once you know the effect of one predictor, the value of knowing the other predictor is rather low. Multicollinearity might arise when a third, unobserved variable has a causal effect on each of the two predictors that are associated with the outcome. In such cases, the actual relationship that matters would be the association between the unobserved variable and the outcome. + +Multicollinearity should not be confused with a raw strong correlation between predictors. What matters is the association between one or more predictor variables, *conditional on the other variables in the model*. + +If multicollinearity is a problem, the model seems to suggest that the predictors in question don't seems to be reliably associated with the outcome (low estimates, high standard errors), although these predictors actually are strongly associated with the outcome, i.e. indeed might have strong effect (_McElreath 2020, chapter 6.1_). + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# multicollinearity +diagnostic_plots[[5]] +``` + +The variance inflation factor (VIF) indicates the magnitude of multicollinearity of model terms. The thresholds for low, moderate and high collinearity are VIF values less than 5, between 5 and 10 and larger than 10, respectively (_James et al. 2013_). Note that these thresholds, although commonly used, are also criticized for being too high. _Zuur et al. (2010)_ suggest using lower values, e.g. a VIF of 3 or larger may already no longer be considered as "low". + +Our model clearly suffers from multicollinearity, as all predictors have high VIF values. + +### How to fix this? + +Usually, predictors with (very) high VIF values should be removed from the model to fix multicollinearity. Some caution is needed for interaction terms. If interaction terms are included in a model, high VIF values are expected. This portion of multicollinearity among the component terms of an interaction is also called "inessential ill-conditioning", which leads to inflated VIF values that are typically seen for models with interaction terms _(Francoeur 2013)_. In such cases, re-fit your model without interaction terms and check this model for collinearity among predictors. + +## Normality of residuals + +In linear regression, residuals should be normally distributed. This can be checked using so-called Q-Q plots (quantile-quantile plot) to compare the shapes of distributions. This plot shows the quantiles of the studentized residuals versus fitted values. + +Usually, dots should fall along the green reference line. If there is some deviation (mostly at the tails), this indicates that the model doesn't predict the outcome well for the range that shows larger deviations from the reference line. In such cases, inferential statistics like the p-value or coverage of confidence intervals can be inaccurate. + +```{r eval=all(successfully_loaded[c("see", "ggplot2")])} +# normally distributed residuals +diagnostic_plots[[6]] +``` + +In our example, we see that most data points are ok, except some observations at the tails. Whether any action is needed to fix this or not can also depend on the results of the remaining diagnostic plots. If all other plots indicate no violation of assumptions, some deviation of normality, particularly at the tails, can be less critical. + +### How to fix this? + +Here are some remedies to fix non-normality of residuals, according to _Pek et al. 2018_. + +1. For large sample sizes, the assumption of normality can be relaxed due to the central limit theorem - no action needed. + +2. Calculating heteroscedasticity-consistent standard errors can help. See section **Homogeneity of variance** for details. + +3. Bootstrapping is another alternative to resolve issues with non-normally residuals. Again, this can be easily done using the **parameters** package, e.g. `parameters::model_parameters(m1, bootstrap = TRUE)` or [`parameters::bootstrap_parameters()`](https://easystats.github.io/parameters/reference/bootstrap_parameters.html). + +# References + +Brooks ME, Kristensen K, Benthem KJ van, Magnusson A, Berg CW, Nielsen A, et al. glmmTMB Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling. The R Journal. 2017;9: 378-400. + +Cook RD. Detection of influential observation in linear regression. Technometrics. 1977;19(1): 15-18. + +Cook RD and Weisberg S. Residuals and Influence in Regression. London: Chapman and Hall, 1982. + +Francoeur RB. Could Sequential Residual Centering Resolve Low Sensitivity in Moderated Regression? Simulations and Cancer Symptom Clusters. Open Journal of Statistics. 2013:03(06), 24-44. + +Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, and Rubin DB. Bayesian data analysis. (Third edition). CRC Press, 2014 + +Gelman A, Greenland S. Are confidence intervals better termed "uncertainty intervals"? BMJ. 2019;l5381. doi:10.1136/bmj.l5381 + +Gelman A, and Hill J. Data analysis using regression and multilevel/hierarchical models. Cambridge; New York. Cambridge University Press, 2007 + +James, G., Witten, D., Hastie, T., and Tibshirani, R. (eds.).An introduction to statistical learning: with applications in R. New York: Springer, 2013 + +Leys C, Delacre M, Mora YL, Lakens D, Ley C. How to Classify, Detect, and Manage Univariate and Multivariate Outliers, With Emphasis on Pre-Registration. International Review of Social Psychology, 2019 + +McElreath, R. Statistical rethinking: A Bayesian course with examples in R and Stan. 2nd edition. Chapman and Hall/CRC, 2020 + +Pek J, Wong O, Wong ACM. How to Address Non-normality: A Taxonomy of Approaches, Reviewed, and Illustrated. Front Psychol (2018) 9:2104. doi: 10.3389/fpsyg.2018.02104 + +Zuur AF, Ieno EN, Elphick CS. A protocol for data exploration to avoid common statistical problems: Data exploration. Methods in Ecology and Evolution (2010) 1:3-14. diff --git a/_pkgdown.yaml b/_pkgdown.yaml index d6e56740c..3ac6db6fc 100644 --- a/_pkgdown.yaml +++ b/_pkgdown.yaml @@ -44,3 +44,17 @@ reference: - title: "Sample Data" contents: - classify_distribution + +articles: + - title: Checking model assumptions and data properties + navbar: ~ + contents: + - check_model + - check_outliers + - simulate_residuals + + - title: Model comparison and testing + navbar: ~ + contents: + - compare + - r2 diff --git a/man/check_model.Rd b/man/check_model.Rd index 4bdea6ebe..d0d7e765b 100644 --- a/man/check_model.Rd +++ b/man/check_model.Rd @@ -212,12 +212,13 @@ inside the error bounds. See \code{\link[=binned_residuals]{binned_residuals()}} \section{Residuals for (Generalized) Linear Models}{ -Plots that check the normality of residuals (Q-Q plot) or the homogeneity of -variance use standardized Pearson's residuals for generalized linear models, -and standardized residuals for linear models. The plots for the normality of -residuals (with overlayed normal curve) and for the linearity assumption use -the default residuals for \code{lm} and \code{glm} (which are deviance residuals for -\code{glm}). +Plots that check the homogeneity of variance use standardized Pearson's +residuals for generalized linear models, and standardized residuals for +linear models. The plots for the normality of residuals (with overlayed +normal curve) and for the linearity assumption use the default residuals +for \code{lm} and \code{glm} (which are deviance residuals for \code{glm}). The Q-Q plots +use simulated residuals (see \code{\link[=simulate_residuals]{simulate_residuals()}}) for non-Gaussian +models and standardized residuals for linear models. } \section{Troubleshooting}{ diff --git a/man/check_overdispersion.Rd b/man/check_overdispersion.Rd index 8586e7833..692f7aa0b 100644 --- a/man/check_overdispersion.Rd +++ b/man/check_overdispersion.Rd @@ -79,7 +79,7 @@ negative binomial, see \emph{Gelman and Hill (2007), pages 115-116}). \section{Tests based on simulated residuals}{ For certain models, resp. model from certain families, tests are based on -simulated residuals (see \code{\link[=simulated_residual]{simulated_residual()}}). These are usually more +simulated residuals (see \code{\link[=simulate_residuals]{simulate_residuals()}}). These are usually more accurate for testing such models than the traditionally used Pearson residuals. However, when simulating from more complex models, such as mixed models or models with zero-inflation, there are several important considerations. @@ -95,19 +95,10 @@ accurate results. } \examples{ -\dontshow{if (getRversion() >= "4.0.0" && require("glmmTMB", quietly = TRUE)) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} - -library(glmmTMB) -data(Salamanders) +\dontshow{if (getRversion() >= "4.0.0" && require("glmmTMB")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +data(Salamanders, package = "glmmTMB") m <- glm(count ~ spp + mined, family = poisson, data = Salamanders) check_overdispersion(m) - -m <- glmmTMB( - count ~ mined + spp + (1 | site), - family = poisson, - data = Salamanders -) -check_overdispersion(m) \dontshow{\}) # examplesIf} } \references{ diff --git a/man/check_residuals.Rd b/man/check_residuals.Rd index 2c5445686..6935c242c 100644 --- a/man/check_residuals.Rd +++ b/man/check_residuals.Rd @@ -61,5 +61,6 @@ check_residuals(res) \dontshow{\}) # examplesIf} } \seealso{ -\code{\link[=simulate_residuals]{simulate_residuals()}} and \code{\link[=check_predictions]{check_predictions()}}. +\code{\link[=simulate_residuals]{simulate_residuals()}}, \code{\link[=check_zeroinflation]{check_zeroinflation()}}, +\code{\link[=check_overdispersion]{check_overdispersion()}} and \code{\link[=check_predictions]{check_predictions()}}. } diff --git a/man/check_zeroinflation.Rd b/man/check_zeroinflation.Rd index a4ddc7135..7a4da3945 100644 --- a/man/check_zeroinflation.Rd +++ b/man/check_zeroinflation.Rd @@ -54,7 +54,7 @@ is internally called if necessary. \section{Tests based on simulated residuals}{ For certain models, resp. model from certain families, tests are based on -simulated residuals (see \code{\link[=simulated_residual]{simulated_residual()}}). These are usually more +simulated residuals (see \code{\link[=simulate_residuals]{simulate_residuals()}}). These are usually more accurate for testing such models than the traditionally used Pearson residuals. However, when simulating from more complex models, such as mixed models or models with zero-inflation, there are several important considerations. diff --git a/man/simulate_residuals.Rd b/man/simulate_residuals.Rd index fd647f411..938f8b0d2 100644 --- a/man/simulate_residuals.Rd +++ b/man/simulate_residuals.Rd @@ -2,9 +2,12 @@ % Please edit documentation in R/simulate_residuals.R \name{simulate_residuals} \alias{simulate_residuals} +\alias{residuals.performance_simres} \title{Simulate randomized quantile residuals from a model} \usage{ simulate_residuals(x, iterations = 250, ...) + +\method{residuals}{performance_simres}(object, quantile_function = NULL, outlier_values = NULL, ...) } \arguments{ \item{x}{A model object.} @@ -12,6 +15,16 @@ simulate_residuals(x, iterations = 250, ...) \item{iterations}{Number of simulations to run.} \item{...}{Arguments passed on to \code{\link[DHARMa:simulateResiduals]{DHARMa::simulateResiduals()}}.} + +\item{object}{A \code{performance_simres} object, as returned by \code{simulate_residuals()}.} + +\item{quantile_function}{A function to apply to the residuals. If \code{NULL}, the +residuals are returned as is. If not \code{NULL}, the residuals are passed to this +function. This is useful for returning normally distributed residuals, for +example: \code{residuals(x, quantile_function = qnorm)}.} + +\item{outlier_values}{A vector of length 2, specifying the values to replace +\code{-Inf} and \code{Inf} with, respectively.} } \value{ Simulated residuals, which can be further processed with @@ -27,7 +40,7 @@ where the residuals are not expected to be normally distributed. This function is a small wrapper around \code{\link[DHARMa:simulateResiduals]{DHARMa::simulateResiduals()}}. It basically only sets \code{plot = FALSE} and adds an additional class attribute (\code{"performance_sim_res"}), which allows using the DHARMa object in own plotting -functions in the \strong{see} package. See also \code{vignette("DHARMa")}. There is a +functions from the \strong{see} package. See also \code{vignette("DHARMa")}. There is a \code{plot()} method to visualize the distribution of the residuals. } \section{Tests based on simulated residuals}{ @@ -52,6 +65,9 @@ accurate results. \dontshow{if (require("DHARMa")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) simulate_residuals(m) + +# extract residuals +head(residuals(simulate_residuals(m))) \dontshow{\}) # examplesIf} } \references{ @@ -64,5 +80,6 @@ of Computational and Graphical Statistics, 5(3), 236. \doi{10.2307/1390802} } } \seealso{ -\code{\link[=check_residuals]{check_residuals()}} and \code{\link[=check_predictions]{check_predictions()}}. +\code{\link[=check_residuals]{check_residuals()}}, \code{\link[=check_zeroinflation]{check_zeroinflation()}}, +\code{\link[=check_overdispersion]{check_overdispersion()}} and \code{\link[=check_predictions]{check_predictions()}}. } diff --git a/tests/testthat/test-check_collinearity.R b/tests/testthat/test-check_collinearity.R index 4811c5afe..042142073 100644 --- a/tests/testthat/test-check_collinearity.R +++ b/tests/testthat/test-check_collinearity.R @@ -24,6 +24,7 @@ test_that("check_collinearity, correct order in print", { test_that("check_collinearity", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not(getRversion() >= "4.0.0") @@ -50,6 +51,7 @@ test_that("check_collinearity", { test_that("check_collinearity", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not(getRversion() >= "4.0.0") diff --git a/tests/testthat/test-check_convergence.R b/tests/testthat/test-check_convergence.R index 8897b785d..0bfadc8a6 100644 --- a/tests/testthat/test-check_convergence.R +++ b/tests/testthat/test-check_convergence.R @@ -29,6 +29,7 @@ test_that("check_convergence", { test_that("check_convergence, glmmTMB", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") data(iris) model <- suppressWarnings(glmmTMB::glmmTMB( diff --git a/tests/testthat/test-check_model.R b/tests/testthat/test-check_model.R index a0aa2304c..32a78e253 100644 --- a/tests/testthat/test-check_model.R +++ b/tests/testthat/test-check_model.R @@ -53,6 +53,7 @@ test_that("`check_model()` works for quantreg", { }) test_that("`check_model()` warnings for tweedie", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("lme4") data(sleepstudy, package = "lme4") diff --git a/tests/testthat/test-check_normality.R b/tests/testthat/test-check_normality.R index efac29db6..44b461aea 100644 --- a/tests/testthat/test-check_normality.R +++ b/tests/testthat/test-check_normality.R @@ -35,6 +35,7 @@ test_that("check_normality | afex", { }) test_that("check_normality | glmmTMB", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not(getRversion() >= "4.0.0") @@ -82,3 +83,34 @@ test_that("check_normality | t-test", { ignore_attr = TRUE ) }) + + +test_that("check_normality | simulated residuals", { + skip_if_not_installed("DHARMa") + m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) + res <- simulate_residuals(m) + out <- check_normality(res) + expect_equal( + as.numeric(out), + 0.2969038, + tolerance = 1e-3, + ignore_attr = TRUE + ) + expect_identical( + capture.output(print(out)), + "OK: residuals appear as normally distributed (p = 0.297)." + ) + + m <- lm(mpg ~ wt + cyl + gear + disp, data = mtcars) + out <- check_normality(m) + expect_equal( + as.numeric(out), + 0.2303071, + tolerance = 1e-3, + ignore_attr = TRUE + ) + expect_identical( + capture.output(print(out)), + "OK: residuals appear as normally distributed (p = 0.230)." + ) +}) diff --git a/tests/testthat/test-check_outliers.R b/tests/testthat/test-check_outliers.R index 016c827af..704b27feb 100644 --- a/tests/testthat/test-check_outliers.R +++ b/tests/testthat/test-check_outliers.R @@ -6,6 +6,7 @@ test_that("zscore negative threshold", { }) test_that("lof illegal threshold", { + skip_if_not_installed("dbscan") expect_error( check_outliers(mtcars$mpg, method = "lof", threshold = -1), "The `threshold` argument" @@ -84,6 +85,7 @@ test_that("mahalanobis_robust which", { }) test_that("mcd which", { + skip_if_not_installed("MASS") # (not clear why method mcd needs a seed) set.seed(42) expect_identical( @@ -199,6 +201,12 @@ test_that("multiple methods which", { # We exclude method ics because it is too slow test_that("all methods which", { skip_if_not_installed("bigutilsr") + skip_if_not_installed("MASS") + skip_if_not_installed("dbscan") + skip_if_not_installed("ICS") + skip_if_not_installed("ICSOutlier") + skip_if_not_installed("loo") + expect_identical( which(check_outliers(mtcars, method = c( @@ -222,6 +230,12 @@ test_that("all methods which", { test_that("multiple methods with ID", { skip_if_not_installed("bigutilsr") + skip_if_not_installed("MASS") + skip_if_not_installed("dbscan") + skip_if_not_installed("ICS") + skip_if_not_installed("ICSOutlier") + skip_if_not_installed("loo") + data <- datawizard::rownames_as_column(mtcars, var = "car") x <- attributes(check_outliers(data, method = c( @@ -280,6 +294,7 @@ test_that("cook multiple methods which", { test_that("pareto which", { skip_if_not_installed("dbscan") + skip_if_not_installed("loo") skip_if_not_installed("rstanarm") set.seed(123) model <- rstanarm::stan_glm(mpg ~ qsec + wt, data = mtcars, refresh = 0) @@ -293,6 +308,7 @@ test_that("pareto which", { test_that("pareto multiple methods which", { skip_if_not_installed("dbscan") + skip_if_not_installed("loo") skip_if_not_installed("rstanarm") set.seed(123) model <- rstanarm::stan_glm(mpg ~ qsec + wt, data = mtcars, refresh = 0) diff --git a/tests/testthat/test-check_overdispersion.R b/tests/testthat/test-check_overdispersion.R index 2930322bb..06cc95dd0 100644 --- a/tests/testthat/test-check_overdispersion.R +++ b/tests/testthat/test-check_overdispersion.R @@ -1,4 +1,5 @@ test_that("check_overdispersion, glmmTMB-poisson", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not(getRversion() >= "4.0.0") data(Salamanders, package = "glmmTMB") @@ -49,6 +50,7 @@ test_that("check_overdispersion, glmmTMB-poisson", { test_that("check_overdispersion, glmmTMB-poisson mixed", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not(getRversion() >= "4.0.0") data(Salamanders, package = "glmmTMB") @@ -76,6 +78,7 @@ test_that("check_overdispersion, glmmTMB-poisson mixed", { test_that("check_overdispersion, zero-inflated and negbin", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("DHARMa") skip_if_not(getRversion() >= "4.0.0") @@ -180,6 +183,7 @@ test_that("check_overdispersion, MASS::negbin", { test_that("check_overdispersion, genpois", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("DHARMa") skip_if_not(getRversion() >= "4.0.0") diff --git a/tests/testthat/test-check_predictions.R b/tests/testthat/test-check_predictions.R index a840894dc..1f7774c8a 100644 --- a/tests/testthat/test-check_predictions.R +++ b/tests/testthat/test-check_predictions.R @@ -35,6 +35,7 @@ test_that("check_predictions", { test_that("check_predictions, glmmTMB", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") data(mtcars) model <- glmmTMB::glmmTMB(vs ~ disp, data = mtcars, family = binomial()) diff --git a/tests/testthat/test-check_residuals.R b/tests/testthat/test-check_residuals.R index 33407abf9..9be029c0c 100644 --- a/tests/testthat/test-check_residuals.R +++ b/tests/testthat/test-check_residuals.R @@ -1,10 +1,61 @@ -test_that("check_singularity, lme4", { +test_that("check_residuals and simulate_residuals", { skip_on_cran() skip_if_not_installed("DHARMa") set.seed(123) dat <- DHARMa::createData(sampleSize = 100, overdispersion = 0.5, family = poisson()) m <- glm(observedResponse ~ Environment1, family = poisson(), data = dat) res <- simulate_residuals(m) + expect_identical( + capture.output(print(res)), + c( + "Simulated residuals from a model of class `glm` based on 250", + " simulations. Use `check_residuals()` to check uniformity of residuals or", + " `residuals()` to extract simulated residuals. It is recommended to refer", + " to `?DHARMa::simulateResiudals` and `vignette(\"DHARMa\")` for more", + " information about different settings in particular situations or for", + " particular models." + ) + ) + # check raw residuals + expect_equal( + head(residuals(res)), + c(0.55349, 0.44012, 0.39826, 0.9825, 0.90753, 0.05809), + tolerance = 1e-4, + ignore_attr = TRUE + ) + expect_equal( + head(residuals(res, quantile_function = stats::qnorm)), + c(0.13448, -0.15068, -0.25785, 2.10826, 1.3257, -1.57097), + tolerance = 1e-4, + ignore_attr = TRUE + ) + # compare to DHARMa + res_d <- DHARMa::simulateResiduals(m, n = 250, plot = FALSE) + expect_equal( + head(residuals(res)), + head(residuals(res_d)), + tolerance = 1e-4, + ignore_attr = TRUE + ) + expect_equal( + head(residuals(res, quantile_function = stats::qnorm)), + head(residuals(res_d, quantileFunction = stats::qnorm)), + tolerance = 1e-4, + ignore_attr = TRUE + ) + # DHARMa args work in residuals.permormance_simres + expect_equal( + residuals(res, quantileFunction = stats::qnorm, outlierValues = c(-3, 3)), + residuals(res_d, quantileFunction = stats::qnorm, outlierValues = c(-3, 3)), + tolerance = 1e-4, + ignore_attr = TRUE + ) + # outlier_values works + expect_identical(sum(is.infinite(residuals(res, quantile_function = stats::qnorm))), 3L) + expect_identical(sum(is.infinite(residuals(res, quantile_function = stats::qnorm, outlier_values = c(-100, 100)))), 0L) # nolint + expect_error(residuals(res, quantile_function = stats::qnorm, outlier_values = 1:3), regex = "`outlier_values` must be") # nolint + + # check_residuals out <- check_residuals(res) expect_equal(out, 0.01884602, ignore_attr = TRUE, tolerance = 1e-4) expect_identical( diff --git a/tests/testthat/test-check_zeroinflation.R b/tests/testthat/test-check_zeroinflation.R index 38a5c7726..a7ff7cbbd 100644 --- a/tests/testthat/test-check_zeroinflation.R +++ b/tests/testthat/test-check_zeroinflation.R @@ -1,4 +1,5 @@ test_that("check_zeroinflation", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") set.seed(123) data(Salamanders, package = "glmmTMB") @@ -21,6 +22,7 @@ test_that("check_zeroinflation", { test_that("check_zeroinflation, glmmTMB with and without zero-inflation component", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("DHARMa") set.seed(123) @@ -107,6 +109,7 @@ test_that("check_zeroinflation, glmer.nb", { test_that("check_zeroinflation, glmmTMB nbinom", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("DHARMa") skip_on_cran() @@ -163,6 +166,7 @@ test_that("check_zeroinflation, MASS::negbin", { test_that("check_zeroinflation, genpois", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") skip_if_not_installed("DHARMa") skip_if_not(getRversion() >= "4.0.0") diff --git a/tests/testthat/test-icc.R b/tests/testthat/test-icc.R index ad4e4b033..68624ff09 100644 --- a/tests/testthat/test-icc.R +++ b/tests/testthat/test-icc.R @@ -90,12 +90,12 @@ test_that("icc", { skip_if_not_installed("lme4") data(sleepstudy, package = "lme4") set.seed(12345) - sleepstudy$grp <- sample(1:5, size = 180, replace = TRUE) + sleepstudy$grp <- sample.int(5, size = 180, replace = TRUE) sleepstudy$subgrp <- NA for (i in 1:5) { filter_group <- sleepstudy$grp == i sleepstudy$subgrp[filter_group] <- - sample(1:30, size = sum(filter_group), replace = TRUE) + sample.int(30, size = sum(filter_group), replace = TRUE) } model <- lme4::lmer( Reaction ~ Days + (1 | grp) + (1 | Subject), @@ -125,3 +125,50 @@ test_that("icc", { expect_equal(out$ICC_adjusted, 0.9104331, tolerance = 0.01) expect_equal(out$ICC_unadjusted, 0.3109478, tolerance = 0.01) }) + + +test_that("icc, glmmTMB 1.1.9+", { + skip_on_cran() + skip_if_not_installed("glmmTMB", minimum_version = "1.1.9") + set.seed(101) + dd <- data.frame( + z = rnorm(1000), + x1 = 1:1000, + x2 = runif(1000, 0, 10), + re = rep(1:20, each = 50) + ) + dd <- transform(dd, x3 = as.factor(ifelse( + x1 <= 500, "Low", sample(c("Middle", "High"), 1000, replace = TRUE) + ))) + dd <- transform(dd, x4 = as.factor(ifelse( + x1 > 500, "High", sample(c("Absent", "Low"), 1000, replace = TRUE) + ))) + dd <- transform(dd, z = z + re * 5) + expect_message({ + mod_TMB <- glmmTMB::glmmTMB( + z ~ x1 + x2 + x3 + x4 + (1 | re), + data = dd, + start = list(theta = 3), + control = glmmTMB::glmmTMBControl(rank_check = "adjust") + ) + }) + expect_equal( + icc(mod_TMB), + data.frame( + ICC_adjusted = 0.995480998331767, + ICC_conditional = 0.244468078371849, + ICC_unadjusted = 0.244468078371849 + ), + ignore_attr = TRUE, + tolerance = 1e-4 + ) + expect_equal( + r2(mod_TMB), + list( + R2_conditional = c(`Conditional R2` = 0.998890233308478), + R2_marginal = c(`Marginal R2` = 0.754422154936629) + ), + ignore_attr = TRUE, + tolerance = 1e-4 + ) +}) diff --git a/tests/testthat/test-r2.R b/tests/testthat/test-r2.R index d287604a1..f521b3aab 100644 --- a/tests/testthat/test-r2.R +++ b/tests/testthat/test-r2.R @@ -46,6 +46,7 @@ skip_if_not_installed("withr") withr::with_environment( new.env(), test_that("r2 glmmTMB, no ranef", { + skip_if(getRversion() > "4.3.3") skip_if_not_installed("glmmTMB") data(Owls, package = "glmmTMB") # linear --------------------------------------------------------------- diff --git a/vignettes/simulate_residuals.Rmd b/vignettes/simulate_residuals.Rmd index 7338aea4e..f1a627147 100644 --- a/vignettes/simulate_residuals.Rmd +++ b/vignettes/simulate_residuals.Rmd @@ -27,7 +27,7 @@ knitr::opts_chunk$set( ) options(digits = 2) -pkgs <- c("DHARMa", "glmmTMB") +pkgs <- c("DHARMa", "glmmTMB", "see") successfully_loaded <- vapply(pkgs, requireNamespace, FUN.VALUE = logical(1L), quietly = TRUE) can_evaluate <- all(successfully_loaded) @@ -61,13 +61,18 @@ simulated_residuals <- simulate_residuals(model) simulated_residuals ``` +The raw residuals can be extracted using `residuals()`: + +```{r} +head(residuals(simulated_residuals)) +``` + + Note that since this inherits the DHARMa class, all the methods implemented in DHARMa just work, including all the tests: ```{r} library(DHARMa) -residuals(simulated_residuals) - DHARMa::testUniformity(simulated_residuals, plot = FALSE) ``` @@ -78,17 +83,26 @@ Finally, run specific checks on the simulated residuals: check_residuals(simulated_residuals) ``` -Or check the entire model. +Further implemented checks are tests for overdispersion, outliers and zero-inflation. + +```{r} +check_overdispersion(simulated_residuals) + +check_zeroinflation(simulated_residuals) + +check_outliers(simulated_residuals) +``` + +The above three functions internally call `simulate_residuals()` for more complex models automatically, so you don't need to call `simulate_residuals()` yourself. Simulated residuals are usually more reliable than the standard residuals, especially for complex models. + +Finally, you can even perform a visual check for the entire model, either by passing the model object directly, or the object returned from `simulate_residuals()`. -```{r, eval=FALSE} -# TODO (not implemented) +```{r, eval=can_evaluate && packageVersion("see") > "0.8.2"} check_model(simulated_residuals) ``` -The `check_model()` function is the main reason we don't want to prematurely extract the residuals in `simulate_residuals()`, because if we do then the `simulated_residuals` won't contain the model fit (`fittedModel` in the output below), so we won't be able to do all of the checks we would want to do using the model (e.g., posterior predictive checks). +The `check_model()` function is the main reason we don't want to prematurely extract the residuals in `simulate_residuals()`, because if we do then the simulated residual won't contain the model fit (`fittedModel` in the output below), so we won't be able to do all of the checks we would want to do using the model (e.g., posterior predictive checks). ```{r} str(simulated_residuals, max.level = 1) ``` - -It would also mean we would need to reimplement some of the tests from DHARMa (e.g., `DHARMa::testOutliers()`) if we're planning to include those checks as well. We probably don't want to do that, since some of them are fairly involved rather than just being wrappers for tests supplied in base R (e.g., ) .