-
-
Notifications
You must be signed in to change notification settings - Fork 87
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #24 from easystats/dev
use insight::get_variances, add overdispersion() and zerocount() funs
- Loading branch information
Showing
38 changed files
with
695 additions
and
634 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file was deleted.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,57 @@ | ||
| #' @title Convergence test for mixed effects models | ||
| #' @name check_convergence | ||
| #' | ||
| #' @description \code{check_convergence()} provides an alternative convergence | ||
| #' test for \code{\link[lme4]{merMod}}-objects. | ||
| #' | ||
| #' @param x A \code{merMod}-object. | ||
| #' @param tolerance Indicates up to which value the convergence result is | ||
| #' accepted. The smaller \code{tolerance} is, the stricter the test | ||
| #' will be. | ||
| #' | ||
| #' @return \code{TRUE} if convergence is fine and \code{FALSE} if convergence | ||
| #' is suspicious. Additionally, the convergence value is returned as return | ||
| #' value's name. | ||
| #' | ||
| #' @details \code{check_convergence()} provides an alternative convergence test for | ||
| #' \code{\link[lme4]{merMod}}-objects, as discussed | ||
| #' \href{https://github.com/lme4/lme4/issues/120}{here} | ||
| #' and suggested by Ben Bolker in | ||
| #' \href{https://github.com/lme4/lme4/issues/120#issuecomment-39920269}{this comment}. | ||
| #' | ||
| #' @examples | ||
| #' library(lme4) | ||
| #' data(cbpp) | ||
| #' set.seed(1) | ||
| #' cbpp$x <- rnorm(nrow(cbpp)) | ||
| #' cbpp$x2 <- runif(nrow(cbpp)) | ||
| #' | ||
| #' model <- glmer( | ||
| #' cbind(incidence, size - incidence) ~ period + x + x2 + (1 + x | herd), | ||
| #' data = cbpp, | ||
| #' family = binomial() | ||
| #' ) | ||
| #' | ||
| #' check_convergence(model) | ||
| #' | ||
| #' @export | ||
| check_convergence <- function(x, tolerance = 0.001) { | ||
| # check for package availability | ||
| if (!requireNamespace("Matrix", quietly = TRUE)) | ||
| stop("Package `Matrix` needed for this function to work. Please install it.", call. = FALSE) | ||
|
|
||
| # is 'x' an lmer object? | ||
| if (!inherits(x, "merMod")) | ||
| warning("`x` must be a `merMod` object.", call. = F) | ||
|
|
||
|
|
||
| relgrad <- with(x@optinfo$derivs, Matrix::solve(Hessian, gradient)) | ||
|
|
||
| # copy logical value, TRUE if convergence is OK | ||
| retval <- max(abs(relgrad)) < tolerance | ||
| # copy convergence value | ||
| names(retval) <- max(abs(relgrad)) | ||
|
|
||
| # return result | ||
| retval | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,97 @@ | ||
| #' @title Check overdispersion of GL(M)M's | ||
| #' @name check_overdispersion | ||
| #' | ||
| #' @description \code{check_overdispersion()} checks generalized linear (mixed) models | ||
| #' for overdispersion. | ||
| #' | ||
| #' @param x Fitted model of class \code{merMod}, \code{glmmTMB}, \code{glm}, | ||
| #' or \code{glm.nb} (package \pkg{MASS}). | ||
| #' @param ... Currently not used. | ||
| #' | ||
| #' @return Results from the on the overdispersion test, like chi-squared statistics, | ||
| #' p-value or dispersion ratio. | ||
| #' | ||
| #' @details A p-value < .05 indicates overdispersion. For \code{merMod}- and | ||
| #' \code{glmmTMB}-objects, \code{check_overdispersion()} is based on the code in | ||
| #' the \href{http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html}{GLMM FAQ}, | ||
| #' section \emph{How can I deal with overdispersion in GLMMs?}. Note that | ||
| #' this function only returns an \emph{approximate} estimate of an | ||
| #' overdispersion parameter, and is probably inaccurate for zero-inflated | ||
| #' mixed models (fitted with \code{glmmTMB}). The same code is also used to | ||
| #' check overdispersion for negative binomial models. | ||
| #' \cr \cr | ||
| #' For Poisson-models, the overdispersion test is based on the code | ||
| #' from \cite{Gelman and Hill (2007), page 115}. | ||
| #' | ||
| #' @references Bolker B et al. (2017): \href{http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html}{GLMM FAQ.} | ||
| #' \cr \cr | ||
| #' Gelman A, Hill J (2007) Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge, New York: Cambridge University Press | ||
| #' | ||
| #' @examples | ||
| #' | ||
| #' @export | ||
| check_overdispersion <- function(x, ...) { | ||
| UseMethod("check_overdispersion") | ||
| } | ||
|
|
||
|
|
||
| #' @importFrom insight get_response | ||
| #' @importFrom stats fitted nobs coef pchisq | ||
| #' @export | ||
| check_overdispersion.glm <- function(x, ...) { | ||
| # check if we have poisson | ||
| if (!stats::family(x)$family %in% c("poisson", "quasipoisson")) | ||
| stop("Model must be from Poisson-family.", call. = F) | ||
|
|
||
| yhat <- stats::fitted(x) | ||
|
|
||
| n <- stats::nobs(x) | ||
| k <- length(stats::coef(x)) | ||
|
|
||
| zi <- (insight::get_response(x) - yhat) / sqrt(yhat) | ||
| chisq <- sum(zi^2) | ||
| ratio <- chisq / (n - k) | ||
| p.value <- stats::pchisq(chisq, df = n - k, lower.tail = FALSE) | ||
|
|
||
| list( | ||
| chisq_statistic = chisq, | ||
| dispersion_ratio = ratio, | ||
| residual_df = n - k, | ||
| p_value = p.value | ||
| ) | ||
| } | ||
|
|
||
|
|
||
| #' @export | ||
| check_overdispersion.negbin <- function(x, ...) { | ||
| check_overdispersion.lme4(x) | ||
| } | ||
|
|
||
|
|
||
| #' @export | ||
| check_overdispersion.merMod <- function(x, ...) { | ||
| check_overdispersion.lme4(x) | ||
| } | ||
|
|
||
|
|
||
| #' @export | ||
| check_overdispersion.glmmTMB <- function(x, ...) { | ||
| check_overdispersion.lme4(x) | ||
| } | ||
|
|
||
|
|
||
| #' @importFrom stats df.residual residuals pchisq | ||
| check_overdispersion.lme4 <- function(x) { | ||
| rdf <- stats::df.residual(x) | ||
| rp <- stats::residuals(x, type = "pearson") | ||
| Pearson.chisq <- sum(rp^2) | ||
| prat <- Pearson.chisq / rdf | ||
| pval <- stats::pchisq(Pearson.chisq, df = rdf, lower.tail = FALSE) | ||
|
|
||
| list( | ||
| chisq = Pearson.chisq, | ||
| ratio = prat, | ||
| rdf = rdf, | ||
| p = pval | ||
| ) | ||
| } |
Oops, something went wrong.