[parashardhapola]

Summary

Most single-cell studies are small, and variancePartition's dream tests them with
Kenward and Roger's degrees of freedom by default (ddf = "adaptive"), where
dream_voom used Satterthwaite's for every study. This adds Kenward and Roger's method
(the adjusted covariance of the fixed effects and its degrees of freedom, for one
contrast and for several) and F-tests of several contrasts at once, as pbkrtest and
lmerTest compute them; dream's choice of degrees of freedom, now dream_voom's
default, with dream's moderated F-tests; and the random effects' variables in
dream_voom's provenance. Everything is checked against pbkrtest 0.5.3, lmerTest 3.1-3
and variancePartition 1.36.2 on lme4 1.1-38.

Changes

• fixed_effect_tests(..., df_method="kenward_roger", kenward_roger_weights=...):
  t-tests with Kenward and Roger's adjusted standard errors and degrees of freedom
  (lmerTest's contest1D(ddf = "Kenward-Roger"), pbkrtest's vcovAdj and
  get_Lb_ddf). kenward_roger_weights says how the adjustment takes precision
  weights: "average" (default), every observation at its feature's average weight,
  as dream computes it with pbkrtest; "precision", each observation's own weight, the
  method for the weighted model. Maximum likelihood fits are refused for Kenward and
  Roger's method, as lmerTest refuses them.
• fixed_effect_f_test and its result FixedEffectFTest (new): the F-test that
  several contrasts are all zero, by Satterthwaite's method as lmerTest's contestMD
  or by Kenward and Roger's as pbkrtest's KRmodcomp. Dependent contrasts test the
  space they span; one contrast gives its t-test.
• FixedEffectTests gains kenward_roger_weights and covariance; MixedModelFit
  keeps its weights and component_levels, which Kenward and Roger's method needs.
• moderated_fixed_effect_tests moderates Kenward-Roger tests too, and
  ModeratedFixedEffectTests gains dream's moderated F-test: f_statistic, f_df1,
  f_df2, f_p_value and standardized_f_statistic (topTable's F.std).
• dream_voom(..., f_tests=None, df_method="adaptive", kenward_roger_weights="average").
  The default changes from Satterthwaite's degrees of freedom in every group to
  dream's choice in each group: Kenward and Roger's for fewer than 20 samples,
  Satterthwaite's otherwise, with dream's fallback (below). df_method="satterthwaite"
  gives the old results. f_tests={"time": ["t1_vs_t0", "t2_vs_t0"]} adds a table per
  F-test and group; a group that cannot estimate one of its contrasts lists the F-test
  as not estimable. The summary gains df_method; each group's record df_method,
  features_satterthwaite_fallback and f_test_df1; the design's record the random
  effects and the values of the variables they group by, as variance_partitioning's
  does (the helper moved to cytepool/_analysis_records.py).
• DifferentialExpressionResult.f_tests, and table(group, name) finds F-tests.
• cytepool/stats/_fixed_effect_inference.py (private): the compiled Kenward-Roger
  computation, the degrees of freedom of contrasts, the F-tests and lmerTest's rule for
  combining degrees of freedom, shared by the three modules.
• Documentation: sections on Kenward and Roger's method and on F-tests in "Mixed
  models, REML and dream", the API pages, a capability matrix row, a Validation section
  and summary row (31 scripts, 2,150 checks), the README, the Overview, a glossary
  entry, and three references with DOIs checked against Crossref (Kenward and Roger
  1997; Halekoh and Hojsgaard 2014; Fai and Cornelius 1996).

How dream behaves (found by running variancePartition 1.36.2)

• ddf = "adaptive" used Kenward and Roger's method for 19 samples and Satterthwaite's
  for 20, on designs where the two differ, although its help page says "<= 20". So the
  threshold is fewer than 20 samples, counted after a group's samples are filtered.
• With Kenward and Roger, dream returns lmerTest's contest1D(ddf = "Kenward-Roger")
  on lme4 fits whose weights are divided by their mean. Its standard errors are
  pbkrtest's adjusted ones, and eBayes adds the prior df to Kenward and Roger's.
• A gene whose adjusted covariance matrix is not positive definite or is ill
  conditioned falls back to Satterthwaite's tests (variancePartition's NEWS for
  1.31.16). On 450 genes built to straddle it, the bound is R's kappa(vcovAdj)
  estimate above 1e6. dream_voom uses the exact condition number, which decided all
  24 fallbacks of the validation's ill-conditioned cell type as dream did; the two can
  differ within about 2% of the bound. The rule depends on the units of continuous
  covariates, which the documentation says, with the remedy (centre and scale them).
• topTable's F with several coefficients is limma's t' R^-1 t / r from the moderated
  t-statistics and each gene's covariance (the adjusted one under Kenward and Roger),
  reproduced to 4e-8.

Kenward and Roger's method with weights

pbkrtest builds the covariance of the observations without lme4's prior weights, while
it starts from the weighted fit's covariance of the coefficients (found by running
pbkrtest 0.5.3). With dream's weights divided by their mean, this treats every
observation as having its gene's average weight. kenward_roger_weights="average" (the
default) reproduces it, so that dream_voom matches dream.

It is not Kenward and Roger's method for the weighted model. That method is
"precision", checked against pbkrtest on the equivalent unweighted model: the
values, fixed effects and random effects multiplied by the square roots of the weights.
With voom's weights the difference is not small. In the validation's weighted cases,
the average treatment gave standard errors 0.976 times the precision treatment's at the
median (0.51 to 1.05) and df 0.91 times (0.26 to 3.7). It also gave standard errors
below the unadjusted ones for 61% and 33% of the genes, which the precision treatment
never does. On the cube of validation/dream.py (cell types of 16 and 13 samples),
dream's treatment gave standard errors smaller than Satterthwaite's for 100% and 89% of
the genes, and the precision treatment larger ones for every gene.

How results change for small studies

Groups of fewer than 20 samples now get Kenward and Roger's tests. Compared with
Satterthwaite's tests, which dream_voom used before:

Data	Samples	p-values over before, median (5% to 95%)	Other
Simulated study	12	1.12 (0.97 to 1.66)	standard errors 1.06 times at the median
Simulated study	18	1.00 (0.91 to 1.07)
Validation cube, one cell type	12	0.99 (0.73 to 1.12)	12 genes below adjusted p 0.05, before and after
validation/dream.py cube, B	16	0.986 (0.81 to 1.01)	total df 15% lower, no gene crosses adjusted p 0.05
validation/dream.py cube, NK	13	0.974 (0.79 to 1.05)	total df 26% lower, no gene crosses adjusted p 0.05

Groups of 20 samples or more are unchanged. No documentation notebook calls
dream_voom, and the outputs of all four checked notebooks are current.

Deliberate differences from dream

• topTable computes every gene's F-test p-value on the total df of the first gene for
  the first coefficient listed (reproduced to 3e-14), so p-values change with the order
  of genes and coefficients. CytePool gives each gene its own: F is the mean of r
  squared uncorrelated moderated t-statistics, each with its contrast's df plus the
  prior df, combined by lmerTest's rule. An F-test of one contrast is then the moderated
  t-test. CytePool's p-values were 0.14 to 11.5 times dream's in the validation.
• F.std uses the number of coefficients where they are dependent; CytePool their rank.
• The fallback to Satterthwaite's tests is decided from the exact condition number.

Validation

validation/kenward_roger_and_f_tests.py (new): 143 checks, 115 pass, 16 information,
12 deliberate differences, none failed. validation/dream.py prints the same 63 checks
as before, now asking dream_voom for df_method="satterthwaite", which its R side
runs.

Comparison	Reference	Largest difference	Tolerance
Kenward and Roger's adjusted covariance, t-tests, df and p-values at CytePool's estimates (9 designs: time course, unbalanced, crossed, nested, boundary, and weights both ways with one and two random effects)	pbkrtest vcovAdj, lmerTest summary and contest1D	1.8e-12	1e-10
Kenward and Roger's F-tests (statistic, df, scaling, p-value)	pbkrtest KRmodcomp, lmerTest contestMD	8.2e-13	1e-10
Satterthwaite's F-statistic	lmerTest contestMD	8.2e-13	1e-10
Satterthwaite's F df and p-values (lmerTest differentiates numerically)	lmerTest contestMD	1.9e-7	1e-5
Kenward-Roger standard errors and df at lme4's own optimum	lmerTest summary	1.2e-7	1e-6
dream with each ddf on studies of 12, 18, 19, 20 and 24 samples: weights, standard errors, df, moderated t, B, z, p-values and F	voomWithDreamWeights, dream, eBayes, topTable	8.0e-8	1e-6
dream_voom with its default on cubes of 12, 15 (ill conditioned) and 24 samples	dream on each cell type	3.3e-7	1e-6
F-test p-values (deliberate)	topTable	0.14 to 11.5 times dream's
Genes falling back to Satterthwaite's tests (deliberate)	dream, R's kappa	24 of 24 decided alike

Reference: R 4.4.3, variancePartition 1.36.2, limma 3.62.1, lme4 1.1-38, lmerTest
3.1-3, pbkrtest 0.5.3, minqa 1.2.8.

RSCRIPT=$REF/r-dream/bin/Rscript PYTHONPATH=. python validation/kenward_roger_and_f_tests.py
RSCRIPT=$REF/r-dream/bin/Rscript PYTHONPATH=. python validation/dream.py

Tests

tests/test_kenward_roger_and_f_tests.py (new, 42 tests, 2 s with compiled code
cached):

• fixtures from pbkrtest and lmerTest for the 9 designs, and from variancePartition
  for dream's moderated tests and dream_voom's default;
• closed forms: in a balanced one-way design Kenward and Roger's df are q - 1 and the
  adjustment vanishes; in a randomised block design the df are the classical ones and
  the F-test of the treatments is the analysis of variance's, by both methods;
  lmerTest's rule for combining df on values derived by hand;
• invariants: a reparametrised design, equivalent contrasts, an F-test of one contrast
  equal to its t-test (also below 2 df, as contestMD gives it, and for dream's
  moderated F), relabelled donors and reordered samples, the scale of the weights,
  equal weights against none, and the number of threads;
• dream_voom: its choice per group, the F tables and records, the fallback, an F-test
  a group cannot estimate, and the error messages.

Updated: test_mixed_models.py, test_dream.py (its dream fixture is
ddf = "Satterthwaite", so it asks for that), test_masked_inputs.py and
test_private_names.py. Fixtures: tests/data/kenward_roger_and_f_tests/, 24 files,
each under 40 KB. The full suite: 4,141 passed, 18 skipped, also with bounds checks and
with the rounding changed.

Costs on one thread for 20,000 genes: Kenward and Roger's tests take 0.65 s at 20
samples and 13.5 s at 120 (Satterthwaite's 0.1 s).

For later

Random slopes (planned for 0.2); F-tests in limma_voom and limma_trend; Kenward and
Roger's method for maximum likelihood fits (refused, as lmerTest does); a faster
adjustment for large studies; sample quality weights for dream.