DHARMa: Residual Diagnostics for HierARchical Models

repository·master·Indexed 42 days ago

https://github.com/florianhartig/dharma

An R package for residual diagnostics of complex models, including hierarchical and generalized linear models (GLMs, GLMMs). It uses a simulation-based approach to transform residuals into standardized quantile residuals. Supports model classes from lme4, glmmTMB, GLMMadaptive, spaMM, mgcv, phylolm, and Bayesian simulations from JAGS, STAN, or BUGS. Provides tools for detecting over/underdispersion, zero-inflation, and spatial, temporal, or phylogenetic autocorrelation.

Tokens
7.5K
Snippets
31
Records
41
Agent score
53%

What's inside DHARMa

  1. Overview of DHARMa residual diagnostics

    master

    DHARMa (Residual Diagnostics for HierARchical Models) provides a simulation-based approach to create scaled (quantile) residuals for fitted generalized linear (mixed) models. These residuals are standardized to values between 0 and 1, making them as intuitively interpretable as residuals from a standard linear regression.

    Supported Model Classes:

    • Linear and GLMMs: lme4 (lmerMod, glmerMod), glmmTMB, GLMMadaptive, and spaMM.
    • Phylogenetic Models: phylolm (phylolm, phyloglm).
    • Additive Models: mgcv (gam).
    • GLMs: glm (including negbin from MASS, but excluding quasi-distributions) and lm.
    • Bayesian Simulations: Externally created simulations (e.g., posterior predictive simulations from JAGS, STAN, or BUGS).

    Key Capabilities:

    • Standardized quantile residuals.
    • Diagnostic tests and plots for model misspecification, including:
      • Over/underdispersion
      • Zero-inflation
      • Residual spatial, temporal, and phylogenetic autocorrelation.
  2. Handle glmmTMB binomial model limitations in DHARMa

    master

    When using glmmTMB with binomial families, simulateResiduals() may trigger a warning regarding model predictions being calculated conditional on fitted random effects. This can result in asymmetries in 'residual vs. predicted' plots.

    If you observe asymmetries in your plots but the residuals otherwise appear correct, you should calculate the model prediction manually using only fixed effects and use those for your plots.

    Additionally, when working with binomial models where the response is provided as an $n/k$ matrix (e.g., cbind(successes, failures)), ensure your data is formatted correctly for glmmTMB before passing it to simulateResiduals().

    # Example of a binomial model with glmmTMB
    newDat$response = cbind(newDat$incidence, newDat$size - newDat$incidence)
    m <- glmmTMB(response ~ period, data = newDat, family = binomial)
    res = simulateResiduals(m)
    plot(res, rank = TRUE)
  3. Configure marginal vs. conditional simulations

    master

    By default, simulateResiduals() performs marginal simulations, re-simulating the entire model structure including all Random Effects (REs).

    In packages like lme4, you can perform conditional simulations (conditioning on the fitted REs) using the re.form argument. Conditioning on REs can increase the power of certain tests, such as dispersion tests.

    # Marginal simulation (default)
    res <- simulateResiduals(fittedModel = fittedModel)
    
    # Conditional simulation (conditioning on REs in lme4)
    res <- simulateResiduals(fittedModel = fittedModel, re.form = NULL)
  4. Conditional vs Unconditional posterior predictive simulations

    master

    When performing residual checks on complex Bayesian models (like GLMMs with random effects or zero-inflation), you can generate two types of simulations:

    Conditional Simulations

    These use the specific fitted values for random effects and zero-inflation parameters estimated for each observation. This tests how well the model describes the data given the estimated latent structure.

    Unconditional Simulations

    These re-sample the random effects and zero-inflation terms from their respective posterior distributions for every simulation. This tests the overall model structure and its ability to generate new data that resembles the original process.

    Note on Unconditional Residuals: When plotting unconditional residuals, the fittedPredictedResponse should be adjusted (e.g., using the mean of the simulations) because the standard conditional predictions may not align with the re-sampled structure.

  5. Understanding the GAM binomial simulation issue in mgcv

    master

    In certain versions of the mgcv package (specifically around version 1.8-14), there is a known issue when fitting binomial Generalized Additive Models (GAMs) using $n/k$ (observed counts) data.

    When using gam() with a binomial family, the fitted() function returns an unnamed vector. This causes the family$simulate() function to produce output with incorrect dimensions (e.g., a dimension of [1] 0 10 instead of the expected [1] 200 10), which can break downstream simulation processes.

    Note: This issue was addressed in DHARMa > 0.1.3 via a compatibility hack. If you encounter dimension mismatches when simulating from a binomial GAM, ensure you are using a recent version of DHARMa.

  6. Run reverse dependency checks for DHARMa

    master

    Before submitting a release to CRAN, you should check how changes affect packages that depend on DHARMa. Use the revdepcheck package to automate this process. This will create a revdep/ directory containing the results of the checks. A high-level summary is generated in cran.md, which can be copied into cran-comments.md for the CRAN submission.

    To set up and run the checks, use the following workflow:

    # Setup the environment for reverse dependency checks
    usethis::use_revdep() 
    
    # Run the checks using multiple workers
    revdepcheck::revdep_check(num_workers=4)
  7. Visualize residuals with the main DHARMa plot

    master

    Use plot(res) to generate the standard DHARMa diagnostic plot. The plot consists of two panels:

    Left Panel (Uniformity & Dispersion):

    • Uniform QQ plot: Checks if residuals follow a uniform distribution.
    • KS-test: Kolmogorov-Smirnov test for uniformity.
    • Dispersion test: Compares observed variance to simulated variance.
    • Outlier tests: Checks if the number of observations outside the simulation envelope is unexpected.

    Right Panel (Dependency):

    • res ~ predicted: Plots residuals against predicted values. A well-specified model should show a uniform distribution in the y-direction.
    • Quantile GAMs: Fits GAMs to the 0.25, 0.5, and 0.75 quantiles. Significant deviations from a straight line indicate potential issues.
    plot(res)
  8. Key changes in DHARMa v0.5.0

    master

    Version 0.5.0 introduces several significant changes to how DHARMa handles residual simulations and supported models:

    • Conditional Simulations for GLMMs: The default simulation method for Generalized Linear Mixed Models (GLMMs) has changed from the supported model's default (which was mostly unconditional) to conditional simulations. This applies to most supported models and packages.
    • brms Support: Added support for the brms package.
    • Plotting Improvements: Various improvements to plotting functions.
    • Bugfixes: Includes various smaller improvements and bugfixes.
  9. Check residuals against predictors

    master

    While the main plot(res) checks residuals against predicted values, it is a best practice to explicitly check residuals against all individual predictors using plotResiduals(). This can reveal issues (like missing predictors) that might not be obvious in the standard plot.

    # Check residuals against a specific predictor
    plotResiduals(res, form = testData$Environment1)
  10. Use DHARMa with glmmTMB models

    master

    To perform residual diagnostics on models fitted with glmmTMB, use the simulateResiduals() function. This works for various model types, including Poisson and Zero-Inflated Negative Binomial (ZINB) models. After generating residuals, you can visualize them using plot().

    Note on glmmTMB limitations: When using glmmTMB, simulateResiduals() calculates model predictions conditional on the fitted random effects. This can sometimes cause asymmetries in plots of residuals vs. predicted values. If you observe such asymmetries but the residuals otherwise appear correct, consider calculating model predictions manually using only fixed effects for plotting purposes.

    # Example: Zero-inflated negative binomial model
    m <- glmmTMB(count ~ 1, zi = ~1, family = nbinom1, Salamanders)
    
    # Generate residuals
    res <- simulateResiduals(m)
    
    # Plot residuals
    plot(res, asFactor = TRUE)
  11. Install DHARMa development release from GitHub

    master

    If you need new features that have not yet been released to CRAN (typically available 1-2 months earlier), you can install the development version directly from the GitHub repository using devtools.

    devtools::install_github(repo = "florianhartig/DHARMa", subdir = "DHARMa", 
                             dependencies = T, build_vignettes = T)