simulateLRT: Simulated likelihood ratio tests for (generalized) linear...

In DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models

Simulated likelihood ratio tests for (generalized) linear mixed models

Description

This function uses the DHARMa model wrappers to generate simulated likelihood ratio tests (LRTs) for (generalized) linear mixed models based on a parametric bootstrap. The motivation for using a simulated LRT rather than a standard ANOVA or AIC for model selection in mixed models is that df for mixed models are not clearly defined, thus standard ANOVA based on Chi2 statistics or AIC are unreliable, in particular for models with large contributions of REs to the likelihood.

Interpretation of the results as in a normal LRT: the null hypothesis is that m0 is correct, the tests checks if the increase in likelihood of m1 is higher than expected, using data simulated from m0.

Usage

simulateLRT(m0, m1, n = 250, seed = 123, plot = TRUE,
  suppressWarnings = TRUE, saveModels = FALSE, ...)

Arguments

m0	null Model.
m1	alternative Model.
n	number of simulations.
seed	random seed.
plot	whether null distribution should be plotted.
suppressWarnings	whether to suppress warnings that occur during refitting the models to simulated data. See details for explanations.
saveModels	Whether to save refitted models.
...	additional parameters to pass on to the simulate function of the model object. See getSimulations for details.

Details

The function performs a simulated LRT, which works as follows:

1. H0: Model 1 is correct.

2. Our test statistic is the log LRT of M1/M2. Empirical value will always be > 1 because in a nested setting, the more complex model cannot have a worse likelihood.

3. To generate an expected distribution of the test statistic under H0, we simulate new response data under M0, refit M0 and M1 on this data, and calculate the LRs.

4. Based on this, calculate p-values etc. in the usual way.

About warnings: warnings such as "boundary (singular) fit: see ?isSingular" will likely occur in this function and are not necessarily the sign of a problem. lme4 warns if RE variances are fit to zero. This is desired / likely in this case, however, because we are simulating data with zero RE variances. Therefore, warnings are turned off per default. For diagnostic reasons, you can turn warnings on, and possibly also inspect fitted models via the parameter saveModels to see if there are any other problems in the re-fitted models.

Data simulations are performed by getSimulations, which is a wrapper for the respective model functions. The default for all packages, wherever possible, is to generate marginal simulations (meaning that REs are re-simulated as well). I see no sensible reason to change this, but if you want to and if supported by the respective regression package, you could do so by supplying the necessary arguments via ...

Note

The logic of an LRT assumes that m0 is nested in m1, which guarantees that the L(M1) > L(M0). The function does not explicitly check if models are nested and will work as long as data can be simulated from M0 that can be refit with M) and M1; however, I would strongly advice against using this for non-nested models unless you have a good statistical reason for doing so.

Also, note that LRTs may be unreliable when fit with REML or some other kind of penalized / restricted ML. Therefore, you should fit model with ML for use in this function.

Author(s)

Florian Hartig

Examples

library(DHARMa)
library(lme4)

# create test data
set.seed(123)
dat <- createData(sampleSize = 200, randomEffectVariance = 1)

# define Null and alternative model (should be nested)
m1 = glmer(observedResponse ~ Environment1 + (1|group), data = dat, family = "poisson")
m0 = glm(observedResponse ~ Environment1 , data = dat, family = "poisson")

## Not run: 
# run LRT - n should be increased to at least 250 for a real study
out = simulateLRT(m0, m1, n = 10)

# To inspect warnings thrown during the refits:
out = simulateLRT(m0, m1, saveModels = TRUE, suppressWarnings = FALSE, n = 10)
summary(out$saveModels[[2]]$refittedM1) # RE SD = 0, no problem
# If there are warnings that seem problematic, 
# could try changing the optimizer or iterations

## End(Not run)

DHARMa documentation built on Oct. 18, 2024, 5:09 p.m.