R/omnibusNULISAseq.R
omnibusNULISAseq.RdTests the overall significance of a model term that spans more than one
coefficient – most commonly an interaction involving a factor with more
than 2 levels, where a per-coefficient t-test doesn't answer "is this term
significant overall?" Works on the model fits already returned by
lmNULISAseq()/lmerNULISAseq() (return_model_fits = TRUE),
computing each requested term's joint test directly from the existing fit
– no refit needed, and no restriction to testing only one term per model.
omnibusNULISAseq(modelFits, terms = NULL, type = 3)Named list of fitted models (one per target) from
lmNULISAseq()/lmerNULISAseq() with
return_model_fits = TRUE. That argument defaults to FALSE
and omits modelFits entirely, so a NULL here means the
fitting call must be re-run with it enabled. Names must be present and
unique – results are keyed by them. All fits are assumed to share a model
formula and class. A non-fit entry is dropped with a per-target warning.
Character vector of term names as they appear in the model
formula (e.g. c("Group:Age", "Sex:Group")); component order need
not match. If NULL (default), every term spanning more than one
coefficient is tested – a single-coefficient term already has a direct
p-value in the coefficient table. Selection reads the design matrix
across all fits, not any one fit's realized degrees of freedom.
Type of sum of squares / hypothesis: 2 or 3
(default 3). Forwarded to car::Anova() for lm,
glm and glmerMod fits, and to lmerTest's
anova() (as "2"/"3") for linear mixed fits.
A data frame with one row per target per tested term:
Target name.
The model term tested.
Numerator degrees of freedom for the term.
Denominator degrees of freedom: the residual df for
lm fits, the term's Satterthwaite estimate for linear mixed
fits, NA on the chi-square paths (glm, glmerMod).
The test statistic – an F statistic for lm and
linear mixed fits, a likelihood-ratio chi-square for glm fits
and a Wald chi-square for glmerMod fits. The
mixed-model F (Satterthwaite) is NOT the likelihood-ratio test that
reduced_modelFormula_fixed's LRTstats produces – the two
are different statistics answering a similar question; do not treat
stat/pval here as directly comparable to LRTstats.
"F" or "Chisq", naming which stat is.
Unadjusted p-value for the term/target.
Benjamini-Hochberg adjusted p-value, computed within each
term across all targets (each term is its own independent family of
tests – consistent with how Fstats/LRTstats compute FDR).
Bonferroni adjusted p-value, same grouping as
pval_FDR.
Targets where the test failed (e.g. a singular fit, or a plain
lme4 fit that could not be converted for the Satterthwaite test)
are dropped with a warning, the same convention lmNULISAseq()/
lmerNULISAseq() use for their own per-target fit failures.
Aliased (rank-deficient) coefficients – routine in unbalanced designs
with empty interaction cells – do not fail the test: the term is tested
on its estimable coefficients only (singular.ok = TRUE on the
car::Anova() paths), with Df reduced to match. These are
the same coefficients the fit itself reports, so the omnibus stays
consistent with the coefficient table. lmer() drops aliased
columns at fit time, so the linear mixed path behaves this way already.
A term with NO estimable coefficient left in a fit is skipped for that
target with a warning.
The test statistic depends on the fit class:
lm fits: exact F test via car::Anova(), with the
residual df as denominator df.
Linear mixed fits (lmerModLmerTest, which is what
lmerNULISAseq() returns): F test with Satterthwaite-estimated
denominator df, via lmerTest's anova() method. The
finite denominator df matters: it reflects the number of independent
units that inform the term (e.g. participants for a
between-participant term, not observations), and a Wald chi-square –
which treats that df as infinite – overstates significance when it
is small. A plain lme4 fit is converted via
lmerTest::as_lmerModLmerTest() first; conversion needs the
fit's deviance function, so it typically fails on a deserialized fit,
and that target is then dropped with a warning like any other
per-target failure.
glm fits: likelihood-ratio chi-square via
car::Anova() – the same test glmNULISAseq_predict()
reports as Chisq_stat.
glmerMod fits: Wald chi-square via car::Anova() –
no Satterthwaite df exists for GLMMs, and a likelihood-ratio test
would need the reduced-model refit this function exists to avoid.
Note this differs from glmerNULISAseq_predict(), which does
refit and so reports an LR chi-square.
reduced_modelFormula/reduced_modelFormula_fixed in
lmNULISAseq()/lmerNULISAseq() answer a similar question by
refitting the model without a term and comparing with anova(). That
approach has two limitations this function does not have:
It tests only ONE term per extra fit – a model with several multi-coefficient terms needs one refit per term.
If the term is nested inside a still-present higher-order term
(e.g. dropping "A:B" while "A:B:C" remains), the
"reduced" model has the same rank as the full model – the comparison
is degenerate (zero residual-df difference, NA p-value),
silently, with no error.
Here each term's marginal contribution is computed directly from the single fit, unaffected by either issue, and every requested term is tested in one call.
if (FALSE) { # \dontrun{
res <- lmNULISAseq(data, sampleInfo, "sampleName",
modelFormula = "Sex * Group + Age",
return_model_fits = TRUE)
omnibusNULISAseq(res$modelFits, terms = c("Sex:Group"))
# Or auto-detect every multi-coefficient term:
omnibusNULISAseq(res$modelFits)
} # }