LRT: Tests of fixed effects (LRTs and ANOVA tables).

Description

* LRT performs a likelihood ratio (LR) test between two model fits, the ``full'' and the ``null'' model fits. If the models differ only in their fixed effects, a standard test based on the asymptotic chi-square distribution is performed, with number of degrees of freedom determined by the function. In addition, parametric bootstrap p-values can be computed, either using the raw bootstrap distribution of the likelihood ratio, or a bootstrap estimate of the Bartlett correction of the LR statistic. This function differs from fixedLRT in its arguments (model fits for LRT, but all arguments required to fit the models for fixedLRT), and in the format of its return value.

If the two models differ in their random effects, a bootstrap test may be performed, and no number of degrees of freedom is deduced, so no asymptotic test is performed. Distinguishing the full and the null model from random-effect specifications by a simple comparison of the model formulas is not always easy, so in some case the model with the lower likelihood is assumed to be the null one (the latter comparison is subject to numerical uncertainties when both fits are equivalent).

If the two models differ in both their fixed and random components, the bootstrap test can also be performed (see Examples), but the procedure further checks that the same model is nested in the other for both components. This requires that a simple comparison of the model formulas is sufficient to assess nestedness for random effects, and may therefore fail.

If the two models differ neither in their fixed nor random components, the residual dispersion models are tentatively compared by then number of degrees of freedom and an analyses may be performed if they are both residual-dispersion models. Check of nestedness are incomplete in this case. Mixed-effect residual models are not yet handled.

* The anova method for fit objects from spaMM has two uses: if a single fit object is provided, ANOVA tables may be returned, with specific procedures for univariate-response LMs, GLMs and LMMs (see Details). Alternatively, if a second fit object is provided (object2 argument), anova performs as an alias for LRT.

Usage

# S3 method for HLfit
anova(object, object2, type = "2", method="", ...)
#
LRT(object, object2, boot.repl = 0, resp_testfn = NULL, 
    simuland = eval_replicate, 
    #     many further arguments can be passed to spaMM_boot via the '...'
    #     These include arguments for parallel computations, such as
    # nb_cores, fit_env,
    #     as well as other named arguments and spaMM_boot's own '...'
    ...)

Value

LRT returns an object of class fixedLRT, actually a list with typical elements (depending on the options)

fullfit: the HLfit object for the full model;
nullfit: the HLfit object for the null model;
basicLRT: A data frame including values of the likelihood ratio chi2 statistic, its degrees of freedom, and the p-value;

and, if a bootstrap was performed:

rawBootLRT

A data frame including values of the likelihood ratio chi2 statistic, its degrees of freedom, and the raw bootstrap p-value;

BartBootLRT

A data frame including values of the Bartlett-corrected likelihood ratio chi2 statistic, its degrees of freedom, and its p-value;

bootInfo

a list with the following elements:

bootreps: A table of fitted likelihoods for bootstrap replicates;

meanbootLRT

The mean likelihood ratio chi-square statistic for bootstrap replicates;

When ANOVA tables are computed, the return format is that of the function called (lmerTest::anova for LMMs) or emulated (for LMs or GLMs).

Arguments

object: Fit object returned by a spaMM fitting function.
object2: Optional second model fit to be be compared to the first (their order does not matter).
type: ANOVA type for LMMs. Note that the default (single-term deletion ANOVA) differs from that of lmerTest.

boot.repl: the number of bootstrap replicates.
resp_testfn: See argument resp_testfn of spaMM_boot.
simuland: a function, passed to spaMM_boot. See argument eval_replicate for default value and requirements.
method: Only non-default value is "t.Chisq" which forces evaluation of a table of chi-squared tests for each fixed-effect term, using the classical “Wald” test (see Details).

...: Further arguments, passed to spaMM_boot (e.g., for parallelization) in the case of LRTs. For ANOVA tables, arguments of functions anova.lm anova.glm, and as_LMLT, respectively for LMs, GLMs and LMMs, may be handled (e.g. the test argument for anova.glm).

Details

* Bootstrap LRTs: A raw bootstrap p-value can be computed from the simulated distribution as (1+sum(t >= t0))/(N+1) where t0 is the original likelihood ratio, t the vector of bootstrap replicates and N its length. See Davison & Hinkley (1997, p. 141) for discussion of the adjustments in this formula. However, a computationally more economical use of the bootstrap is to provide a Bartlett correction for the likelihood ratio test in small samples. According to this correction, the mean value \(m\) of the likelihood ratio statistic under the null hypothesis is computed (here estimated by a parametric bootstrap) and the original LR statistic is multiplied by \(n/m\) where \(n\) is the number of degrees of freedom of the test.

If random effects are tested, only the raw p-value is computed. Its null distribution may include a probability mass in 1 (the discussion in Details of get_RLRsim_args applies).

* The ANOVA-table functionality has been included here mainly to provide access to F tests (including, for LMMs, the “Satterthwaite method” as developed by Fai and Cornelius, 1996), using pre-existing procedures as template or backend for expediency and familiarity:

ANOVA tables for LMs and GLMs are conceived to replicate the functionality, output format and details of base R anova, and therefore replicate some of their limitations, e.g., they only perform sequential analysis (“type 1”) in the same way as anova.lm and anova.glm. However, a difference occurs for Gamma GLMs, because the dispersion estimates for Gamma GLMs differ between stats::glm and spaMM fits (see Details in method). Therefore, F tests and Mallows' Cp differ too; results from spaMM REML fits being closer than ML fits to those from glm() fits;
For LMMs, ANOVA tables are provided by interfacing lmerTest::anova (with non-default type). This procedure should handle all types of LMMs that can be fitted by spaMM; yet, the displayed information should be inspected to check that some fitted random-effect parameters are not ignored when computing information for the Satterthwaite method.
For fitted models that do not lay within previous categories, such as GLMMs, models with a residual-dispersion submodel, and multivariate-response models, a table of tests for single-term deletions using the classical “Wald” chi-squared test based on coefficient values and their conditional standard error estimates will be returned. LRTs (moreover, with bootstrap correction) are more reliable than such tests and, as calling them requires a second model to be explicitly specified, they may also help users thinking about the hypothesis they are testing.

References

Bartlett, M. S. (1937) Properties of sufficiency and statistical tests. Proceedings of the Royal Society (London) A 160: 268-282.

Davison A.C., Hinkley D.V. (1997) Bootstrap methods and their applications. Cambridge Univ. Press, Cambridge, UK.

Fai AH, Cornelius PL (1996). Approximate F-tests of multiple degree of freedom hypotheses in generalised least squares analyses of unbalanced split-plot experiments. Journal of Statistical Computation and Simulation, 54(4), 363-378. tools:::Rd_expr_doi("10.1080/00949659608811740")

Examples

Run this code

data("wafers")
## Gamma GLMM with log link
m1 <- HLfit(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),family=Gamma(log),
          resid.model = ~ X3+I(X3^2) ,data=wafers,method="ML")
m2 <- update(m1,formula.= ~ . -I(X2^2))
#
anova(m1,m2)
try(anova(m1)) # fails because the 'resid.model' is not handled.

## ANOVA table for GLM
# Gamma example, from McCullagh & Nelder (1989, pp. 300-2), as in 'glm' doc:
clotting <- data.frame(
    u = c(5,10,15,20,30,40,60,80,100),
    lot1 = c(118,58,42,35,27,25,21,19,18),
    lot2 = c(69,35,26,21,18,16,13,12,12))
spglm <- fitme(lot1 ~ log(u), data = clotting, family = Gamma, method="REML")
anova(spglm, test = "F") 
anova(spglm, test = "Cp") 
anova(spglm, test = "Chisq")
anova(spglm, test = "Rao") 

## ANOVA table for LMM
if(requireNamespace("lmerTest", quietly=TRUE)) {
  lmmfit <- fitme(y ~X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch),data=wafers)
  print(anova(lmmfit)) # => Satterthwaite method, here giving p-values 
                       #   quite close to traditional t-tests given by:
  summary(lmmfit, details=list(p_value=TRUE))
}

## 'anova' (Wald chi-squared tests...) for GLMM 
wfit <- fitme(y ~ X1+X2+X1*X3+X2*X3+I(X2^2)+(1|batch), family=Gamma(log),
              rand.family=inverse.Gamma(log), resid.model = ~ X3+I(X3^2) , data=wafers)
anova(wfit)              

## Using resp_testfn argument for bootstrap LRT:
if (FALSE) {
set.seed(1L)
d <- data.frame(success = rbinom(10, size = 1, prob = 0.9), x = 1:10)
xx <- cbind(1,d$x)
table(d$success)
m_x <- fitme(success ~ x, data = d, family = binomial())
m_0 <- fitme(success ~ 1, data = d, family = binomial())
#
# Bootstrap LRTs:
anova(m_x, m_0, boot.repl = 100,
      resp_testfn=function(y) {! is_separated(xx,as.numeric(y),verbose=FALSE)})
}

## Models differing both in fixed and random effects:
if (spaMM.getOption("example_maxtime")>11) { 
 set.seed(123)
 dat <- data.frame(g = rep(1:10, e = 10), x = (x<-rnorm(100)), 
                   y = 0.1 * x + rnorm(100))
 m <- fitme(y ~ x + (1|g), data=dat)
 m0 <- fitme(y ~ 1, data=dat) 
 (bootpval <- LRT(m,m0, boot.repl = 199L)$rawBootLRT$p_value)
 ## See help("get_RLRsim_args") for a fast and accurate test procedure
}

Run the code above in your browser using DataLab