LRT performs a likelihood ratio (LR) test between two model fits. It differs from another function with the same effect, fixedLRT, by its arguments (model fits for LRT, but all arguments required to fit the models for fixedLRT), and by the format of its return value. LRT determines which model is the more complete one by comparing model components including the fixed-effect, random-effect, residual-dispersion model specifications, and response families (offsets are ignored). Then, a standard test based on the asymptotic chi-square distribution is performed. 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.
These different tests perform diffferently depending on the differences between the two models:
* If the models differ only by their fixed effects, the asymptotic LRT may be anticonservative, but the Bartlett-corrected one is generally well-calibrated.
* If the two models differ by their random effects, tests based on the chi-square distribution (including their Bartlett-corrected version) may be poorly behaved, as such tests assume unbounded parameters, as contrasted to, e.g., necessarily positive variances.
In such cases the raw boostrap test may be the only reliable test. The procedure aims to detect and report such issues, but may not report all problems: users remain responsible for applying the tests in valid conditions (see Caveats in Details section). In simple cases (such as comparing a fixed-effect to a mixed-effect model with the same fixed-effect term), the chi-square tests may not be reported. In other cases (see Examples) they may otherwise be reported, with a warning when the procedure detects some cases of estimates at the boundary for the full model, or detects cases where the LR statistic of bootstrap replicates is often zero (also suggesting that estimates are at the boundary in such replicates).
* If the fits differ by the fixed effects terms of their residual-dispersion models (but not by any random effect specification), tests based on the chi-square distribution are reported. A bootstrap can be performed as in other cases.
* Tests for some cases of nested response families (e.g., the Poisson versus its extensions) are tentatively allowed.
* In some cases the full and the null models cannot be identified and the basic LRT based on the chi-square distribution will not be returned, but a bootstrap test may still be performed.
* The case where residual-dispersion models of either fit include random effects is problematic as, for such fits, the fitting procedure does not maximize the reported likelihood. The basic LRT is not returned when the two fits differ by their random effects, but is still performed otherwise (see Examples); and a bootstrap test may still be performed in both cases.
LRT(object, object2, boot.repl = 0L, resp_testfn = NULL,
simuland = eval_replicate, include="call",
# 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 '...'
...)LRT returns an object of class fixedLRT, unless some exception occurs.
This is a list with typical elements (depending on the options)
the fit objects for the full and null models, or only their calls, depending on the include argument;
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:
A data frame including values of the likelihood ratio chi2 statistic, its degrees of freedom, and the raw bootstrap p-value;
A data frame including values of the Bartlett-corrected likelihood ratio chi2 statistic, its degrees of freedom, and its p-value;
a list with the following elements:
A table of fitted likelihoods for bootstrap replicates;
The mean likelihood ratio chi-square statistic for bootstrap replicates;
the state of .Random.seed at the beginning of the simulation of the bootstrap samples (present only is such simulation was run).
Fit object returned by a spaMM fitting function.
Optional second model fit to be compared to the first (their order does not matter, except in non-standard cases where the second model is taken as teh null one and a message is issued).
the number of bootstrap replicates.
See argument resp_testfn of spaMM_boot.
a function, passed to spaMM_boot. See argument eval_replicate for default value and requirements.
character string, either "fit" or "call" (default; other values are ignored).
This argument controls whether the fits compared, or only their calls, are included in the return value.
Further arguments, passed to spaMM_boot (e.g., for parallelization, or a random seed).
* Identifying a nested model: spaMM tries to guess which is the nested model, but this is not always possible.
In particular, the fit with (substantially) lower logLik is not necessarily the one with fewer degrees of freedom, if the models are not nested, or even for some form of “nestedness” with constrained parameters. The latter case may occur if one first fits a simple model, then fits a more complex model containing the first model but one constrains estimates of shared parameters to the values given by the first fit: the second fit should then have a higher likelihood but it may still have fewer estimated parameters than the first one. A likelihood ratio test is invalid anyway when such constraints are used.
Second, although nestedness is relatively easy to assess for fixed effects, equivalent random-effect models can be specified by diverse syntaxes, so a simple textual comparison of the random-effect terms may not be enough, and model specifications that hinder such a comparison should be avoided. When differences in random effects are tested, the null distribution of the LR may include a probability mass in 1: the discussion in Details of get_RLRsim_args applies.
* data checks: LRTs are invalid if the two models are fitted on different informative data. Informative data may differ even if the input data are identical, because the fitting procedures ignore lines of data with missing information. The retained data may then differ for the two models, if additional variables of the more complete model are missing from some lines of data retained in the fit of the less complete model.
LRT perform some checks based on comparison of the data used by the two models. However, a perfect check all cases where the informative data differ is not easily defined. For example, data may appear different but still contain the same information if lines of the data are reordered. Hence, users remain ultimately responsible for checking whether informative data are equivalent in the two models. They should keep in mind that data processing is applied on a per-submodel basis for multivariate-response fits, and pois4mlogit has further criteria for definition of informative data for each model.
* 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.
Bootstrap samples are not returned, and they are by default simulated internally (by a call to spaMM_boot itself calling simulate4boot). There are several ways to recover the samples if you need them: (1) simulate them beforehand and provide them as the boot_samples argument, passed to spaMM_boot through the ... argument; (2) provide a seed value (passed down to simulate4boot through the ... argument) and re-run simulate4boot with this seed and with other relevant arguments of the original LRT call; (3) same as (2), but using the set.seed function insted of the seed argument; (4) given only an LRT output where the previous solutions were ignored, you can still use the RNGstates value contained in the object returned by LRT to reset the generator state and re-run simulate4boot, as shown in the Examples for that function.
* Caveats: (1) An evaluated log-likelihood ratio can be slightly negative, e.g. when a fixed-effect model is compared to a mixed one, or a spatial random effect to a block effect, if parameters of the more complete model are estimated within bounds (e.g., variance>1e-06, or Matern smoothness>0.005) designed to avoid numerical singularities, while the less complete model corresponds to a boundary case (e.g., variance=0, or smoothness=0). The bootstrap procedure tries to identify these cases and then corrects slightly negative logL ratios to 0. (2) The Bartlett correction is applicable when the true distribution of the LRT departs smoothly from the chi-square distribution, but not in cases where it has a probability mass in zero (at typically occurs in the same boundary cases).
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.
See also fixedLRT and anova.HLfit
## 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)})
}
#### Various cases were asymptotic tests may be unreliable:
set.seed(123)
dat <- data.frame(g = rep(1:10, e = 10), x = (x<-rnorm(100)),
y = 0.1 * x + rnorm(100))
m0 <- fitme(y ~ 1, data=dat)
## (1) Models differing both by fixed and random effects:
#
# (note the warning for variance at boundary):
#
if (spaMM.getOption("example_maxtime")>11) {
m <- fitme(y ~ x + (1|g), data=dat)
LRT(m,m0, boot.repl = 199L)
}
## See help("get_RLRsim_args") for a fast and accurate test procedure
## (2) Models differing also by residual-dispersion models:
#
if (spaMM.getOption("example_maxtime")>25) {
m <- fitme(y ~ x + (1|g), data=dat, resid.model= ~x)
LRT(m,m0, boot.repl = 99L)
}
## (3) Models differing (also) by their random-effects in resid.model:
#
m <- fitme(y ~ x, data=dat, resid.model= ~1+(1|g))
LRT(m,m0) # no test performed
Run the code above in your browser using DataLab