Skip to contents

Boundary-corrected likelihood-ratio test for variance components

Source: R/lrt-boundary.R
lrt-boundary.Rd

chibar_pvalue() and lrt_boundary() are the chi-bar-square (\(\bar\chi^2\)) boundary-corrected likelihood-ratio machinery for testing q variance components = 0, ported term-for-term from DRM.jl's src/chibar.jl (chibar_pvalue, lines 83-95; lrt_boundary, lines 134-140, at DRM.jl pin 430ef64cc). Testing a variance at zero is a boundary problem: the null value sits on the edge of the parameter space, so the usual \(\chi^2(q)\) reference for 2 * (logLik(full) - logLik(reduced)) is wrong and conservative (its p-values are too large, so the test loses power). Under the regularity conditions of Self and Liang (1987) and Stram and Lee (1994) the statistic follows a mixture of \(\chi^2\) distributions instead:

Usage

chibar_pvalue(statistic, q = 1L)

lrt_boundary(full, reduced, q = 1L)

Arguments

statistic

Likelihood-ratio statistic(s), 2 * (logLik(full) - logLik(reduced)). Vectorised; NA propagates.

q

Number of variance components tested at zero (the mixture order). Only 1 and 2 are supported; any other value aborts, as in DRM.jl.

full

A drmTMB fit carrying the q variance component(s).

reduced

The same model with those component(s) removed.

Value

chibar_pvalue() returns a numeric vector of upper-tail p-values the length of statistic.

lrt_boundary() returns a list of class drm_lrt_boundary with elements statistic (2 * (logLik(full) - logLik(reduced))), q, pvalue (chibar_pvalue(statistic, q)), pvalue_naive (the \(\chi^2(q)\) upper tail at max(statistic, 0), for comparison; always pvalue <= pvalue_naive), and df (df(full) - df(reduced), the number of extra parameters in full; DRM.jl does not report it). The two fits must be on the same observations; when df differs from q a drmTMB_lrt_boundary_df_mismatch warning says the mixture may not apply.

Details

  • q = 1 (one boundary parameter, e.g. dropping one (1 | g)): the null is \(0.5\,\chi^2_0 + 0.5\,\chi^2_1\), so p = 0.5 * P(chisq_1 > stat). At stat = 0 this is 0.5; for stat > 0 it is exactly half the naive \(\chi^2_1\) p-value.

  • q = 2 (two independent boundary parameters): the null is \(0.25\,\chi^2_0 + 0.5\,\chi^2_1 + 0.25\,\chi^2_2\), so p = 0.5 * P(chisq_1 > stat) + 0.25 * P(chisq_2 > stat). At stat = 0 this is 0.75.

drmTMB's anova() deliberately implements no likelihood-ratio comparison, and confint(method = "wald") only flags a variance component at its boundary (conf.status = "wald_at_boundary"); this pair is the first drmTMB tool that computes a p-value that is correct there.

Assumptions

  • The q parameters dropped between full and reduced are variances tested at the boundary 0, and reduced is full with exactly those components removed (nested models, same data).

  • For q = 2 the two components are independent (zero information correlation). With correlated components the 0.25/0.5/0.25 weights depend on the information matrix and this mixture is only approximate.

  • All other parameters are interior, so the standard asymptotic expansion holds for them.

  • Maximum-likelihood fits. REML log-likelihoods are comparable only across variance structures with the same fixed effects, so lrt_boundary() aborts when either fit is REML and the two mean structures differ, or when one fit is REML and the other ML. A penalized (penalty = drm_phylo_penalty(...), MAP) fit aborts too: its variance components are shrunk by the prior, so the likelihood ratio has no chi-square reference. Experimental MSPL fits abort as for every likelihood method.

A negative statistic (the reduced model fit better, which signals non-nested models or a fit that did not converge) is returned as-is but is clamped to 0 inside both p-values, which then take their boundary values (0.5 for q = 1, 0.75 for q = 2; 1 for the naive p-value). Inspect statistic directly in that case.

References

Self, S. G. and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82, 605-610.

Stram, D. O. and Lee, J. W. (1994). Variance components testing in the longitudinal mixed effects model. Biometrics, 50, 1171-1177.

Examples

chibar_pvalue(3.5, q = 1)
#> [1] 0.03068441
chibar_pvalue(3.5, q = 1) == 0.5 * pchisq(3.5, 1, lower.tail = FALSE)
#> [1] TRUE
chibar_pvalue(0, q = 1)   # 0.5, the boundary point mass
#> [1] 0.5
chibar_pvalue(0, q = 2)   # 0.75
#> [1] 0.75

# Random intercept vs no random effect: dropping (1 | g) removes ONE
# variance, so q = 1.
set.seed(20260610)
G <- 30; m <- 12
g <- factor(rep(seq_len(G), each = m))
x <- rnorm(G * m)
b <- 0.8 * rnorm(G)
dat <- data.frame(y = 0.5 - 0.4 * x + b[g] + 0.7 * rnorm(G * m), x = x, g = g)
full <- drmTMB(bf(y ~ x + (1 | g), sigma ~ 1), family = gaussian(), data = dat)
reduced <- drmTMB(bf(y ~ x, sigma ~ 1), family = gaussian(), data = dat)
lrt_boundary(full, reduced, q = 1)
#> Boundary-corrected likelihood-ratio test (chi-bar-square mixture)
#>   statistic = 142.4926 on 1 boundary variance component (df = 1)
#>   p-value (chi-bar-square) = < 2.2e-16
#>   p-value (naive chisq_1)   = < 2.2e-16