diff --git a/README.md b/README.md index 07f1446..1b28eed 100644 --- a/README.md +++ b/README.md @@ -19,9 +19,9 @@ what is planned. | Normalization and filtering | TMM, upper-quartile and RLE factors, DESeq2 size factors, CPM, log-CPM, `filter_by_expression` | edgeR, DESeq2 | | Designs | Model matrices from R formulas, named contrasts, design checks | R, limma, lme4 | | Linear models | Ordinary, weighted and generalised least squares for every gene (in parallel when each gene needs its own fit), contrasts (named contrasts are matched to the coefficients by name) | limma | -| Mixed models | `fit_mixed_model`: linear mixed models with random intercepts of one or more grouping factors (nested or crossed), with observation weights, fitted to every gene by REML or maximum likelihood, with exact zeros for variance components at the boundary, fitted values with the predicted random effects, and approximate residual degrees of freedom; `fixed_effect_tests`: t-tests of contrasts with Satterthwaite's degrees of freedom, as lmerTest computes them for dream | lme4, lmerTest, variancePartition | +| Mixed models | `fit_mixed_model`: linear mixed models with random intercepts of one or more grouping factors (nested or crossed), with observation weights, fitted to every gene by REML or maximum likelihood, with exact zeros for variance components at the boundary, fitted values with the predicted random effects, and approximate residual degrees of freedom; `fixed_effect_tests`: t-tests of contrasts with Satterthwaite's degrees of freedom, as lmerTest computes them, or Kenward-Roger's adjusted standard errors and degrees of freedom, as pbkrtest computes them, with precision weights taken as dream takes them (every observation at its gene's average weight) or as the weighted model has them; `fixed_effect_f_test`: F-tests of several contrasts at once by either method (lmerTest's `contestMD`, pbkrtest's `KRmodcomp`) | lme4, lmerTest, pbkrtest, variancePartition | | Precision weights for counts | `voom`: log-CPM values and observation weights from the mean-variance trend, whose LOWESS span is chosen from the number of genes as in limma, with prior weights for samples, observations or genes, for a weighted linear model (limma-voom); `voom_with_quality_weights`, which combines them with sample quality weights | limma | -| dream | `voom_with_dream_weights`: voom's precision weights from the mean-variance trend of every gene's mixed model, fitted by maximum likelihood, as variancePartition's `voomWithDreamWeights`; `moderated_fixed_effect_tests`: Satterthwaite t-tests of a mixed model moderated by empirical Bayes as variancePartition's `eBayes` moderates `dream` fits, with `topTable`'s z-statistic | variancePartition | +| dream | `voom_with_dream_weights`: voom's precision weights from the mean-variance trend of every gene's mixed model, fitted by maximum likelihood, as variancePartition's `voomWithDreamWeights`; `moderated_fixed_effect_tests`: Satterthwaite or Kenward-Roger t-tests of a mixed model moderated by empirical Bayes as variancePartition's `eBayes` moderates `dream` fits, with `topTable`'s z-statistic, and the moderated F-test of several contrasts | variancePartition | | Sample quality weights | `sample_quality_weights`: a weight for each sample, or each group of samples, from how closely its values follow the linear model, by REML or the gene-by-gene update, as limma's `arrayWeights` | limma | | Negative binomial models | `fit_negative_binomial_model`: a negative binomial log-linear model for every gene, with known dispersions, as edgeR's `glmFit` (one-way layouts group by group, other designs by Newton's method with Levenberg damping, in parallel), with log fold changes shrunk by prior counts; exact unit deviances and deviances | edgeR, mpmath | | Quasi-likelihood F-tests | `fit_quasi_likelihood` and `quasi_likelihood_f_test`: edgeR 4's `glmQLFit` and `glmQLFTest` with their defaults: the common negative binomial dispersion of the most highly expressed genes by Cox-Reid adjusted profile likelihood (`common_dispersion`), deviances adjusted for small counts with exact moments of the unit deviance (`negative_binomial_unit_deviance_moments`), the average quasi-dispersion, limma's prior for unequal degrees of freedom with a trend in the average log CPM (`average_log_cpm`) or robust, prior weights and offsets, and F-tests of coefficients and named contrasts with edgeR's columns (logFC, logCPM, F, PValue, FDR) | edgeR, limma, mpmath | @@ -29,7 +29,7 @@ what is planned. | Empirical Bayes moderation | Moderated t and F statistics, B-statistics, TREAT, with constant, trend and robust priors, estimated as limma estimates them by default: by the method of moments, or, when the residual degrees of freedom differ between genes, by weighted moments and profile likelihood (limma's `fitFDistUnequalDF1`); either can be chosen | limma | | Gene-set tests | `camera` (competitive, allowing for the correlation between the genes of a set, given or estimated per set; parametric or rank-based) and `fry` (self-contained, directional and mixed p-values), for one contrast of a linear model, with voom or other precision weights and trend or robust moderation; the rotation tests `roast` and `mroast` (self-contained, with the mean, floor-mean, mean-50 and mean-square set statistics, contribution weights per gene or per set, and the proportions of active genes) and `romer` (competitive, by the ranks of the genes), with rotations from a seeded generator; gene sets named by gene id or position. Gene sets read from and written to GMT files (MSigDB, Reactome, WikiPathways, the Gene Ontology), with the source, version and licence of the collection, matched to an analysis and filtered by size; over-representation analysis of the genes up, down or either way in a differential-expression result, a table or a list of genes (limma's `kegga` and `goana`), with its dot plot | limma, R, decoupler | | Multiple testing | Bonferroni, Holm, Hochberg, Hommel, Benjamini-Hochberg, Benjamini-Yekutieli | R | -| Differential expression | `limma_voom` and `limma_trend`: limma-voom and limma-trend on every group (cell type) of a cube, and `dream_voom`: dream, with the design's random intercepts, such as a donor's, in every group's fits, with the design matched to the cube by sample id, samples with too few cells left out, rank-deficient groups fitted on their column space with the contrasts they can estimate, and a provenance record for every result. Options for pseudobulk data, off by default: sample quality weights, prior weights from the number of cells, a minimum total count per sample, and p-values adjusted across cell types | limma, edgeR, variancePartition; on simulated data and on published single-cell data | +| Differential expression | `limma_voom` and `limma_trend`: limma-voom and limma-trend on every group (cell type) of a cube, and `dream_voom`: dream, with the design's random intercepts, such as a donor's, in every group's fits, dream's choice of Kenward-Roger's or Satterthwaite's degrees of freedom by the number of samples in each group, and F-tests of several contrasts, with the design matched to the cube by sample id, samples with too few cells left out, rank-deficient groups fitted on their column space with the contrasts they can estimate, and a provenance record for every result. Options for pseudobulk data, off by default: sample quality weights, prior weights from the number of cells, a minimum total count per sample, and p-values adjusted across cell types | limma, edgeR, variancePartition; on simulated data and on published single-cell data | | Variance partitioning | `variance_fractions`: the fraction of every gene's variance that each variable of a design explains, as variancePartition's `fitExtractVarPartModel`: a mixed model of every gene by maximum likelihood (or REML), with the categorical variables as random intercepts, the continuous ones as fixed effects and precision weights such as voom's, or a linear model without random effects; `variance_partitioning`: the same on every cell type of a cube, with dream's precision weights, the samples each fit used and a provenance record | variancePartition, lme4 | | Meta-analysis | `meta_analysis`: for every gene, the fixed-effect (inverse-variance) or random-effects estimate of a common effect from the effects and standard errors of several studies, with the between-study variance by DerSimonian-Laird, REML or Paule-Mandel, Cochran's Q, I² and H², z-tests or the Knapp-Hartung adjustment, and genes measured in some studies only; `differential_expression_meta_analysis`: the same for the `limma_voom` or `limma_trend` results of several cubes, matched by cell type, contrast and gene id, with p-values adjusted over the genes and a provenance record | metafor | | Composition tests | `propeller`: differences in cell-type proportions between conditions, from the cell numbers of a cube or a table of them, with the logit or arcsine square-root transform, a design matched by sample id, t-tests of its contrasts and the F-test of a factor, robust empirical Bayes moderation, speckle's columns (baseline and mean proportions, proportion ratio, statistic, p-value, FDR) and a provenance record | speckle, limma | @@ -45,13 +45,10 @@ These are not implemented yet: group of a cube, and their `legacy = TRUE` methods - DESeq2's shrunken log fold changes (`lfcShrink` with apeglm, ashr or normal priors), its replacement of outlier counts, and its local and glmGamPoi dispersion trends -- The rest of dream (variancePartition): Kenward-Roger degrees of freedom (dream's - default for 20 samples or fewer), F-tests of several contrasts at once, and random - slopes, with which variancePartition also partitions the variance within subsets of - the samples. dream with random intercepts and Satterthwaite's - degrees of freedom is implemented (see the table above); `fit_linear_model` refuses - a design whose random effects it would drop, and its message points to - `fit_mixed_model`. +- Random slopes in mixed models, with which variancePartition also partitions the + variance within subsets of the samples. dream with random intercepts is implemented + (see the table above); `fit_linear_model` refuses a design whose random effects it + would drop, and its message points to `fit_mixed_model`. - Pathway or transcription-factor activity - Composition tests with precision weights or random effects (crumblr) - Ligand-receptor scoring and differential interactions @@ -76,7 +73,11 @@ These are not implemented yet: - A mixed model stores, for every gene, the covariance matrix of the coefficients and its derivative with respect to each variance parameter, for Satterthwaite's degrees of freedom: about 8 x genes x (components + 2) x coefficients² bytes, so 23 MB for - 20,000 genes, two random intercepts and 6 coefficients. `fit_mixed_model` fits random + 20,000 genes, two random intercepts and 6 coefficients, and its weights, as many + floats as the values. Kenward-Roger's tests invert each gene's n x n covariance + matrix of the observations, as pbkrtest does: O(n³) operations per gene for n + samples, 0.7 s for 20,000 genes of 20 samples and 13 s for 120 samples on one + thread, where Satterthwaite's take 0.1 s. `fit_mixed_model` fits random intercepts only, by REML or maximum likelihood, to complete data with positive weights: missing values or a zero weight make it raise an error, as no gene is fitted to fewer samples. Its fitted values take as much memory as the values. It also @@ -231,9 +232,16 @@ These are not implemented yet: limma's `contrasts.fit` refuses, and they report the contrasts that cannot be estimated, which limma either refuses or, when they weigh no aliased column, computes with a changed meaning. -- `dream_voom` and `moderated_fixed_effect_tests` use Satterthwaite's degrees of - freedom, as `dream(ddf = "Satterthwaite")`; dream's default, `ddf = "adaptive"`, - uses Kenward-Roger's for 20 samples or fewer, which are not available yet, nor are +- `dream_voom` chooses the degrees of freedom as `dream(ddf = "adaptive")` does: + Kenward-Roger's for fewer than 20 samples in a group (with 20, dream uses + Satterthwaite's, although its help page says "<= 20"), computed as dream computes + them with pbkrtest, with every observation at its gene's average weight, which can + make standard errors smaller than the unadjusted ones + (`kenward_roger_weights="precision"` uses each observation's weight), and + Satterthwaite's for a gene whose adjusted covariance + matrix is ill conditioned, which dream decides from R's estimate of the condition + number and CytePool from the exact one. dream's `topTable` gives the F-test of every + gene the degrees of freedom of the first gene; CytePool gives each gene its own. No sample quality weights for dream. As variancePartition, they moderate the residual variance of an observation of each gene's average weight, which keeps the mean-variance trend that voom's weights remove; `trend=True` lets the prior follow diff --git a/cytepool/_analysis_records.py b/cytepool/_analysis_records.py index 4f43992..5b2c35a 100644 --- a/cytepool/_analysis_records.py +++ b/cytepool/_analysis_records.py @@ -59,6 +59,18 @@ def design_record(design: Design, contrasts: Contrasts) -> dict[str, Any]: } +def mixed_design_record(design: Design, contrasts: Contrasts) -> dict[str, Any]: + """``design_record`` with the random effects, as ``"random"``, and the values of + the variables they group by for every sample, among the ``"samples"``.""" + record = design_record(design, contrasts) + record["random"] = [str(effect) for effect in design.random] + samples = design.samples + for effect in design.random: + for name in effect.grouping: + record["samples"].setdefault(name, samples[name].tolist()) + return record + + def software_versions() -> dict[str, str]: """The versions of CytePool, numpy, numba and Python that run an analysis.""" return { diff --git a/cytepool/differential_expression.py b/cytepool/differential_expression.py index 3613463..dcc2f4b 100644 --- a/cytepool/differential_expression.py +++ b/cytepool/differential_expression.py @@ -118,15 +118,16 @@ dream: :func:`dream_voom` takes the steps of limma-voom up to the normalization, then - those of variancePartition's ``voomWithDreamWeights``, ``dream(ddf = - "Satterthwaite")``, ``eBayes`` and ``topTable``: precision weights from the - mean-variance trend of the design's mixed model, fitted to every feature by - maximum likelihood (:func:`~cytepool.stats.dream.voom_with_dream_weights`); the - mixed model of every feature fitted by REML with those weights - (:func:`~cytepool.stats.mixed_models.fit_mixed_model`); Satterthwaite's (1946) - t-tests of the estimable contrasts - (:func:`~cytepool.stats.mixed_models.fixed_effect_tests`); - and their empirical Bayes moderation + those of variancePartition's ``voomWithDreamWeights``, ``dream``, ``eBayes`` and + ``topTable``: precision weights from the mean-variance trend of the design's + mixed model, fitted to every feature by maximum likelihood + (:func:`~cytepool.stats.dream.voom_with_dream_weights`); the mixed model of every + feature fitted by REML with those weights + (:func:`~cytepool.stats.mixed_models.fit_mixed_model`); t-tests of the estimable + contrasts (:func:`~cytepool.stats.mixed_models.fixed_effect_tests`) with the + degrees of freedom dream chooses by default (``ddf = "adaptive"``), Kenward and + Roger's (1997) in a group of fewer than 20 samples and Satterthwaite's (1946) + otherwise; F-tests of sets of contrasts; and their empirical Bayes moderation (:func:`~cytepool.stats.dream.moderated_fixed_effect_tests`). The random effects of a group are those of its samples: a grouping factor's levels without samples in the group drop out. :mod:`cytepool.stats.dream` documents the conventions found by @@ -233,6 +234,9 @@ Satterthwaite FE (1946). An approximate distribution of estimates of variance components. Biometrics Bulletin 2, 110-114. doi:10.2307/3002019 + Kenward MG, Roger JH (1997). Small sample inference for fixed effects from + restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558 + Heumos L, Schaar AC, Lance C, Litinetskaya A, Drost F, Zappia L, Luecken MD, Strobl DC, Henao J, Curion F, Single-cell Best Practices Consortium, Schiller HB, Theis FJ (2023). Best practices for single-cell analysis across modalities. Nature @@ -256,6 +260,7 @@ cube_record, design_record, matched_samples, + mixed_design_record, software_versions, ) from ._columns import text_column @@ -278,7 +283,7 @@ ) from ._jit import thread_limit from ._masked_arrays import reject_masked -from ._messages import PACKAGE_DIRECTORY, listing +from ._messages import PACKAGE_DIRECTORY, listing, quoted from ._options import ( checked_choice, checked_flag, @@ -291,13 +296,20 @@ from .design import Contrasts, Design, make_contrasts from .design._design_matrices import is_random_intercept, refuse_random_effects from .normalization import log_cpm +from .stats._fixed_effect_inference import ( + InferenceParts, + chosen_parts, + contrast_tests, + kenward_roger_parts, + satterthwaite_parts, +) from .stats._mean_variance_trend import scaled_prior_weights, trend_span, voom_inputs from .stats._mixed_model_designs import random_intercepts -from .stats.dream import _dream_weights, moderated_fixed_effect_tests +from .stats.dream import _dream_weights, _moderated_tests from .stats.empirical_bayes import moderated_statistics from .stats.linear_models import contrasts_fit, fit_linear_model +from .stats.mixed_models import FixedEffectTests, MixedModelFit from .stats.mixed_models import _fit as _fit_mixed_model -from .stats.mixed_models import fixed_effect_tests from .stats.multiple_testing import adjust_p_values from .stats.sample_weights import ( _estimate, @@ -325,13 +337,18 @@ "p_value_adjustment": "BH", } # dream's fits, recorded in the provenance: maximum likelihood for the weights, as -# variancePartition's voomWithDreamWeights fits, REML and Satterthwaite's degrees of -# freedom for the tests (dream's REML = TRUE, ddf = "Satterthwaite"). -_DREAM_SETTINGS = { - "weights_fit": "ml", - "fit": "reml", - "df_method": "satterthwaite", -} +# variancePartition's voomWithDreamWeights fits, and REML for the tests (dream's REML = +# TRUE). +_DREAM_SETTINGS = {"weights_fit": "ml", "fit": "reml"} +_DF_METHODS = ("adaptive", "satterthwaite", "kenward_roger") +# dream's ddf = "adaptive" uses Kenward and Roger's method for fewer samples than this. +# variancePartition 1.36.2 documents "<= 20 samples", but gives Kenward and Roger's +# degrees of freedom for 19 samples and Satterthwaite's for 20 (found by probing). +_KENWARD_ROGER_SAMPLES = 20 +# dream tests a feature by Satterthwaite's method when Kenward and Roger's adjusted +# covariance matrix of the coefficients is not positive definite or has a condition +# number above this (variancePartition's NEWS for 1.31.16; the bound found by probing). +_CONDITION_LIMIT = 1e6 _METHODS = ("limma_voom", "limma_trend", "dream_voom") @@ -360,10 +377,11 @@ class DifferentialExpressionResult: (``B``). Tables from :func:`dream_voom` also give ``df``, the degrees of freedom of each - test (Satterthwaite's plus the prior's, ``df.total``), after ``t``, and - ``z_statistic``, the z-statistic with the test's p-value (``topTable``'s - ``z.std``), after ``b_statistic``; a feature whose mixed model has no minimum has - NaN statistics. + test (Satterthwaite's or Kenward and Roger's plus the prior's, ``df.total``), + after ``t``, and ``z_statistic``, the z-statistic with the test's p-value + (``topTable``'s ``z.std``), after ``b_statistic``; a feature whose mixed model has + no minimum has NaN statistics. Its F-tests have tables of their own, whose + columns :func:`dream_voom` lists. The summary has one row per group: ``group_id``; ``status``, ``"fitted"`` or why the group was not fitted (``"no_samples"``, ``"no_residual_df"``, @@ -373,10 +391,13 @@ class DifferentialExpressionResult: ``n_samples_used`` and ``n_samples_left_out``; ``n_features_tested``; ``df_residual``, for dream the median over the features of their approximate residual degrees of freedom (in a group not fitted, the residual degrees of - freedom of the fixed effects, as for limma); ``prior_df``, the prior degrees of freedom (their - median over features with robust moderation, where each feature has its own); and - ``not_estimable``, the contrasts the group's samples cannot estimate, separated - by ``", "``. Numbers are NaN, or 0, where a group was not fitted that far. + freedom of the fixed effects, as for limma); ``prior_df``, the prior degrees of + freedom (their median over features with robust moderation, where each feature + has its own); for dream, ``df_method``, the degrees of freedom of the group's + tests (``"satterthwaite"`` or ``"kenward_roger"``; empty where the group was not + fitted); and ``not_estimable``, the contrasts (and F-tests) the group's samples + cannot estimate, separated by ``", "``. Numbers are NaN, or 0, where a group was + not fitted that far. The sample table has one row per group and design sample, in the cube's sample order: ``group_id``, ``sample_id``, ``n_cells`` (when the cube records cell @@ -509,9 +530,17 @@ def samples(self) -> Table: @property def tables(self) -> Mapping[tuple[str, str], Table]: - """The table of every group and contrast that was tested, by (group, contrast).""" + """The table of every group and contrast, or F-test, that was tested, by + (group, contrast) or (group, F-test).""" return self._tables + @property + def f_tests(self) -> dict[str, tuple[str, ...]]: + """The F-tests of :func:`dream_voom`, each with the contrasts it tests at + once, as a new dict; empty for the other methods or without F-tests.""" + tests = json.loads(self._record)["settings"].get("f_tests") or {} + return {name: tuple(members) for name, members in tests.items()} + @property def settings(self) -> dict[str, Any]: """The settings of the analysis, as a new dict (also in the provenance).""" @@ -525,17 +554,21 @@ def provenance(self) -> dict[str, Any]: It holds the method; the CytePool, numpy, numba and Python versions; the cube (``"cube"``): its fingerprint, the digests it is built from and its own provenance; the layer read; the design (formula, columns, levels, the values - of its variables for every sample, and the weights of every contrast); the - settings; and, for every group, its status, the samples that entered its - design, in the order they were fitted, those left out with the reason, the - design columns its samples cannot separate, the contrasts it cannot - estimate with the reason, the number of features tested, for limma-voom and - dream the span of its trend (``voom_span``, null where voom did not run), for - dream the number of features whose mixed model has no minimum - (``features_not_fitted``) and of those whose REML fit stopped without - converging (``features_not_converged``; their tests use the estimates where - it stopped), the - sample weights its fit used, in the order of its samples + of its variables for every sample, and the weights of every contrast; for + dream also its random effects, ``"random"``, and the values of the variables + they group by); the settings; and, for every group, its status, the samples + that entered its design, in the order they were fitted, those left out with + the reason, the design columns its samples cannot separate, the contrasts it + cannot estimate with the reason, the number of features tested, for + limma-voom and dream the span of its trend (``voom_span``, null where voom + did not run), for dream the number of features whose mixed model has no + minimum (``features_not_fitted``) and of those whose REML fit stopped + without converging (``features_not_converged``; their tests use the + estimates where it stopped), the degrees of freedom of its tests + (``df_method``), the number of features whose Kenward-Roger tests fell back + to Satterthwaite's (``features_satterthwaite_fallback``) and the numerator + degrees of freedom of each F-test (``f_test_df1``), the sample weights its + fit used, in the order of its samples (``prior_weights`` with ``prior_weights``, ``quality_weights`` with ``quality_weights``; null where the group was not fitted), and the warnings that the estimation of its sample quality weights gave (``warnings``, such as @@ -550,11 +583,12 @@ def provenance(self) -> dict[str, Any]: return record def table(self, group: str, contrast: str) -> Table: - """Return the table of one group and contrast. + """Return the table of one group and contrast, or F-test. Args: group: A group id. - contrast: A contrast name. + contrast: A contrast name, or the name of an F-test of + :func:`dream_voom`. Raises: KeyError: If the group or contrast was not analysed, or the group was @@ -569,10 +603,12 @@ def table(self, group: str, contrast: str) -> Table: f"no group {group!r} was analysed; groups are: " f"{listing(list(record)) or '(none)'}" ) - if contrast not in self._contrasts.names: + f_tests = list(self.f_tests) + if contrast not in self._contrasts.names and contrast not in f_tests: raise KeyError( f"no contrast {contrast!r} was tested; contrasts are: " f"{listing(list(self._contrasts.names))}" + + (f"; F-tests are: {listing(f_tests)}" if f_tests else "") ) entry = record[group] if entry["status"] != FITTED: @@ -888,12 +924,15 @@ def dream_voom( *, groups: str | Sequence[str] | None = None, contrasts: Mapping[str, Any] | Sequence[str] | str | None = None, + f_tests: Mapping[str, Sequence[str]] | None = None, min_cells: int = 10, min_counts: float = 0, expression_filter: bool | Mapping[str, float] = True, normalization: Literal["tmm"] | None = "tmm", span: float | None = None, prior_weights: Literal["n_cells"] | None = None, + df_method: Literal["adaptive", "satterthwaite", "kenward_roger"] = "adaptive", + kenward_roger_weights: Literal["average", "precision"] = "average", robust: bool = False, trend: bool = False, adjust_across_groups: bool = False, @@ -910,9 +949,9 @@ def dream_voom( :func:`~cytepool.stats.dream.voom_with_dream_weights`); every feature's mixed model is fitted by REML with those weights (``dream``, :func:`~cytepool.stats.mixed_models.fit_mixed_model`); its estimable contrasts - are tested with Satterthwaite's degrees of freedom - (:func:`~cytepool.stats.mixed_models.fixed_effect_tests`); and the tests are - moderated by empirical Bayes as variancePartition's ``eBayes`` moderates them + are tested (:func:`~cytepool.stats.mixed_models.fixed_effect_tests`) with the + degrees of freedom of ``df_method``; and the tests are moderated by empirical + Bayes as variancePartition's ``eBayes`` moderates them (:func:`~cytepool.stats.dream.moderated_fixed_effect_tests`). The module documentation gives the steps and the differences from variancePartition. @@ -923,9 +962,35 @@ def dream_voom( whose mixed model has no minimum, which the fixed and random effects reproduce, is tested with NaN statistics; its group's record counts them. - Satterthwaite's degrees of freedom are used for every group, as ``dream(..., - ddf = "Satterthwaite")``; dream's default, ``ddf = "adaptive"``, uses - Kenward-Roger's for 20 samples or fewer, which are not available yet. + Degrees of freedom: + As dream's ``ddf``, by default (``"adaptive"``) Kenward and Roger's method + in a group with fewer than 20 samples and Satterthwaite's otherwise. + variancePartition documents "20 samples or fewer"; with 20 its ``dream`` + uses Satterthwaite's (found by probing variancePartition 1.36.2). By + default (``kenward_roger_weights="average"``), Kenward and Roger's + adjustment is computed as dream computes it with pbkrtest, with every + observation at its feature's average weight; with voom's weights that + adjustment can make standard errors smaller than the unadjusted ones, + which ``"precision"``, the method for the weighted model, never does (see + :func:`~cytepool.stats.mixed_models.fixed_effect_tests`). As dream does, a + feature whose adjusted covariance matrix is not positive definite or has a + condition number above 1e6 is tested by Satterthwaite's method instead; + the group's record counts such features. That number depends on the units + of continuous covariates: centring and scaling them avoids it. + + F-tests: + ``f_tests`` names sets of contrasts to test at once, such as every + difference among three time points: each is tested by dream's moderated + F-statistic, which ``topTable`` gives for several coefficients (see + :mod:`cytepool.stats.dream`), in every group that can estimate all its + contrasts. Its table has, for every tested feature, ``feature_id`` (and + ``symbol``), ``log2_fold_change_`` for each contrast, + ``average_log2_cpm``, ``f_statistic``, ``df`` (the denominator degrees of + freedom), ``p_value``, ``adjusted_p_value`` (and + ``adjusted_p_value_across_groups``), and ``standardized_f_statistic``, the + F-statistic on infinite denominator degrees of freedom with the same + p-value (``F.std``). The numerator degrees of freedom are the number of + linearly independent contrasts, which the group's record gives. Args: cube: A cube with a ``counts`` layer. @@ -935,6 +1000,9 @@ def dream_voom( groups: A group id or a list of them; all groups by default. contrasts: The contrasts to test, matched to the design's columns by name; by default the design's named contrasts. + f_tests: A mapping of an F-test's name to the names of the contrasts it + tests at once, among those tested, such as ``{"time": ["t1_vs_t0", + "t2_vs_t0"]}``; None for none. A name must differ from every contrast's. min_cells: The fewest cells a sample needs in a group to enter its fit; ignored for cubes without cell numbers. min_counts: The fewest counts a sample needs in a group, over every @@ -949,6 +1017,11 @@ def dream_voom( prior_weights: ``"n_cells"`` to give each sample its number of cells in the group as a prior weight, as dreamlet weights samples, or None for none. Needs a cube that records cell numbers. + df_method: ``"adaptive"`` (dream's default, above), ``"satterthwaite"`` or + ``"kenward_roger"`` (with dream's fallback, above) for every group. + kenward_roger_weights: ``"average"``, as dream computes Kenward and Roger's + adjustment, or ``"precision"``, with each observation's own weight; + see above. robust: Whether to estimate the prior of the moderation robustly, as ``eBayes(robust = TRUE)`` does. trend: Whether the prior variance of the moderation follows the mean log2 @@ -963,20 +1036,24 @@ def dream_voom( keeps numba's setting. The results do not depend on it. Returns: - The tables of every group and estimable contrast, a summary of every group, - the samples used and left out, the settings and a provenance record, as for - :func:`limma_voom`; the tables also give each test's degrees of freedom - (``df``) and z-statistic (``z_statistic``). + The tables of every group and estimable contrast, and of every F-test, a + summary of every group, the samples used and left out, the settings and a + provenance record, as for :func:`limma_voom`; the tables of contrasts also + give each test's degrees of freedom (``df``) and z-statistic + (``z_statistic``), and the summary the degrees of freedom of each group's + tests (``df_method``). Raises: TypeError: If ``cube`` is not a Cube, ``design`` not a Design, an option has - the wrong type, or ``contrasts``, ``groups``, ``min_cells``, - ``expression_filter`` or ``n_threads`` is or holds a numpy masked array. + the wrong type, or ``contrasts``, ``f_tests``, ``groups``, + ``min_cells``, ``expression_filter`` or ``n_threads`` is or holds a + numpy masked array. ValueError: If the design has no random effects or one that is not a random intercept, has no ``sample_id`` column, names samples the cube does not have or names a sample twice, if there are no contrasts to - test, an option is invalid, or ``prior_weights`` is ``"n_cells"`` and - the cube does not record cell numbers. + test, an option is invalid, an F-test names a contrast that is not + tested or shares a contrast's name, or ``prior_weights`` is + ``"n_cells"`` and the cube does not record cell numbers. KeyError: If a group is not in the cube. """ return _analyse( @@ -992,6 +1069,11 @@ def dream_voom( "normalization": normalization_setting(normalization), "span": span_setting(span), "prior_weights": prior_weight_setting(prior_weights), + "df_method": checked_choice(df_method, "df_method", _DF_METHODS), + "kenward_roger_weights": checked_choice( + kenward_roger_weights, "kenward_roger_weights", ("average", "precision") + ), + "f_tests": _f_test_setting(f_tests), "robust": checked_flag(robust, "robust"), "trend": checked_flag(trend, "trend"), "adjust_across_groups": checked_flag( @@ -1012,6 +1094,9 @@ class _GroupFit(GroupRecord): prior_df: float = float("nan") quality_weights: np.ndarray | None = None tables: dict[str, Table] = field(default_factory=dict) + df_method: str = "" # dream: the degrees of freedom of the group's tests + features_satterthwaite_fallback: int = 0 # dream: Kenward-Roger tests not used + f_test_df1: dict[str, int] = field(default_factory=dict) # dream: per F-test def _analyse( @@ -1067,6 +1152,8 @@ def _analyse( data_ids="the cube's sample ids", ) tested = _contrasts_to_test(design, contrasts, method) + if settings.get("f_tests"): + _check_f_tests(settings["f_tests"], tested) positions = group_positions(cube, groups) inputs = GroupInputs(method, cube, order, cube_rows, tested, settings) fit_group = _fit_dream_group if method == "dream_voom" else _fit_group @@ -1081,12 +1168,13 @@ def _analyse( ) fits.append(fit) if settings["adjust_across_groups"]: - _adjust_across_groups(fits, tested.names) + _adjust_across_groups(fits, [*tested.names, *(settings.get("f_tests") or {})]) sample_ids = cube.samples["sample_id"][cube_rows] + dream = method == "dream_voom" return DifferentialExpressionResult( method, tested, - _summary_table(fits), + _summary_table(fits, dream), sample_table(cube, cube_rows, sample_ids, positions, fits, settings), { (fit.group_id, contrast): table @@ -1098,7 +1186,7 @@ def _analyse( "versions": software_versions(), "cube": cube_record(cube), "layer": COUNTS_LAYER, - "design": design_record(design, tested), + "design": (mixed_design_record if dream else design_record)(design, tested), "settings": { **settings, **_FIXED_SETTINGS, @@ -1247,25 +1335,28 @@ def _fit_dream_group(inputs: GroupInputs, design: Design, g: int) -> _GroupFit: return fit.stop("failed", f"voom with dream weights: {error}") model = _fit_mixed_model(weighted.log_cpm, weighted.weights, False, data.mixed) fit.voom_span = weighted.span + fitted = np.isfinite(model.residual_variance) # Differential expression always has contrasts to test. contrasts = cast(Contrasts, inputs.contrasts) names = tuple(contrasts.names[p] for p in data.estimable) - # An estimable contrast is its kept part applied to the fit on the kept columns. - tests = fixed_effect_tests( - model, - Contrasts( - contrasts.matrix[np.ix_(data.kept, data.estimable)], names, data.columns - ), + # An estimable contrast is its kept part applied to the fit on the kept columns, + # which are the model's coefficients. + matrix = contrasts.matrix[np.ix_(data.kept, data.estimable)] + use_kenward_roger = settings["df_method"] == "kenward_roger" or ( + settings["df_method"] == "adaptive" and n_used < _KENWARD_ROGER_SAMPLES ) + average = settings["kenward_roger_weights"] == "average" + with thread_limit(settings["n_threads"]): + parts = _dream_parts(model, use_kenward_roger, average, fitted, fit) + tests = _contrast_tests(model, matrix, names, parts, average) try: - moderated = moderated_fixed_effect_tests( - tests, robust=settings["robust"], trend=settings["trend"] + moderated = _moderated_tests( + tests, parts, robust=settings["robust"], trend=settings["trend"] ) except ValueError as error: # for example, no feature has a fit return fit.stop("failed", f"empirical Bayes moderation: {error}") adjusted = adjust_p_values(moderated.p_value, "bh", axis=0) - fitted = np.isfinite(model.residual_variance) fit.features_not_fitted = int(np.count_nonzero(~fitted)) fit.features_not_converged = int(np.count_nonzero(fitted & ~model.converged)) fit.df_residual = float(np.median(model.residual_df[fitted])) @@ -1291,9 +1382,103 @@ def _fit_dream_group(inputs: GroupInputs, design: Design, g: int) -> _GroupFit: "z_statistic": moderated.z[:, column], } ) + for f_name, members in (settings["f_tests"] or {}).items(): + missing = [member for member in members if member not in names] + if missing: + fit.not_estimable[f_name] = ( + "not tested: the samples of this group cannot estimate its " + f"contrast{'' if len(missing) == 1 else 's'} {quoted(missing)}" + ) + continue + columns = [names.index(member) for member in members] + subset = _contrast_tests( + model, matrix[:, columns], tuple(members), parts, average + ) + f_moderated = _moderated_tests( + subset, parts, robust=settings["robust"], trend=settings["trend"] + ) + fit.f_test_df1[f_name] = f_moderated.f_df1 + fit.tables[f_name] = Table( + { + **identity, + **{ + f"log2_fold_change_{member}": tests.estimate[:, column] + for member, column in zip(members, columns, strict=True) + }, + "average_log2_cpm": model.mean_expression, + "f_statistic": f_moderated.f_statistic, + "df": f_moderated.f_df2, + "p_value": f_moderated.f_p_value, + "adjusted_p_value": adjust_p_values(f_moderated.f_p_value, "bh"), + "standardized_f_statistic": f_moderated.standardized_f_statistic, + } + ) return fit +def _dream_parts( + model: MixedModelFit, + kenward_roger: bool, + average_weight: bool, + fitted: np.ndarray, + fit: _GroupFit, +) -> InferenceParts: + """The parts of the group's tests, with dream's choice of degrees of freedom. + + With Kenward and Roger's method, a feature whose adjusted covariance matrix is + not positive definite or has a condition number above 1e6 takes Satterthwaite's + parts, as dream does; ``fit`` records the method and how many features did. + Call it within ``thread_limit``. + """ + satterthwaite = satterthwaite_parts(model) + if not kenward_roger: + fit.df_method = "satterthwaite" + return satterthwaite + adjusted = kenward_roger_parts(model, average_weight) + usable = np.all(np.isfinite(adjusted.covariance), axis=(1, 2)) + condition = np.full(usable.size, np.inf) + if usable.any(): + values = np.linalg.eigvalsh(adjusted.covariance[usable]) + with np.errstate(divide="ignore", invalid="ignore"): + condition[usable] = np.where( + values[:, 0] > 0.0, values[:, -1] / values[:, 0], np.inf + ) + fallback = fitted & ~(condition <= _CONDITION_LIMIT) + fit.df_method = "kenward_roger" + fit.features_satterthwaite_fallback = int(np.count_nonzero(fallback)) + return chosen_parts(adjusted, satterthwaite, fallback) + + +def _contrast_tests( + model: MixedModelFit, + matrix: np.ndarray, + names: tuple[str, ...], + parts: InferenceParts, + average_weight: bool, +) -> FixedEffectTests: + """The t-tests of the contrasts in the columns of ``matrix``, from ``parts``.""" + estimate, standard_error, t, df, p_value = contrast_tests( + model.coefficients, parts, matrix + ) + covariance = np.einsum("jc,gjk,kd->gcd", matrix, parts.covariance, matrix) + kenward_roger = parts.method == "kenward_roger" + return FixedEffectTests( + fit=model, + contrasts=matrix, + contrast_names=names, + df_method=parts.method, + kenward_roger_weights=( + ("average" if average_weight else "precision") if kenward_roger else None + ), + estimate=estimate, + standard_error=standard_error, + t=t, + df=df, + p_value=p_value, + covariance=(covariance + np.swapaxes(covariance, 1, 2)) / 2.0, + ) + + def _voom_weights( y: np.ndarray, x: np.ndarray, @@ -1389,6 +1574,60 @@ def _contrasts_to_test(design: Design, contrasts: Any, function: str) -> Contras return chosen +def _f_test_setting(value: Any) -> dict[str, list[str]] | None: + """The F-tests: a mapping of name to the names of its contrasts, or None.""" + reject_masked(value, "f_tests") + if value is None: + return None + if not isinstance(value, Mapping): + raise TypeError( + "f_tests must be a mapping of an F-test's name to the names of the " + "contrasts it tests at once, such as {'time': ['t1_vs_t0', 't2_vs_t0']}, " + f"or None, got {type(value).__name__}" + ) + tests: dict[str, list[str]] = {} + for name, members in value.items(): + if not isinstance(name, str): + raise TypeError(f"f_tests must be keyed by names (str), got {name!r}") + reject_masked(members, f"f_tests[{name!r}]") + if isinstance(members, str) or not isinstance(members, Sequence | np.ndarray): + raise TypeError( + f"f_tests[{name!r}] must be a list of contrast names, got " + f"{type(members).__name__}" + ) + listed = list(members) + if not listed or not all(isinstance(member, str) for member in listed): + raise ValueError( + f"f_tests[{name!r}] must name at least one contrast, as text, got " + f"{listed!r}" + ) + repeated = sorted({member for member in listed if listed.count(member) > 1}) + if repeated: + raise ValueError( + f"f_tests[{name!r}] names the contrast(s) {quoted(repeated)} more than " + "once" + ) + tests[name] = listed + return tests or None + + +def _check_f_tests(f_tests: Mapping[str, Sequence[str]], tested: Contrasts) -> None: + """Refuse F-tests of contrasts that are not tested, or named as a contrast.""" + for name, members in f_tests.items(): + if name in tested.names: + raise ValueError( + f"the F-test {name!r} has the name of a contrast; tables are found by " + "name, so give it another" + ) + unknown = [member for member in members if member not in tested.names] + if unknown: + raise ValueError( + f"the F-test {name!r} names {quoted(unknown)}, which " + f"{'is' if len(unknown) == 1 else 'are'} not among the contrasts " + f"tested: {quoted(list(tested.names))}" + ) + + def _sample_weight_settings( prior_weights: Any, quality_weights: Any, method: str ) -> dict[str, Any]: @@ -1446,29 +1685,33 @@ def _quality_settings(value: Any, method: str) -> dict[str, Any] | None: # Results. -def _summary_table(fits: list[_GroupFit]) -> Table: - return Table( - { - "group_id": text_column([fit.group_id for fit in fits]), - "status": text_column([fit.status for fit in fits]), - "reason": text_column([fit.reason for fit in fits]), - "n_samples_used": np.array( - [np.count_nonzero(fit.used) for fit in fits], dtype=np.int64 - ), - "n_samples_left_out": np.array( - [np.count_nonzero(~fit.used) for fit in fits], dtype=np.int64 - ), - "n_features_tested": np.array( - [fit.tested.size if fit.status == FITTED else 0 for fit in fits], - dtype=np.int64, - ), - "df_residual": np.array([fit.df_residual for fit in fits], dtype=float), - "prior_df": np.array([fit.prior_df for fit in fits], dtype=float), - "not_estimable": text_column( - [", ".join(fit.not_estimable) for fit in fits] - ), - } +def _summary_table(fits: list[_GroupFit], dream: bool) -> Table: + columns: dict[str, Any] = { + "group_id": text_column([fit.group_id for fit in fits]), + "status": text_column([fit.status for fit in fits]), + "reason": text_column([fit.reason for fit in fits]), + "n_samples_used": np.array( + [np.count_nonzero(fit.used) for fit in fits], dtype=np.int64 + ), + "n_samples_left_out": np.array( + [np.count_nonzero(~fit.used) for fit in fits], dtype=np.int64 + ), + "n_features_tested": np.array( + [fit.tested.size if fit.status == FITTED else 0 for fit in fits], + dtype=np.int64, + ), + "df_residual": np.array([fit.df_residual for fit in fits], dtype=float), + "prior_df": np.array([fit.prior_df for fit in fits], dtype=float), + } + if dream: + # The degrees of freedom of the group's tests; empty where none were made. + columns["df_method"] = text_column( + [fit.df_method if fit.status == FITTED else "" for fit in fits] + ) + columns["not_estimable"] = text_column( + [", ".join(fit.not_estimable) for fit in fits] ) + return Table(columns) # The provenance record. @@ -1502,6 +1745,14 @@ def _group_record( record["features_not_converged"] = ( fit.features_not_converged if fitted else None ) + # The degrees of freedom of the tests, the features whose Kenward-Roger + # tests fell back to Satterthwaite's, and the numerator degrees of freedom + # of each F-test; null when the group was not fitted. + record["df_method"] = fit.df_method if fitted else None + record["features_satterthwaite_fallback"] = ( + fit.features_satterthwaite_fallback if fitted else None + ) + record["f_test_df1"] = dict(fit.f_test_df1) if fitted else None # The sample weights of the fit, in the order of "samples"; null when not fitted. for key in ("prior_weights", "quality_weights"): if settings.get(key) is not None: diff --git a/cytepool/stats/_fixed_effect_inference.py b/cytepool/stats/_fixed_effect_inference.py new file mode 100644 index 0000000..5f8af8e --- /dev/null +++ b/cytepool/stats/_fixed_effect_inference.py @@ -0,0 +1,748 @@ +"""Degrees of freedom and F-tests of the fixed effects of mixed model fits. + +The tests of :mod:`cytepool.stats.mixed_models`, dream's moderated tests +(:mod:`cytepool.stats.dream`) and ``dream_voom`` share: + +- ``InferenceParts``, made by ``satterthwaite_parts`` and ``kenward_roger_parts``: per + feature, the covariance matrix of the fixed effects that the statistics use, the + unadjusted one, its derivatives with respect to the variance parameters, and the + covariance matrix of those parameters, from which the degrees of freedom of any + contrast follow. ``chosen_parts`` takes each feature's parts from one of two sets. +- ``contrast_tests``: the t-tests of contrasts. +- ``f_test``: the F-test that several contrasts are all zero, by Satterthwaite's + method as lmerTest's ``contestMD`` computes it or by Kenward and Roger's as + pbkrtest's ``KRmodcomp`` does. +- ``combined_df``: the denominator degrees of freedom of an F-statistic that is the + mean of the squares of independent t-statistics, from theirs. +- ``contrast_rank``: the number of linearly independent contrasts. + +Kenward and Roger (1997): + The covariance matrix of the observations is linear in the variance parameters + phi = (sigma_1^2, ..., sigma_K^2, s^2), Sigma = sum_i phi_i G_i, with G_k = + Z_k Z_k' for the random effects and G = R, a diagonal matrix, for the residuals. + With Phi the covariance matrix of the fixed effects, P_i = -X' Sigma^-1 G_i + Sigma^-1 X, Q_ij = X' Sigma^-1 G_i Sigma^-1 G_j Sigma^-1 X and W the inverse of the + expected information of phi under REML, (1/2) tr(P G_i P G_j) with P = Sigma^-1 - + Sigma^-1 X Phi X' Sigma^-1, the adjusted covariance matrix is Phi + 2 Lambda with + + Lambda = Phi (sum_ij W_ij (Q_ij - P_i Phi P_j)) Phi + + (their R_ij vanish, the covariance being linear in phi). For a contrast c the + degrees of freedom are 2 (c' Phi c)^2 / sum_ij W_ij (c' E_i c)(c' E_j c) with E_i = + -Phi P_i Phi, those of Satterthwaite's formula with the expected information; for + l contrasts they come with a scaling of the F-statistic (Kenward and Roger 1997, + section 4; Halekoh and Hojsgaard 2014, section 3). + + pbkrtest builds Sigma for an lme4 fit without its prior weights, R = I with s^2 + the fit's residual variance, while Phi is the fit's own covariance matrix, which + the weights enter (found by probing pbkrtest 0.5.3). dream divides each feature's + weights by their mean before it fits, so that s^2 is the residual variance of an + observation of average weight: ``kenward_roger_parts`` computes this with + ``average_weight=True``. With ``average_weight=False``, R = diag(w_mean / w), and + Sigma is the covariance matrix of the weighted model, whose Phi it is: Kenward and + Roger's method for that model. The two agree when the weights are equal. + + Each feature takes O(n^3) operations for n samples, as in pbkrtest, whose + ``vcovAdj`` also inverts the n x n covariance matrix Sigma. The fixed effects + enter in an orthonormal basis of their columns, so that badly scaled covariates do + not cost digits, and the results are transformed back. + +References: + Fai AH-T, Cornelius PL (1996). Approximate F-tests of multiple degree of freedom + hypotheses in generalized least squares analyses of unbalanced split-plot + experiments. Journal of Statistical Computation and Simulation 54, 363-378. + doi:10.1080/00949659608811740 + + Halekoh U, Hojsgaard S (2014). A Kenward-Roger approximation and parametric + bootstrap methods for tests in linear mixed models: the R package pbkrtest. + Journal of Statistical Software 59(9), 1-32. doi:10.18637/jss.v059.i09 + + Kenward MG, Roger JH (1997). Small sample inference for fixed effects from + restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558 + + Kuznetsova A, Brockhoff PB, Christensen RHB (2017). lmerTest package: tests in + linear mixed effects models. Journal of Statistical Software 82(13), 1-26. + doi:10.18637/jss.v082.i13 +""" + +import math +from concurrent.futures import ThreadPoolExecutor +from typing import Any, NamedTuple + +import numba +import numpy as np +import numpy.typing as npt + +from .._jit import jit +from .._linear_algebra import cholesky, column_space, symmetric +from .distributions import f_sf, t_sf + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.int64] + +# lmerTest's get_Fstat_ddf: degrees of freedom within this of each other are averaged. +_EQUAL_DF = 1e-8 + + +class InferenceParts(NamedTuple): + """What the tests of contrasts of a mixed model fit compute their results from. + + Rows are features; the variance parameters are the variance components and the + residual variance, in that order. + """ + + covariance: FloatArray # (G, p, p) covariance of the fixed effects for the tests + unadjusted: FloatArray # (G, p, p) the fit's covariance matrix + derivatives: FloatArray # (G, K + 1, p, p) of the covariance, or -Phi P_i Phi + parameter_covariance: FloatArray # (G, K + 1, K + 1) of the variance parameters + method: str # "satterthwaite" or "kenward_roger" + + +class FTests(NamedTuple): + """F-tests that several contrasts are all zero, one per feature.""" + + statistic: FloatArray # (G,) F-statistics + df1: int # numerator degrees of freedom: the number of independent contrasts + df2: FloatArray # (G,) denominator degrees of freedom + p_value: FloatArray # (G,) + scaling: FloatArray # (G,) Kenward and Roger's scaling of the statistic (or 1) + + +def satterthwaite_parts(fit: Any) -> InferenceParts: + """The parts of Satterthwaite's tests: the fit's covariance matrix, its exact + derivatives and the inverse of the observed information of the variance + parameters, from which lmerTest computes them.""" + return InferenceParts( + fit.cov_coefficients, + fit.cov_coefficients, + fit.cov_coefficients_derivatives, + fit.variance_parameter_covariance, + "satterthwaite", + ) + + +def kenward_roger_parts(fit: Any, average_weight: bool) -> InferenceParts: + """The parts of Kenward and Roger's tests of a REML fit; see the module docstring. + + Args: + fit: A ``MixedModelFit`` by REML, with its weights and the levels of its + grouping factors. + average_weight: Whether the adjustment and the degrees of freedom take every + observation at its feature's average weight, as pbkrtest and dream do, + rather than at its own. + + Returns: + The adjusted covariance matrix of the fixed effects, the fit's, the matrices + -Phi P_i Phi and the inverse of the expected information, for every feature; + NaN for a feature that has no fit, or whose expected information is singular. + """ + x = np.asarray(fit.design, dtype=np.float64) + n_samples, p = x.shape + n_features = fit.coefficients.shape[0] + levels = np.ascontiguousarray(fit.component_levels, dtype=np.int64) + n_components = levels.shape[0] + n_levels = (levels.max(axis=1) + 1).astype(np.int64) + # The fixed effects in an orthonormal basis: x = q r. + q, r = np.linalg.qr(x) + inverse = np.linalg.solve(r, np.eye(p)) + covariance = np.ascontiguousarray( + np.einsum("ij,gjk,lk->gil", r, fit.cov_coefficients, r, optimize=True) + ) + weights = fit.weights + if weights is None or average_weight: + scales = np.ones((1, n_samples)) + elif np.ndim(weights) == 1: + scales = (np.mean(weights) / np.asarray(weights))[np.newaxis, :] + else: + scales = (np.mean(weights, axis=0) / np.asarray(weights)).T + scales = np.ascontiguousarray(scales, dtype=np.float64) + # The residual variance of an observation of average weight, the s^2 of both. + residual = np.ascontiguousarray( + fit.residual_variance / fit.mean_weights, dtype=np.float64 + ) + components = np.ascontiguousarray(fit.variance_components, dtype=np.float64) + parameters = n_components + 1 + adjusted = np.full((n_features, p, p), np.nan) + derivatives = np.full((n_features, parameters, p, p), np.nan) + parameter_covariance = np.full((n_features, parameters, parameters), np.nan) + arguments = ( + np.ascontiguousarray(q), + levels, + n_levels, + scales, + components, + residual, + covariance, + ) + outputs = (adjusted, derivatives, parameter_covariance) + threads = min(numba.get_num_threads(), n_features) + if threads > 1: + # Runs of consecutive features on Python threads; the compiled routine + # releases the GIL, and each feature is computed by the same operations. + n_runs = min(n_features, 4 * threads) + bounds = [ + (run * n_features // n_runs, (run + 1) * n_features // n_runs) + for run in range(n_runs) + ] + with ThreadPoolExecutor(max_workers=threads) as pool: + for _ in pool.map( + lambda run: _kenward_roger_features(*arguments, *run, *outputs), + bounds, + ): + pass + else: + _kenward_roger_features(*arguments, 0, n_features, *outputs) + # Back to the coefficients of x: beta = r^-1 beta_q. + adjusted = np.einsum("ij,gjk,lk->gil", inverse, adjusted, inverse, optimize=True) + derivatives = np.einsum( + "ij,gajk,lk->gail", inverse, derivatives, inverse, optimize=True + ) + return InferenceParts( + symmetric(adjusted), + fit.cov_coefficients, + symmetric(derivatives), + parameter_covariance, + "kenward_roger", + ) + + +def chosen_parts( + first: InferenceParts, second: InferenceParts, use_second: npt.NDArray[np.bool_] +) -> InferenceParts: + """Each feature's parts from ``second`` where ``use_second``, else ``first``, + under the method of ``first``.""" + pick = use_second.reshape(-1, 1, 1) + return InferenceParts( + np.where(pick, second.covariance, first.covariance), + np.where(pick, second.unadjusted, first.unadjusted), + np.where(pick[..., np.newaxis], second.derivatives, first.derivatives), + np.where(pick, second.parameter_covariance, first.parameter_covariance), + first.method, + ) + + +def contrast_tests( + coefficients: FloatArray, parts: InferenceParts, matrix: FloatArray +) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray, FloatArray]: + """t-tests of the contrasts in the columns of ``matrix``. + + Returns: + The estimates, standard errors, t-statistics, degrees of freedom and two-sided + p-values, ``(G, n_contrasts)`` each. + """ + estimate = coefficients @ matrix + variance = np.einsum("jc,gjk,kc->gc", matrix, parts.covariance, matrix) + df = contrast_df(parts, matrix) + with np.errstate(divide="ignore", invalid="ignore"): + standard_error = np.sqrt(variance) + t = estimate / standard_error + p_value = 2.0 * np.asarray(t_sf(np.abs(t), df), dtype=np.float64) + return estimate, standard_error, t, df, p_value + + +def contrast_df(parts: InferenceParts, matrix: FloatArray) -> FloatArray: + """The degrees of freedom of the contrasts in the columns of ``matrix``. + + 2 v^2 / (g' A g), with v = c' Phi c for the unadjusted Phi, g_i = c' D_i c and A + the covariance matrix of the variance parameters: Satterthwaite's, or Kenward and + Roger's for a single contrast. NaN where g' A g is not positive. + """ + variance = np.einsum("jc,gjk,kc->gc", matrix, parts.unadjusted, matrix) + gradient = np.einsum("jc,gajk,kc->gca", matrix, parts.derivatives, matrix) + spread = np.einsum( + "gca,gab,gcb->gc", gradient, parts.parameter_covariance, gradient + ) + with np.errstate(divide="ignore", invalid="ignore"): + df: FloatArray = np.where(spread > 0.0, 2.0 * variance**2 / spread, np.nan) + return df + + +def contrast_rank(matrix: FloatArray) -> int: + """The number of linearly independent contrasts among the columns of ``matrix``, + by the rank rule of R's ``qr()``.""" + return column_space(matrix).rank + + +def f_test( + coefficients: FloatArray, + parts: InferenceParts, + matrix: FloatArray, + rank: int, +) -> FTests: + """The F-test that the contrasts in the columns of ``matrix`` are all zero. + + The contrasts span ``rank`` dimensions. Each feature's are replaced by the + ``rank`` eigenvectors of their unadjusted covariance matrix with the largest + eigenvalues, c_m = L p_m, which span the same hypothesis; an exactly dependent + contrast adds an eigenvalue of zero. + + - Satterthwaite's method, as lmerTest's ``contestMD``: the estimates of c_m are + uncorrelated, with variances the eigenvalues d_m; F is the mean of their + squared t-statistics, and its denominator degrees of freedom combine their + Satterthwaite degrees of freedom by ``combined_df`` (Fai and Cornelius 1996; + Kuznetsova et al. 2017). + - Kenward and Roger's, as pbkrtest's ``KRmodcomp``: the Wald statistic with the + adjusted covariance matrix, divided by the rank and scaled by lambda, on + (rank, m) degrees of freedom (Kenward and Roger 1997, section 4). It does not + depend on how the hypothesis is spanned. + """ + covariance = np.einsum("jc,gjk,kd->gcd", matrix, parts.unadjusted, matrix) + values, vectors = _largest_eigenpairs(covariance, rank) + directions = np.einsum("jc,gcm->gjm", matrix, vectors) + estimates = np.einsum("gj,gjm->gm", coefficients, directions) + n_features = coefficients.shape[0] + if parts.method == "satterthwaite": + gradient = np.einsum( + "gjm,gajk,gkm->gma", directions, parts.derivatives, directions + ) + spread = np.einsum( + "gma,gab,gmb->gm", gradient, parts.parameter_covariance, gradient + ) + with np.errstate(divide="ignore", invalid="ignore"): + statistic = np.sum(estimates**2 / values, axis=1) / rank + df = np.where(spread > 0.0, 2.0 * values**2 / spread, np.nan) + df2 = combined_df(df) + scaling = np.ones(n_features) + else: + statistic, df2, scaling = _kenward_roger_f( + estimates, directions, values, parts, rank + ) + p_value = np.asarray(f_sf(statistic, float(rank), df2), dtype=np.float64) + return FTests(statistic, rank, df2, p_value, scaling) + + +def correlation_directions( + parts: InferenceParts, matrix: FloatArray, covariance: FloatArray, rank: int +) -> tuple[FloatArray, FloatArray, FloatArray]: + """The independent directions of t-statistics of contrasts, and their df. + + For the t-statistics of the contrasts in the columns of ``matrix``, whose + estimates have the covariance matrix ``covariance`` (``(G, c, c)``), the + eigenvectors v_m of their correlation matrix R with the ``rank`` largest + eigenvalues lambda_m give uncorrelated statistics v_m' t / lambda_m^(1/2), the + t-statistics of the contrasts L D^-1/2 v_m (D the variances of the estimates), + as limma's F-statistic t' R^-1 t / rank decomposes. Their degrees of freedom are + those of these contrasts by the method of ``parts``. + + Returns: + The eigenvalues ``(G, rank)``, the eigenvectors ``(G, c, rank)`` and the + degrees of freedom ``(G, rank)``, the largest eigenvalue first. + """ + with np.errstate(divide="ignore", invalid="ignore"): + scale = np.sqrt(np.diagonal(covariance, axis1=1, axis2=2)) + correlation = covariance / (scale[:, :, np.newaxis] * scale[:, np.newaxis, :]) + values, vectors = _largest_eigenpairs(correlation, rank) + with np.errstate(divide="ignore", invalid="ignore"): + directions = np.einsum( + "jc,gc,gcm->gjm", matrix, 1.0 / scale, vectors, optimize=True + ) + variance = np.einsum("gjm,gjk,gkm->gm", directions, parts.unadjusted, directions) + gradient = np.einsum("gjm,gajk,gkm->gma", directions, parts.derivatives, directions) + spread = np.einsum( + "gma,gab,gmb->gm", gradient, parts.parameter_covariance, gradient + ) + with np.errstate(divide="ignore", invalid="ignore"): + df = np.where(spread > 0.0, 2.0 * variance**2 / spread, np.nan) + return values, vectors, df + + +def _largest_eigenpairs( + matrices: FloatArray, rank: int +) -> tuple[FloatArray, FloatArray]: + """The ``rank`` largest eigenvalues of each symmetric matrix, the largest first as + R's ``eigen`` orders them, and their eigenvectors; NaN for a matrix that is not + finite, such as that of a feature without a fit.""" + values = np.full((matrices.shape[0], rank), np.nan) + vectors = np.full((matrices.shape[0], matrices.shape[1], rank), np.nan) + complete = np.all(np.isfinite(matrices), axis=(1, 2)) + if complete.any(): + found, basis = np.linalg.eigh(matrices[complete]) + values[complete] = found[:, ::-1][:, :rank] + vectors[complete] = basis[:, :, ::-1][:, :, :rank] + return values, vectors + + +def combined_df(df: FloatArray) -> FloatArray: + """The denominator degrees of freedom of the mean of m squared independent + t-statistics, from their degrees of freedom ``df`` (``(G, m)``), as lmerTest's + ``get_Fstat_ddf`` documents them. + + One keeps its own, so that the F-test of one contrast is its t-test, as lmerTest's + ``contestMD`` gives it also below 2 (found by running it). For several: 2 when any + is at most 2; their mean when consecutive ones differ by less than 1e-8; otherwise + 2 E / (E - m) with E = sum nu / (nu - 2), which matches the mean of the + F-statistic (Fai and Cornelius 1996). An infinite df counts 1 in E; all infinite + give infinity. NaN when any is NaN. + """ + m = df.shape[1] + if m == 1: + own: FloatArray = df[:, 0].copy() + return own + with np.errstate(divide="ignore", invalid="ignore"): + ratio = np.where(np.isinf(df), 1.0, df / (df - 2.0)) + expected = np.sum(ratio, axis=1) + combined = 2.0 * expected / (expected - m) + combined = np.where(expected == m, np.inf, combined) + with np.errstate(invalid="ignore"): # infinity minus infinity + steps = np.abs(np.diff(df, axis=1)) + equal = np.all(steps < _EQUAL_DF, axis=1) + combined = np.where(equal, np.mean(df, axis=1), combined) + combined = np.where(np.any(df <= 2.0, axis=1), 2.0, combined) + result: FloatArray = np.where(np.isnan(df).any(axis=1), np.nan, combined) + return result + + +def _kenward_roger_f( + estimates: FloatArray, + directions: FloatArray, + values: FloatArray, + parts: InferenceParts, + rank: int, +) -> tuple[FloatArray, FloatArray, FloatArray]: + """Kenward and Roger's F-statistic, denominator df and scaling, per feature. + + With Theta = L' (L Phi L')^-1 L for the unadjusted Phi, A1 = sum_ij W_ij + tr(Theta E_i) tr(Theta E_j) and A2 = sum_ij W_ij tr(Theta E_i Theta E_j), where the + rows of L are the directions, whose L Phi L' is diag(values). + """ + adjusted = np.einsum("gjm,gjk,gkn->gmn", directions, parts.covariance, directions) + terms = np.einsum("gjm,gajk,gkn->gamn", directions, parts.derivatives, directions) + with np.errstate(divide="ignore", invalid="ignore"): + scaled = terms / np.sqrt(values[:, np.newaxis, :, np.newaxis]) + scaled = scaled / np.sqrt(values[:, np.newaxis, np.newaxis, :]) + traces = np.einsum("gamm->ga", scaled) + pairs = np.einsum("gamn,gbnm->gab", scaled, scaled) + w = parts.parameter_covariance + a1 = np.einsum("gab,ga,gb->g", w, traces, traces) + a2 = np.einsum("gab,gab->g", w, pairs) + n = float(rank) + b = (a1 + 6.0 * a2) / (2.0 * n) + g = ((n + 1.0) * a1 - (n + 4.0) * a2) / ((n + 2.0) * a2) + denominator = 3.0 * n + 2.0 * (1.0 - g) + c1 = g / denominator + c2 = (n - g) / denominator + c3 = (n + 2.0 - g) / denominator + expected = 1.0 / (1.0 - a2 / n) + variance = (2.0 / n) * (1.0 + c1 * b) / ((1.0 - c2 * b) ** 2 * (1.0 - c3 * b)) + rho = variance / (2.0 * expected**2) + df2 = 4.0 + (n + 2.0) / (n * rho - 1.0) + scaling = df2 / (expected * (df2 - 2.0)) + statistic = scaling * _wald(estimates, adjusted) / n + return statistic, df2, scaling + + +def _wald(estimates: FloatArray, covariance: FloatArray) -> FloatArray: + """The Wald statistics e' C^-1 e of the features, NaN where C is singular or not + finite.""" + statistic = np.full(estimates.shape[0], np.nan) + complete = np.all(np.isfinite(covariance), axis=(1, 2)) & np.all( + np.isfinite(estimates), axis=1 + ) + rows = np.flatnonzero(complete) + try: + solved = np.linalg.solve(covariance[rows], estimates[rows][..., np.newaxis]) + statistic[rows] = np.sum(estimates[rows] * solved[..., 0], axis=1) + except np.linalg.LinAlgError: # one is singular: feature by feature + for row in rows: + try: + solved = np.linalg.solve(covariance[row], estimates[row]) + except np.linalg.LinAlgError: + continue + statistic[row] = float(estimates[row] @ solved) + return statistic + + +# The adjustment, compiled. + + +@jit(nogil=True) +def _kenward_roger_features( + basis: FloatArray, + codes: IntArray, + n_levels: IntArray, + residual_scales: FloatArray, + components: FloatArray, + residual_variance: FloatArray, + covariance: FloatArray, + start: int, + stop: int, + adjusted: FloatArray, + derivatives: FloatArray, + parameter_covariance: FloatArray, +) -> None: + """Kenward and Roger's parts of features ``start:stop``, in the basis. + + Everything is computed with Sigma / s^2 = R + sum_k gamma_k Z_k Z_k' (gamma_k = + sigma_k^2 / s^2), whose inverse is M, and U = M q: then X' Sigma^-1 G_i Sigma^-1 X + = U' G_i U / s^4 =: H_i / s^4, the matrices of Q_ij are U' G_i M G_j U / s^6, and + the expected information is (1/2) (tr(M G_i M G_j) - 2 tr(Phi_V U' G_i M G_j U) + + tr(Phi_V H_i Phi_V H_j)) / s^4 with Phi_V = Phi / s^2; Lambda is s^2 Phi_V + (sum_ij W_V,ij (U' G_i M G_j U - H_i Phi_V H_j)) Phi_V for W_V the inverse of the + information without its factor 1 / s^4. Only the rows of these features are + written, so runs of features can be computed on several threads at once. + + Args: + basis: ``(n, p)`` orthonormal basis of the fixed effects. + codes: ``(K, n)`` level of each sample for every grouping factor. + n_levels: The number of levels of each grouping factor. + residual_scales: ``(G, n)`` the diagonal of R, or one row for every feature. + components: ``(G, K)`` variance components. + residual_variance: ``(G,)`` s^2. + covariance: ``(G, p, p)`` Phi in the basis. + start, stop: The run of features. + adjusted, derivatives, parameter_covariance: The outputs, NaN until written: + Phi + 2 Lambda, -Phi P_i Phi = Phi_V H_i Phi_V, and W. + """ + n, p = basis.shape + n_components = codes.shape[0] + m = n_components + 1 + largest = 1 + for k in range(n_components): + largest = max(largest, n_levels[k]) + v = np.empty((n, n)) + lower_inverse = np.empty((n, n)) + inverse = np.empty((n, n)) + u = np.empty((n, p)) + sums = np.empty((largest, p)) + b = np.empty((m, n, p)) + c = np.empty((m, n, p)) + h = np.empty((m, p, p)) + scaled = np.empty((m, p, p)) + products = np.empty((m, m, p, p)) + level_columns = np.empty((max(n_components, 1), n, largest)) + blocks = np.empty((largest, largest)) + traces = np.empty((m, m)) + information = np.empty((m, m)) + information_inverse = np.empty((m, m)) + elimination = np.empty((m, 2 * m)) + phi = np.empty((p, p)) + inner = np.empty((p, p)) + half = np.empty((p, p)) + shared = residual_scales.shape[0] == 1 + for g in range(start, stop): + s2 = residual_variance[g] + if not (s2 > 0.0 and math.isfinite(s2)): + continue + usable = True + for k in range(n_components): + if not (components[g, k] >= 0.0 and math.isfinite(components[g, k])): + usable = False + for a in range(p): + for e in range(p): + if not math.isfinite(covariance[g, a, e]): + usable = False + if not usable: + continue + r = residual_scales[0] if shared else residual_scales[g] + # V = R + sum_k gamma_k Z_k Z_k', its Cholesky factor and its inverse M. + for i in range(n): + for j in range(i + 1): + value = r[i] if i == j else 0.0 + for k in range(n_components): + if codes[k, i] == codes[k, j]: + value += components[g, k] / s2 + v[i, j] = value + if not cholesky(v, n, 0.0): + continue + for column in range(n): + for row in range(column): + lower_inverse[row, column] = 0.0 + lower_inverse[column, column] = 1.0 / v[column, column] + for row in range(column + 1, n): + total = 0.0 + for t in range(column, row): + total += v[row, t] * lower_inverse[t, column] + lower_inverse[row, column] = -total / v[row, row] + for i in range(n): + for j in range(i + 1): + total = 0.0 + for t in range(i, n): + total += lower_inverse[t, i] * lower_inverse[t, j] + inverse[i, j] = total + inverse[j, i] = total + for i in range(n): + for a in range(p): + total = 0.0 + for t in range(n): + total += inverse[i, t] * basis[t, a] + u[i, a] = total + for a in range(p): + for e in range(p): + phi[a, e] = covariance[g, a, e] / s2 + # B_i = G_i U and H_i = U' G_i U: level sums for the random effects. + for k in range(n_components): + for level in range(n_levels[k]): + for a in range(p): + sums[level, a] = 0.0 + for t in range(n): + for a in range(p): + sums[codes[k, t], a] += u[t, a] + for t in range(n): + for a in range(p): + b[k, t, a] = sums[codes[k, t], a] + for a in range(p): + for e in range(a + 1): + total = 0.0 + for level in range(n_levels[k]): + total += sums[level, a] * sums[level, e] + h[k, a, e] = total + h[k, e, a] = total + for t in range(n): + for a in range(p): + b[n_components, t, a] = r[t] * u[t, a] + for a in range(p): + for e in range(a + 1): + total = 0.0 + for t in range(n): + total += r[t] * u[t, a] * u[t, e] + h[n_components, a, e] = total + h[n_components, e, a] = total + # C_i = M B_i, and the products B_i' C_j = U' G_i M G_j U. + for i in range(m): + for s in range(n): + for a in range(p): + total = 0.0 + for t in range(n): + total += inverse[s, t] * b[i, t, a] + c[i, s, a] = total + for i in range(m): + for j in range(m): + for a in range(p): + for e in range(p): + total = 0.0 + for t in range(n): + total += b[i, t, a] * c[j, t, e] + products[i, j, a, e] = total + # tr(M G_i M G_j), from M Z_k, the column sums of M by level. + for k in range(n_components): + for s in range(n): + for level in range(n_levels[k]): + level_columns[k, s, level] = 0.0 + for t in range(n): + level_columns[k, s, codes[k, t]] += inverse[s, t] + for k in range(n_components): + for other in range(k + 1): + for level in range(n_levels[k]): + for level_other in range(n_levels[other]): + blocks[level, level_other] = 0.0 + for s in range(n): + for level_other in range(n_levels[other]): + blocks[codes[k, s], level_other] += level_columns[ + other, s, level_other + ] + total = 0.0 + for level in range(n_levels[k]): + for level_other in range(n_levels[other]): + total += blocks[level, level_other] ** 2 + traces[k, other] = total + traces[other, k] = total + total = 0.0 + for s in range(n): + squares = 0.0 + for level in range(n_levels[k]): + squares += level_columns[k, s, level] ** 2 + total += r[s] * squares + traces[k, n_components] = total + traces[n_components, k] = total + total = 0.0 + for s in range(n): + for t in range(n): + total += r[s] * r[t] * inverse[s, t] ** 2 + traces[n_components, n_components] = total + # The expected information, without its factor 1 / s^4. + for i in range(m): + for a in range(p): + for e in range(p): + total = 0.0 + for f in range(p): + total += phi[a, f] * h[i, f, e] + scaled[i, a, e] = total + for i in range(m): + for j in range(i + 1): + cross = 0.0 + square = 0.0 + for a in range(p): + for e in range(p): + cross += phi[a, e] * products[i, j, e, a] + square += scaled[i, a, e] * scaled[j, e, a] + value = 0.5 * (traces[i, j] - 2.0 * cross + square) + information[i, j] = value + information[j, i] = value + if not _invert(information, m, information_inverse, elimination): + continue + # Lambda = Phi_V (sum_ij W_ij (B_i' C_j - H_i Phi_V H_j)) Phi_V, times s^2. + for a in range(p): + for e in range(p): + total = 0.0 + for i in range(m): + for j in range(m): + correction = 0.0 + for f in range(p): + correction += h[i, a, f] * scaled[j, f, e] + total += information_inverse[i, j] * ( + products[i, j, a, e] - correction + ) + inner[a, e] = total + for a in range(p): + for e in range(p): + total = 0.0 + for f in range(p): + total += phi[a, f] * inner[f, e] + half[a, e] = total + for a in range(p): + for e in range(p): + total = 0.0 + for f in range(p): + total += half[a, f] * phi[f, e] + adjusted[g, a, e] = s2 * phi[a, e] + 2.0 * s2 * total + for i in range(m): + for a in range(p): + for e in range(p): + total = 0.0 + for f in range(p): + total += scaled[i, a, f] * phi[f, e] + derivatives[g, i, a, e] = total + for i in range(m): + for j in range(m): + parameter_covariance[g, i, j] = s2 * s2 * information_inverse[i, j] + + +@jit(nogil=True) +def _invert( + matrix: FloatArray, size: int, inverse: FloatArray, work: FloatArray +) -> bool: + """The inverse of the small matrix ``matrix[:size, :size]`` into ``inverse``. + + Gauss-Jordan elimination with partial pivoting, which serves the symmetric + matrices that need not be positive definite. Returns False if a pivot is zero or + not finite. + """ + for i in range(size): + for j in range(size): + work[i, j] = matrix[i, j] + work[i, size + j] = 1.0 if i == j else 0.0 + for column in range(size): + pivot = column + for row in range(column + 1, size): + if abs(work[row, column]) > abs(work[pivot, column]): + pivot = row + value = work[pivot, column] + if not (value != 0.0 and math.isfinite(value)): + return False + if pivot != column: + for j in range(2 * size): + swap = work[column, j] + work[column, j] = work[pivot, j] + work[pivot, j] = swap + for j in range(2 * size): + work[column, j] /= value + for row in range(size): + if row != column: + factor = work[row, column] + if factor != 0.0: + for j in range(2 * size): + work[row, j] -= factor * work[column, j] + for i in range(size): + for j in range(size): + inverse[i, j] = work[i, size + j] + return True diff --git a/cytepool/stats/dream.py b/cytepool/stats/dream.py index 2dc9e1b..f067949 100644 --- a/cytepool/stats/dream.py +++ b/cytepool/stats/dream.py @@ -4,10 +4,10 @@ the linear model, for designs with repeated measures: several samples of a donor, tissues of one patient, time courses. The donor (or another grouping factor) gets a random intercept, whose variance is estimated for every gene, and contrasts of the -fixed effects are tested with Satterthwaite's (1946) degrees of freedom. dreamlet -(Hoffman et al. 2026) runs it on the pseudobulk profiles of every cell type. The -steps, as variancePartition runs them (``voomWithDreamWeights``, ``dream``, -``eBayes`` and ``topTable``): +fixed effects are tested with Kenward and Roger's (1997) degrees of freedom in small +studies and Satterthwaite's (1946) in larger ones. dreamlet (Hoffman et al. 2026) runs +it on the pseudobulk profiles of every cell type. The steps, as variancePartition runs +them (``voomWithDreamWeights``, ``dream``, ``eBayes`` and ``topTable``): 1. Precision weights, :func:`voom_with_dream_weights`. voom's (Law et al. 2014): the log2 counts per million are fitted, the square roots of the residual standard @@ -16,15 +16,20 @@ (:mod:`cytepool.stats.voom`). The fit is the mixed model of every feature, by maximum likelihood, and the fitted values include the predicted random effects. 2. Fits and tests: :func:`~cytepool.stats.mixed_models.fit_mixed_model` by REML with - those weights, and :func:`~cytepool.stats.mixed_models.fixed_effect_tests`, whose - t-statistics have Satterthwaite's degrees of freedom. + those weights, and :func:`~cytepool.stats.mixed_models.fixed_effect_tests`, with + the degrees of freedom of dream's ``ddf``: by default (``"adaptive"``) Kenward and + Roger's for fewer than 20 samples, with their adjusted standard errors, and + Satterthwaite's otherwise. 3. Moderation, :func:`moderated_fixed_effect_tests`. The residual variance of every feature, that of an observation of its average weight, is shrunk towards a prior by limma's empirical Bayes (Smyth 2004) on the approximate residual degrees of freedom of the mixed model, tr((I - H)' (I - H)) for its hat matrix H (``MixedModelFit.residual_df``). The moderated t-statistic divides an estimate by its standard error with the posterior variance in place of the residual - variance, and has Satterthwaite's degrees of freedom plus the prior's. + variance, and has the degrees of freedom of the test plus the prior's. Several + contrasts are tested at once by the moderated F-statistic t' R^-1 t / r, limma's, + for their moderated t-statistics and the correlation matrix R of their estimates, + with r the number of linearly independent contrasts. Conventions that the documentation leaves open were found by running variancePartition 1.36.2, with limma 3.62.1 and lme4 1.1-38, as a black box: @@ -62,8 +67,24 @@ the two-sided p-value: the z-statistic with the p-value of the t-statistic, so that tests on different degrees of freedom can be compared. It computes it from the p-value floored at 1e-300. - -Deliberate differences from variancePartition, reported by ``validation/dream.py``: +- ``dream(ddf = "adaptive")``, the default, uses Kenward and Roger's method for 19 + samples and Satterthwaite's for 20; the help page says "<= 20 samples". With + ``ddf = "Kenward-Roger"`` it tests with lmerTest's ``contest``, which calls + pbkrtest: Kenward and Roger's adjustment is computed with every observation at the + average weight (see :mod:`cytepool.stats.mixed_models`), the standard errors are + the adjusted ones, and ``eBayes`` moderates the tests as it does Satterthwaite's. A gene whose adjusted + covariance matrix of the coefficients is not positive definite or has a condition + number above 1e6, by R's estimate ``kappa``, is tested by Satterthwaite's method + instead (variancePartition's NEWS for 1.31.16 states the rule; the bound was found + by probing). The condition number depends on the units of continuous covariates. +- ``topTable`` with several coefficients gives the coefficients, ``AveExpr``, ``F``, + ``P.Value``, ``adj.P.Val`` and ``F.std``: F is limma's moderated F-statistic, from + the moderated t-statistics and each gene's covariance matrix of the coefficients, + Kenward and Roger's adjusted one under their method, and ``F.std`` the F-statistic + on infinite denominator degrees of freedom with the same p-value. + +Deliberate differences from variancePartition, reported by ``validation/dream.py`` +and ``validation/kenward_roger_and_f_tests.py``: - The LOWESS span is chosen from the number of features by default, as :func:`~cytepool.stats.voom.voom` and limma's ``voom`` (since limma 3.66) choose it. @@ -75,10 +96,22 @@ it is left out of the trend, its weights are read at its values, which the effects reproduce, and its tests are NaN. lme4 fits it with a residual standard deviation near zero, which enters the trend, and ``dream`` tests it. - -Kenward-Roger's degrees of freedom (``dream``'s default up to 20 samples), random -slopes and F-tests of several contrasts are not implemented yet; variancePartition's -variance partitioning is :mod:`cytepool.variance_partitioning`. +- ``topTable`` computes the p-values of an F-test of every gene on the total degrees + of freedom of the first gene for the first coefficient listed (found by probing), + so that they change with the order of the genes and of the coefficients. Here each + gene has its own: F is the mean of the squares of r uncorrelated moderated + t-statistics, those of the eigenvectors of R, whose degrees of freedom are those of + their contrasts plus the prior's, and the denominator degrees of freedom combine + them as lmerTest combines those of its F-tests (2 E / (E - r), with E = sum nu_m / + (nu_m - 2); see :mod:`cytepool.stats.mixed_models`). With one contrast, the F-test + is the moderated t-test. ``F.std`` uses the number of coefficients where it + differs from their rank; here the rank, the F-test's numerator degrees of freedom. +- ``dream_voom`` decides which genes fall back to Satterthwaite's method from the + exact condition number, where dream uses R's estimate: they agree except within + about 2% of 1e6. + +Random slopes are not implemented yet; variancePartition's variance partitioning is +:mod:`cytepool.variance_partitioning`. References: Hoffman GE, Roussos P (2021). Dream: powerful differential expression analysis for @@ -95,6 +128,9 @@ analysis tools for RNA-seq read counts. Genome Biology 15, R29. doi:10.1186/gb-2014-15-2-r29 + Kenward MG, Roger JH (1997). Small sample inference for fixed effects from + restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558 + Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027 @@ -126,12 +162,27 @@ voom_log_cpm, voom_mean_log_count, ) +from cytepool.stats._fixed_effect_inference import ( + InferenceParts, + combined_df, + contrast_rank, + correlation_directions, + kenward_roger_parts, + satterthwaite_parts, +) from cytepool.stats._mixed_model_designs import ( MixedDesign, maximum_likelihood, mixed_design, ) -from cytepool.stats.distributions import norm_isf, norm_logsf, t_logsf, t_sf +from cytepool.stats.distributions import ( + f_isf, + f_sf, + norm_isf, + norm_logsf, + t_logsf, + t_sf, +) from cytepool.stats.empirical_bayes import moderated_statistics from cytepool.stats.linear_models import LinearModelFit from cytepool.stats.mixed_models import FixedEffectTests, _fit @@ -187,6 +238,19 @@ class ModeratedFixedEffectTests: proportion: The assumed proportion of non-zero contrasts behind ``lods``. prior_estimator: The estimator of the prior: ``"moments"`` or ``"likelihood"``. + f_statistic: ``(n_features,)`` moderated F-statistics testing every + contrast of ``tests`` at once: t' R^-1 t / ``f_df1``, for the moderated + t-statistics t and the correlation matrix R of their estimates + (``topTable``'s ``F`` for these contrasts). + f_df1: The numerator degrees of freedom of the F-tests, the number of + linearly independent contrasts. + f_df2: ``(n_features,)`` their denominator degrees of freedom; see the + module docstring. + f_p_value: ``(n_features,)`` p-values of ``f_statistic`` on ``f_df1`` and + ``f_df2`` degrees of freedom. + standardized_f_statistic: ``(n_features,)`` the F-statistic on ``f_df1`` + and infinite degrees of freedom with the same p-value (``topTable``'s + ``F.std``), so that tests on different degrees of freedom compare. """ tests: FixedEffectTests @@ -203,6 +267,11 @@ class ModeratedFixedEffectTests: coefficient_prior_variance: FloatArray proportion: float prior_estimator: str + f_statistic: FloatArray + f_df1: int + f_df2: FloatArray + f_p_value: FloatArray + standardized_f_statistic: FloatArray def __repr__(self) -> str: return ( @@ -312,12 +381,15 @@ def moderated_fixed_effect_tests( shrunk towards a prior on the approximate residual degrees of freedom by limma's empirical Bayes (:func:`~cytepool.stats.empirical_bayes.moderated_statistics`, with its estimators and options), each contrast's t-statistic uses the posterior - variance, and its degrees of freedom are Satterthwaite's plus the prior's. See - the module docstring. + variance, and its degrees of freedom are those of the tests (Satterthwaite's or + Kenward and Roger's) plus the prior's. The moderated F-statistic tests every + contrast of ``tests`` at once, as ``topTable`` does for several coefficients. + See the module docstring. Args: tests: Tests from :func:`~cytepool.stats.mixed_models.fixed_effect_tests`, - of a REML fit as ``dream`` makes them by default. + of a REML fit as ``dream`` makes them by default. To F-test some of the + contrasts, test those alone. proportion: Assumed proportion of non-zero contrasts, in (0, 1). stdev_coef_lim: Limits for the standard deviation of a non-zero contrast, as for ``moderated_statistics``. @@ -342,6 +414,38 @@ def moderated_fixed_effect_tests( f"tests must be a FixedEffectTests, from fixed_effect_tests, got " f"{type(tests).__name__}" ) + if tests.df_method == "kenward_roger": + parts = kenward_roger_parts(tests.fit, tests.kenward_roger_weights == "average") + else: + parts = satterthwaite_parts(tests.fit) + return _moderated_tests( + tests, + parts, + proportion=proportion, + stdev_coef_lim=stdev_coef_lim, + trend=trend, + robust=robust, + winsor_tail_p=winsor_tail_p, + prior_estimator=prior_estimator, + ) + + +def _moderated_tests( + tests: FixedEffectTests, + parts: InferenceParts, + *, + proportion: float = 0.01, + stdev_coef_lim: Sequence[float] = (0.1, 4.0), + trend: bool | npt.ArrayLike = False, + robust: bool = False, + winsor_tail_p: float | Sequence[float] = (0.05, 0.1), + prior_estimator: Literal["moments", "likelihood"] | None = None, +) -> ModeratedFixedEffectTests: + """The steps of :func:`moderated_fixed_effect_tests`, with the parts of the tests. + + ``parts`` are those from which ``tests`` were computed, feature by feature: the + degrees of freedom of the directions of the F-test come from them. + """ fit = tests.fit matrix = tests.contrasts if matrix.shape[1] == 0: @@ -362,10 +466,9 @@ def moderated_fixed_effect_tests( sigma = np.sqrt(variance) with np.errstate(divide="ignore", invalid="ignore"): stdev_unscaled = tests.standard_error / sigma[:, np.newaxis] - per_feature = ( - np.einsum("jc,gjk,kd->gcd", matrix, fit.cov_coefficients, matrix) - / variance[:, np.newaxis, np.newaxis] - ) + # The covariance matrix that the tests use: dream's cov.coefficients.list, + # Kenward and Roger's adjusted one under their method (found by probing). + per_feature = tests.covariance / variance[:, np.newaxis, np.newaxis] cross = fit.design.T @ fit.design linear = LinearModelFit( coefficients=tests.estimate, @@ -393,6 +496,12 @@ def moderated_fixed_effect_tests( df_total = tests.df + prior_df[:, np.newaxis] t = moderated.t p_value = 2.0 * np.asarray(t_sf(np.abs(t), df_total), dtype=np.float64) + f_value, f_df1, f_df2 = _moderated_f(t, tests, parts, prior_df) + with np.errstate(invalid="ignore"): + f_p_value = np.asarray(f_sf(f_value, float(f_df1), f_df2), dtype=np.float64) + standardized = np.asarray( + f_isf(f_p_value, float(f_df1), np.inf), dtype=np.float64 + ) return ModeratedFixedEffectTests( tests=tests, residual_variance_at_mean_weight=_frozen(variance), @@ -408,9 +517,43 @@ def moderated_fixed_effect_tests( coefficient_prior_variance=moderated.coefficient_prior_variance, proportion=moderated.proportion, prior_estimator=moderated.prior_estimator, + f_statistic=_frozen(f_value), + f_df1=f_df1, + f_df2=_frozen(f_df2), + f_p_value=_frozen(f_p_value), + standardized_f_statistic=_frozen(standardized), ) +def _moderated_f( + t: FloatArray, tests: FixedEffectTests, parts: InferenceParts, prior_df: FloatArray +) -> tuple[FloatArray, int, FloatArray]: + """dream's moderated F-statistic of every contrast, and its degrees of freedom. + + F = t' R^-1 t / r, for the moderated t-statistics and the correlation matrix R of + the estimates, is the mean of the squares of r uncorrelated moderated + t-statistics, those of the eigenvectors of R (``correlation_directions``). Each + has the degrees of freedom of its contrast plus the prior's, as a moderated + t-statistic of dream does, and the denominator degrees of freedom combine them + as lmerTest combines those of an F-test (``combined_df``). + + Returns: + The F-statistics, the numerator and the denominator degrees of freedom. + """ + rank = contrast_rank(tests.contrasts) + n_features = t.shape[0] + if rank == 0: + return np.full(n_features, np.nan), 0, np.full(n_features, np.nan) + values, vectors, df = correlation_directions( + parts, tests.contrasts, tests.covariance, rank + ) + with np.errstate(divide="ignore", invalid="ignore"): + scores = np.einsum("gc,gcm->gm", t, vectors) / np.sqrt(values) + f_value: FloatArray = np.sum(scores**2, axis=1) / rank + f_df2 = combined_df(df + prior_df[:, np.newaxis]) + return f_value, rank, f_df2 + + def _dream_weights( y: FloatArray, sizes: FloatArray, diff --git a/cytepool/stats/mixed_models.py b/cytepool/stats/mixed_models.py index 64c8861..acfb62d 100644 --- a/cytepool/stats/mixed_models.py +++ b/cytepool/stats/mixed_models.py @@ -12,12 +12,14 @@ ``(1 | donor) + (1 | batch)``. The variance parameters are estimated by restricted maximum likelihood (REML; Patterson and Thompson 1971), or by maximum likelihood (``method="ml"``), and the fixed effects by generalised least squares given them, as -lme4's ``lmer`` with ``REML = TRUE`` (or ``FALSE``) does (Bates et al. 2015). Tests of -contrasts of the fixed effects use Satterthwaite's degrees of freedom, as lmerTest -computes them (Kuznetsova et al. 2017), which variancePartition's ``dream`` uses for -repeated measures in gene expression (Hoffman and Roussos 2021); -:mod:`cytepool.stats.dream` builds dream's precision weights and moderated tests on -these fits. +lme4's ``lmer`` with ``REML = TRUE`` (or ``FALSE``) does (Bates et al. 2015). Contrasts +of the fixed effects are tested one at a time by t-tests (:func:`fixed_effect_tests`), +or several at once by F-tests (:func:`fixed_effect_f_test`), with Satterthwaite's +degrees of freedom, as lmerTest computes them (Kuznetsova et al. 2017), or Kenward and +Roger's, as pbkrtest computes them (Halekoh and Hojsgaard 2014). variancePartition's +``dream`` uses both for repeated measures in gene expression (Hoffman and Roussos +2021); :mod:`cytepool.stats.dream` builds dream's precision weights and moderated +tests on these fits. The REML criterion: With theta_k = sigma_k / sigma, the relative standard deviations that lme4 @@ -108,9 +110,58 @@ by probing lmerTest). lmerTest computes the derivatives numerically; here they are exact. +Kenward and Roger's method: + With few samples the variance parameters are estimated imprecisely, and the + covariance matrix C of the estimates, computed as if they were known, + understates the variance of the fixed effects. Kenward and Roger (1997) add to C + a correction for its own bias and for the variability of the variance parameters, + of order 1 / n, to second order in their estimation error, whose covariance they + take as the inverse of the expected information of REML; the t-statistic of a + contrast uses the adjusted matrix, and its degrees of freedom are Satterthwaite's + formula with that expected information in place of the observed one. For an + F-test of l contrasts they scale the Wald statistic and match its first two + moments to an F distribution on (l, m) degrees of freedom. The covariance matrix + of the observations is linear in the variance parameters, so the second + derivatives that their formula also allows for vanish. A variance component + estimated at zero takes part, as in pbkrtest, so that the degrees of freedom are + smaller than Satterthwaite's, which leave it out. In balanced designs, such as a + randomised block design, the adjustment vanishes and the method gives the tests + of the analysis of variance, which are exact there. + + With weights, pbkrtest builds the covariance matrix of the observations without + the prior weights of an lme4 fit, with the residual variance of an observation + of weight one, while it takes C from the fit, which the weights enter (found by + probing pbkrtest 0.5.3); lmerTest's Kenward-Roger tests are pbkrtest's. dream + divides each feature's weights by their mean before it fits, so that its tests + are those of every observation at the average weight, which do not depend on the + scale of the weights. ``kenward_roger_weights="average"`` computes these. + ``"precision"`` uses each observation's own weight, the covariance + matrix of the model fitted: Kenward and Roger's method for it, which pbkrtest + gives for the same model fitted without weights to the values, fixed effects and + random effects multiplied by the square roots of the weights. With voom's + weights the two can differ by tens of percent in the degrees of freedom and by + several percent in the standard errors, and the first can make standard errors + smaller than the unadjusted ones, which the second never does (the Validation + page gives sizes). Each feature costs O(n^3) operations for n samples, as in + pbkrtest. + +F-tests: + :func:`fixed_effect_f_test` tests that the contrasts L beta are all zero. With + Satterthwaite's method, as lmerTest's ``contestMD``, the contrasts are replaced by + the eigenvectors of their covariance matrix L C L' with positive eigenvalues + (in the order of R's ``eigen``, the largest first), whose estimates are + uncorrelated; F is the mean of their squared t-statistics, and its denominator + degrees of freedom are 2 E / (E - l) with E = sum nu_m / (nu_m - 2) over their + Satterthwaite degrees of freedom nu_m, which matches the mean of F (Fai and + Cornelius 1996): their mean when they differ by less than 1e-8, and 2 when one is + at most 2, as lmerTest's ``get_Fstat_ddf`` documents. These depend on how the + contrasts span the hypothesis, as lmerTest's do. Kenward and Roger's F-test, as + pbkrtest's ``KRmodcomp``, does not. The F-test of one contrast is its t-test, on + the same degrees of freedom even below 2, as ``contestMD`` gives it (found by + running it). + Limits of this version: - Random slopes, Kenward-Roger degrees of freedom (dream's default for 20 samples - or fewer) and F-tests of several contrasts at once are not implemented yet. + Random slopes are not implemented yet. :mod:`cytepool.variance_partitioning` partitions the variance of these fits. Every feature needs a value for every sample, with a positive weight. Designs whose fixed and random effects together span every sample, and so reproduce @@ -125,6 +176,11 @@ constraints. SIAM Journal on Control and Optimization 20, 221-246. doi:10.1137/0320018 + Fai AH-T, Cornelius PL (1996). Approximate F-tests of multiple degree of freedom + hypotheses in generalized least squares analyses of unbalanced split-plot + experiments. Journal of Statistical Computation and Simulation 54, 363-378. + doi:10.1080/00949659608811740 + Giesbrecht FG, Burns JC (1985). Two-stage analysis based on a mixed model: large-sample asymptotic theory and small-sample simulation results. Biometrics 41, 477-486. doi:10.2307/2530872 @@ -132,6 +188,10 @@ Golub GH, Van Loan CF (2013). Matrix Computations, 4th edition. Johns Hopkins University Press, Baltimore. + Halekoh U, Hojsgaard S (2014). A Kenward-Roger approximation and parametric + bootstrap methods for tests in linear mixed models: the R package pbkrtest. + Journal of Statistical Software 59(9), 1-32. doi:10.18637/jss.v059.i09 + Harville DA (1977). Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association 72, 320-338. doi:10.1080/01621459.1977.10480998 @@ -143,6 +203,9 @@ repeated measures designs. Bioinformatics 37, 192-201. doi:10.1093/bioinformatics/btaa687 + Kenward MG, Roger JH (1997). Small sample inference for fixed effects from + restricted maximum likelihood. Biometrics 53, 983-997. doi:10.2307/2533558 + Kuznetsova A, Brockhoff PB, Christensen RHB (2017). lmerTest package: tests in linear mixed effects models. Journal of Statistical Software 82(13), 1-26. doi:10.18637/jss.v082.i13 @@ -168,17 +231,32 @@ from cytepool._jit import jit, thread_limit from cytepool._linear_algebra import RANK_TOLERANCE, cholesky, symmetric from cytepool._masked_arrays import real_array, reject_masked, unmasked_array +from cytepool._options import checked_choice from cytepool._summation import observed_means from cytepool.design._contrast_weights import contrast_matrix from cytepool.design.design import Design +from cytepool.stats._fixed_effect_inference import ( + InferenceParts, + contrast_rank, + contrast_tests, + f_test, + kenward_roger_parts, + satterthwaite_parts, +) from cytepool.stats._mixed_model_designs import ( MixedDesign, maximum_likelihood, mixed_design, ) -from cytepool.stats.distributions import t_sf -__all__ = ["FixedEffectTests", "MixedModelFit", "fit_mixed_model", "fixed_effect_tests"] +__all__ = [ + "FixedEffectFTest", + "FixedEffectTests", + "MixedModelFit", + "fit_mixed_model", + "fixed_effect_f_test", + "fixed_effect_tests", +] FloatArray = npt.NDArray[np.float64] IntArray = npt.NDArray[np.int64] @@ -268,6 +346,10 @@ class MixedModelFit: of an observation of the feature's average weight. mean_expression: ``(n_features,)`` unweighted mean of each feature's values (limma's ``Amean``, R's ``rowMeans`` to the last bit). + weights: The precision weights of the fit, as given: ``(n_samples,)`` sample + weights, ``(n_samples, n_features)`` observation weights, or None. + component_levels: ``(n_components, n_samples)`` the level of each sample in + the grouping factor of each variance component, numbered from 0. """ coefficients: FloatArray @@ -289,6 +371,8 @@ class MixedModelFit: residual_df: FloatArray mean_weights: FloatArray mean_expression: FloatArray + weights: FloatArray | None + component_levels: IntArray def __repr__(self) -> str: return ( @@ -311,28 +395,39 @@ class FixedEffectTests: contrasts: ``(n_coefficients, n_contrasts)`` weights of each contrast on the fit's coefficients. contrast_names: Names of the contrasts. - df_method: How the degrees of freedom were computed: ``"satterthwaite"``. + df_method: How the degrees of freedom were computed: ``"satterthwaite"`` or + ``"kenward_roger"``. + kenward_roger_weights: For Kenward and Roger's method, how its adjustment + took the weights: ``"average"`` or ``"precision"`` (see + :func:`fixed_effect_tests`); None for Satterthwaite's. estimate: ``(n_features, n_contrasts)`` estimates c' beta. standard_error: ``(n_features, n_contrasts)`` standard errors - (c' C c)^(1/2), with C the covariance matrix of the coefficients. + (c' C c)^(1/2), with C the covariance matrix of the coefficients: + the fit's (``cov_coefficients``) for Satterthwaite's method, Kenward + and Roger's adjusted one for theirs. t: ``(n_features, n_contrasts)`` t-statistics, the estimates over their standard errors. - df: ``(n_features, n_contrasts)`` Satterthwaite's degrees of freedom; NaN - where the asymptotic covariance of the variance parameters is not - defined. + df: ``(n_features, n_contrasts)`` degrees of freedom; NaN where they are + not defined, as where the asymptotic covariance of the variance + parameters is not. p_value: ``(n_features, n_contrasts)`` two-sided p-values of ``t`` on ``df`` degrees of freedom. + covariance: ``(n_features, n_contrasts, n_contrasts)`` the covariance + matrix of the estimates, L' C L for the contrasts L and the C of + ``standard_error``. """ fit: MixedModelFit contrasts: FloatArray contrast_names: tuple[str, ...] df_method: str + kenward_roger_weights: str | None estimate: FloatArray standard_error: FloatArray t: FloatArray df: FloatArray p_value: FloatArray + covariance: FloatArray def __repr__(self) -> str: return ( @@ -341,6 +436,53 @@ def __repr__(self) -> str: ) +@dataclass(frozen=True, eq=False, repr=False) +class FixedEffectFTest: + """F-tests that several contrasts of the fixed effects are all zero. + + One test per feature, of the hypothesis that every contrast is zero; the arrays + are read-only. + + Attributes: + fit: The mixed model fit. + contrasts: ``(n_coefficients, n_contrasts)`` weights of each contrast on the + fit's coefficients. + contrast_names: Names of the contrasts. + df_method: ``"satterthwaite"`` or ``"kenward_roger"``. + kenward_roger_weights: ``"average"`` or ``"precision"`` for Kenward and + Roger's method, as in :class:`FixedEffectTests`; None for + Satterthwaite's. + f_statistic: ``(n_features,)`` F-statistics; Kenward and Roger's are + scaled by ``scaling``. + df1: The numerator degrees of freedom, the number of linearly independent + contrasts. + df2: ``(n_features,)`` denominator degrees of freedom. + p_value: ``(n_features,)`` p-values of ``f_statistic`` on ``df1`` and + ``df2`` degrees of freedom. + scaling: ``(n_features,)`` the factor by which Kenward and Roger scale the + Wald statistic (pbkrtest's ``F.scaling``); 1 for Satterthwaite's + method. + """ + + fit: MixedModelFit + contrasts: FloatArray + contrast_names: tuple[str, ...] + df_method: str + kenward_roger_weights: str | None + f_statistic: FloatArray + df1: int + df2: FloatArray + p_value: FloatArray + scaling: FloatArray + + def __repr__(self) -> str: + return ( + f"FixedEffectFTest(n_features={self.f_statistic.shape[0]}, " + f"contrast_names={self.contrast_names!r}, df1={self.df1}, " + f"df_method={self.df_method!r})" + ) + + def fit_mixed_model( y: npt.ArrayLike, design: Design, @@ -425,7 +567,7 @@ def _fit( x, names, codes, n_levels, components = parts n_samples, n_coefficients = x.shape values = _response_matrix(y, n_samples) - weight_rows = _positive_weights(weights, values.shape) + weight_rows, given_weights = _positive_weights(weights, values.shape) n_features = values.shape[1] # The fixed effects in an orthonormal basis of the columns of x, x = q r: the @@ -533,6 +675,8 @@ def _fit( residual_df=_frozen(residual_df), mean_weights=_frozen(mean_weights), mean_expression=_frozen(observed_means(values)), + weights=None if given_weights is None else _frozen(given_weights), + component_levels=_frozen(codes, np.int64), ) @@ -541,15 +685,21 @@ def fixed_effect_tests( contrasts: Any = None, *, names: Sequence[str] | None = None, - df_method: Literal["satterthwaite"] = "satterthwaite", + df_method: Literal["satterthwaite", "kenward_roger"] = "satterthwaite", + kenward_roger_weights: Literal["average", "precision"] = "average", ) -> FixedEffectTests: - """t-tests of contrasts of the fixed effects with Satterthwaite's degrees of freedom. - - For a contrast c, the t-statistic is c' beta / (c' C c)^(1/2) and its degrees of - freedom are Satterthwaite's, as lmerTest's ``contest(fit, L, joint = FALSE)`` - and ``summary`` compute them; see the module docstring. Without contrasts, each - coefficient is tested, as in lmerTest's coefficient table. Kenward-Roger's - degrees of freedom, and F-tests of several contrasts at once, are planned. + """t-tests of contrasts of the fixed effects, by Satterthwaite or Kenward-Roger. + + For a contrast c, the t-statistic is c' beta / (c' C c)^(1/2). With + ``df_method="satterthwaite"``, C is the fit's covariance matrix and the degrees + of freedom are Satterthwaite's, as lmerTest's ``contest(fit, L, joint = FALSE)`` + and ``summary`` compute them. With ``"kenward_roger"``, C is Kenward and Roger's + adjusted covariance matrix, which allows for the estimation of the variance + parameters, and the degrees of freedom are theirs, as lmerTest's ``contest1D(fit, + L, ddf = "Kenward-Roger")`` computes them with pbkrtest's ``vcovAdj`` and + ``get_Lb_ddf``. The module docstring gives both methods. Without contrasts, each + coefficient is tested, as in lmerTest's coefficient table. + :func:`fixed_effect_f_test` tests several contrasts at once. Contrasts that name the coefficients they weigh are matched to ``fit.coefficient_names`` by name, as by @@ -559,14 +709,24 @@ def fixed_effect_tests( its index. Weights without names are taken in the order of the coefficients. Args: - fit: A fit from :func:`fit_mixed_model`. + fit: A fit from :func:`fit_mixed_model`; by REML for Kenward and Roger's + method. contrasts: The contrasts, one per column of a ``(n_coefficients, n_contrasts)`` matrix: a ``Contrasts``, a pandas DataFrame, a mapping of contrast name to weights, a matrix or a single vector. None tests every coefficient. names: Contrast names, when ``contrasts`` is not a mapping or a ``Contrasts``; see ``contrasts_fit``. - df_method: ``"satterthwaite"``, the only method so far. + df_method: ``"satterthwaite"`` (lmerTest's default) or + ``"kenward_roger"``. + kenward_roger_weights: How Kenward and Roger's adjustment and degrees of + freedom take a fit's weights: ``"average"`` computes them with every + observation at its feature's average weight, as variancePartition's + ``dream`` does with pbkrtest, which leaves out the weights of an lme4 fit + (lmerTest's ``ddf = "Kenward-Roger"`` gives the same for weights of mean + one); ``"precision"`` with each observation's own weight, the covariance + of the model fitted. The two agree without weights or with equal + weights; see the module docstring. Ignored for Satterthwaite's method. Returns: The estimates, standard errors, t-statistics, degrees of freedom and @@ -575,52 +735,139 @@ def fixed_effect_tests( Raises: TypeError: If ``fit`` is not a ``MixedModelFit``, or ``contrasts`` or ``names`` is not valid or is or holds a numpy masked array. - ValueError: If the contrasts do not match the coefficients, are not finite, - or ``df_method`` is unknown. + ValueError: If the contrasts do not match the coefficients or are not + finite, ``df_method`` or ``kenward_roger_weights`` is unknown, or + Kenward and Roger's method is asked of a maximum likelihood fit. """ - if not isinstance(fit, MixedModelFit): - raise TypeError(f"fit must be a MixedModelFit, got {type(fit).__name__}") - if df_method != "satterthwaite": - raise ValueError( - f"df_method must be 'satterthwaite', got {df_method!r}; Kenward-Roger's " - "degrees of freedom are not available yet" - ) - if contrasts is None: - reject_masked(names, "names") - if names is not None: - raise ValueError("names are given only with contrasts") - matrix = np.eye(len(fit.coefficient_names)) - contrast_names = fit.coefficient_names - else: - matrix, contrast_names = contrast_matrix( - contrasts, names, fit.coefficient_names - ) - estimate = fit.coefficients @ matrix - variance = np.einsum("jc,gjk,kc->gc", matrix, fit.cov_coefficients, matrix) - gradient = np.einsum( - "jc,gajk,kc->gca", matrix, fit.cov_coefficients_derivatives, matrix - ) - spread = np.einsum( - "gca,gab,gcb->gc", gradient, fit.variance_parameter_covariance, gradient + fit = _checked_fit(fit) + matrix, contrast_names = _contrasts_of(fit, contrasts, names) + parts, variant = _inference_parts(fit, df_method, kenward_roger_weights) + estimate, standard_error, t, df, p_value = contrast_tests( + fit.coefficients, parts, matrix ) - with np.errstate(divide="ignore", invalid="ignore"): - standard_error = np.sqrt(variance) - t = estimate / standard_error - df = np.where(spread > 0.0, 2.0 * variance**2 / spread, np.nan) - p_value = 2.0 * np.asarray(t_sf(np.abs(t), df)) return FixedEffectTests( fit=fit, contrasts=_frozen(matrix), contrast_names=tuple(contrast_names), - df_method=df_method, + df_method=parts.method, + kenward_roger_weights=variant, estimate=_frozen(estimate), standard_error=_frozen(standard_error), t=_frozen(t), df=_frozen(df), p_value=_frozen(p_value), + covariance=_frozen( + symmetric(np.einsum("jc,gjk,kd->gcd", matrix, parts.covariance, matrix)) + ), ) +def fixed_effect_f_test( + fit: MixedModelFit, + contrasts: Any, + *, + names: Sequence[str] | None = None, + df_method: Literal["satterthwaite", "kenward_roger"] = "satterthwaite", + kenward_roger_weights: Literal["average", "precision"] = "average", +) -> FixedEffectFTest: + """F-tests that several contrasts of the fixed effects are all zero. + + One test per feature of the hypothesis L beta = 0, for the contrasts in the + columns of L, such as every difference among three time points or the levels of + an interaction. With ``df_method="satterthwaite"`` it is lmerTest's + ``contestMD(fit, L)``: the Wald statistic over the number of independent + contrasts, whose denominator degrees of freedom combine those of independent + contrasts that span the hypothesis (Fai and Cornelius 1996; Kuznetsova et al. + 2017). With ``"kenward_roger"`` it is pbkrtest's ``KRmodcomp(fit, L)`` (and + lmerTest's ``contestMD(fit, L, ddf = "Kenward-Roger")``): the Wald statistic with + the adjusted covariance matrix, scaled, on Kenward and Roger's degrees of freedom. + The module docstring gives both. Contrasts that are linear combinations of the + others add nothing to the hypothesis: the numerator degrees of freedom are the + number of independent contrasts, by the rank rule of R's ``qr()``. + + Contrasts are matched to the coefficients as by :func:`fixed_effect_tests`. + + Args: + fit: A fit from :func:`fit_mixed_model`; by REML for Kenward and Roger's + method. + contrasts: The contrasts, one per column of a ``(n_coefficients, + n_contrasts)`` matrix, in any form :func:`fixed_effect_tests` takes. + names: Contrast names, when ``contrasts`` is not a mapping or a + ``Contrasts``. + df_method: ``"satterthwaite"`` or ``"kenward_roger"``. + kenward_roger_weights: ``"average"`` or ``"precision"``, as for + :func:`fixed_effect_tests`. + + Returns: + The F-statistic, its degrees of freedom and its p-value for every feature. + + Raises: + TypeError: As :func:`fixed_effect_tests`. + ValueError: As :func:`fixed_effect_tests`, or if every contrast is zero. + """ + fit = _checked_fit(fit) + matrix, contrast_names = _contrasts_of(fit, contrasts, names) + rank = contrast_rank(matrix) + if rank == 0: + raise ValueError( + "the contrasts are all zero, so they state no hypothesis to test; give " + "at least one contrast with a non-zero weight" + ) + parts, variant = _inference_parts(fit, df_method, kenward_roger_weights) + tests = f_test(fit.coefficients, parts, matrix, rank) + return FixedEffectFTest( + fit=fit, + contrasts=_frozen(matrix), + contrast_names=tuple(contrast_names), + df_method=parts.method, + kenward_roger_weights=variant, + f_statistic=_frozen(tests.statistic), + df1=tests.df1, + df2=_frozen(tests.df2), + p_value=_frozen(tests.p_value), + scaling=_frozen(tests.scaling), + ) + + +def _checked_fit(fit: Any) -> MixedModelFit: + if not isinstance(fit, MixedModelFit): + raise TypeError(f"fit must be a MixedModelFit, got {type(fit).__name__}") + return fit + + +def _contrasts_of( + fit: MixedModelFit, contrasts: Any, names: Sequence[str] | None +) -> tuple[FloatArray, tuple[str, ...]]: + """The contrast matrix on the fit's coefficients and the contrasts' names.""" + if contrasts is None: + reject_masked(names, "names") + if names is not None: + raise ValueError("names are given only with contrasts") + return np.eye(len(fit.coefficient_names)), fit.coefficient_names + matrix, contrast_names = contrast_matrix(contrasts, names, fit.coefficient_names) + return matrix, tuple(contrast_names) + + +def _inference_parts( + fit: MixedModelFit, df_method: Any, kenward_roger_weights: Any +) -> tuple[InferenceParts, str | None]: + """The parts of the tests by the method asked for, and how Kenward and Roger's + adjustment took the weights (None for Satterthwaite's method).""" + method = checked_choice(df_method, "df_method", ("satterthwaite", "kenward_roger")) + variant = checked_choice( + kenward_roger_weights, "kenward_roger_weights", ("average", "precision") + ) + if method == "satterthwaite": + return satterthwaite_parts(fit), None + if fit.method != "reml": + raise ValueError( + "Kenward and Roger's method is defined for REML fits, and this fit " + "maximised the likelihood; fit with method='reml', as lmerTest and " + "pbkrtest require (pbkrtest refits such a model by REML)" + ) + return kenward_roger_parts(fit, variant == "average"), variant + + # Arguments. @@ -651,10 +898,11 @@ def _response_matrix(y: npt.ArrayLike, n_samples: int) -> FloatArray: def _positive_weights( weights: npt.ArrayLike | None, shape: tuple[int, int] -) -> FloatArray: - """The weights as rows of samples, one per feature or one for all.""" +) -> tuple[FloatArray, FloatArray | None]: + """The weights as rows of samples, one per feature or one for all, and as given + (None without weights).""" if weights is None: - return np.ones((1, shape[0])) + return np.ones((1, shape[0])), None w = real_array(unmasked_array(weights, "weights"), "weights") if w.shape == (shape[0],): rows = w[np.newaxis, :] @@ -675,7 +923,7 @@ def _positive_weights( "not finite. A mixed model is fitted to every sample of a feature: lme4 " "gives an infinite REML criterion for a zero weight" ) - return _writable_c_array(rows, np.float64) + return _writable_c_array(rows, np.float64), np.asarray(w, dtype=np.float64) def _writable_c_array(values: npt.NDArray[Any], dtype: type) -> npt.NDArray[Any]: diff --git a/cytepool/variance_partitioning.py b/cytepool/variance_partitioning.py index 19b6fe3..9c3b4dc 100644 --- a/cytepool/variance_partitioning.py +++ b/cytepool/variance_partitioning.py @@ -121,8 +121,8 @@ from ._analysis_records import ( cube_record, - design_record, matched_samples, + mixed_design_record, software_versions, ) from ._columns import text_column @@ -912,7 +912,7 @@ def variance_partitioning( "versions": software_versions(), "cube": cube_record(cube), "layer": COUNTS_LAYER, - "design": _design_record(design, no_contrasts), + "design": mixed_design_record(design, no_contrasts), "settings": {**settings, **_FIXED_SETTINGS}, "groups": [_group_record(fit, sample_ids, settings) for fit in fits], }, @@ -1011,17 +1011,6 @@ def _summary_table(fits: list[_GroupFit]) -> Table: ) -def _design_record(design: Design, contrasts: Contrasts) -> dict[str, Any]: - """The design's record, with its random effects and the variables they group by.""" - record = design_record(design, contrasts) - record["random"] = [str(effect) for effect in design.random] - samples = design.samples - for effect in design.random: - for name in effect.grouping: - record["samples"].setdefault(name, samples[name].tolist()) - return record - - def _group_record( fit: _GroupFit, sample_ids: np.ndarray, settings: Mapping[str, Any] ) -> dict[str, Any]: diff --git a/docs/source/index.md b/docs/source/index.md index ffcba76..0eab43a 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -39,7 +39,7 @@ what CytePool offers for each today. | Does it replicate, and how should studies be combined? | Meta-analysis of every gene across studies, with the between-study variance, Q and I². Not yet: replication measures and shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md) | | Can patients be stratified from their profiles? | Not yet; cubes convert to AnnData and MuData for other tools | [Capability matrix](reference/capabilities.md) | | Which cell-cell interactions change? | Not yet | [Related software](reference/related_software.md) | -| Do effects depend on genotype, time, sex or context? | Interactions and F-tests in designs, likelihood ratio tests, random intercepts for repeated measures. Not yet: random slopes. QTL mapping is left to dedicated tools | [Designs, linear models and contrasts](statistics/linear_models.md), [Count models](statistics/count_models.md) | +| Do effects depend on genotype, time, sex or context? | Interactions and F-tests in designs, likelihood ratio tests, random intercepts for repeated measures, with F-tests and small-sample (Kenward-Roger) tests in mixed models. Not yet: random slopes. QTL mapping is left to dedicated tools | [Designs, linear models and contrasts](statistics/linear_models.md), [Count models](statistics/count_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | | How do cells respond to drugs and doses? | A dose as a factor or a covariate of a design. Not yet: dose-response fits, benchmark doses, synergy | [Capability matrix](reference/capabilities.md) | | Can public studies be pooled into one cube? | Not yet; cubes import the pseudobulk layouts of scanpy, decoupler and LIANA, and meta-analysis combines results across cubes | [Export a cube from another tool](guides/export_cubes.ipynb), [Meta-analysis](statistics/meta_analysis.md) | diff --git a/docs/source/reference/api/dream.md b/docs/source/reference/api/dream.md index 8a8f461..43c03f7 100644 --- a/docs/source/reference/api/dream.md +++ b/docs/source/reference/api/dream.md @@ -1,7 +1,7 @@ --- myst: html_meta: - description: Voom precision weights from the mean-variance trend of mixed model fits, and empirical Bayes moderated Satterthwaite t-tests of their fixed effects with topTable's z-statistic, as variancePartition's dream computes them, in CytePool. + description: Voom precision weights from the mean-variance trend of mixed model fits, and empirical Bayes moderated Satterthwaite or Kenward-Roger t-tests and F-tests of their fixed effects with topTable's z-statistic, as variancePartition's dream computes them, in CytePool. --- # dream @@ -10,10 +10,12 @@ myst: from mixed models of every gene, fitted by maximum likelihood, as variancePartition's `voomWithDreamWeights` does. {func}`~cytepool.stats.dream.moderated_fixed_effect_tests` moderates the -Satterthwaite t-tests of {func}`~cytepool.stats.mixed_models.fixed_effect_tests` by -empirical Bayes, as variancePartition's `eBayes` does for `dream` fits, and adds -`topTable`'s z-statistic. {func}`~cytepool.differential_expression.dream_voom` runs -these steps on every cell type of a cube. +Satterthwaite or Kenward-Roger t-tests of +{func}`~cytepool.stats.mixed_models.fixed_effect_tests` by empirical Bayes, as +variancePartition's `eBayes` does for `dream` fits, adds `topTable`'s z-statistic, and +tests the contrasts at once by the moderated F-statistic. +{func}`~cytepool.differential_expression.dream_voom` runs these steps on every cell +type of a cube, with dream's choice of the degrees of freedom. ```{eval-rst} .. automodule:: cytepool.stats.dream diff --git a/docs/source/reference/api/mixed_models.md b/docs/source/reference/api/mixed_models.md index 1aec696..f2a093b 100644 --- a/docs/source/reference/api/mixed_models.md +++ b/docs/source/reference/api/mixed_models.md @@ -1,7 +1,7 @@ --- myst: html_meta: - description: Fit a linear mixed model with random intercepts to every gene by REML or maximum likelihood, as lme4's lmer does, with fitted values and approximate residual degrees of freedom, and test contrasts with Satterthwaite's degrees of freedom, as lmerTest does. + description: Fit a linear mixed model with random intercepts to every gene by REML or maximum likelihood, as lme4's lmer does, with fitted values and approximate residual degrees of freedom, and test contrasts one at a time or several at once with Satterthwaite's degrees of freedom, as lmerTest does, or Kenward-Roger's, as pbkrtest does. --- # Mixed models diff --git a/docs/source/reference/bibliography.md b/docs/source/reference/bibliography.md index 4d46003..996d291 100644 --- a/docs/source/reference/bibliography.md +++ b/docs/source/reference/bibliography.md @@ -33,7 +33,8 @@ limma and its methods [@smyth2004; @smyth2005; @ritchie2015; @phipson2016; [@bates2015; @patterson1971; @harville1977; @henderson1975; @satterthwaite1946; @giesbrecht1985; @kuznetsova2017; @bertsekas1982; @hoffman2016; @hoffman2021; -@hoffman2026a] +@hoffman2026a]; Kenward and Roger's tests and F-tests [@kenward1997; @halekoh2014; +@fai1996] ## Gene sets, pathways and interactions diff --git a/docs/source/reference/capabilities.md b/docs/source/reference/capabilities.md index 98cfbf5..0835634 100644 --- a/docs/source/reference/capabilities.md +++ b/docs/source/reference/capabilities.md @@ -28,6 +28,7 @@ size; its sections give every check. | Model matrices from R formulas, named contrasts, design checks | {api}`Design`, {api}`model_matrix`, {api}`make_contrasts`, {api}`Design.check` | R 4.4.3, limma 3.62.1, lme4 2.0.6 | exact | [`design.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/design.py) | | Linear models of every gene, contrasts (limma's `lmFit`, `contrasts.fit`) | {api}`fit_linear_model`, {api}`contrasts_fit` | limma 3.62.1 | exact | [`linear_models.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/linear_models.py) | | Mixed models with random intercepts by REML or ML, Satterthwaite's t-tests | {api}`fit_mixed_model`, {api}`fixed_effect_tests` | lme4 2.0.6, lmerTest 3.2.1 | close | [`mixed_models.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/mixed_models.py) | +| Kenward-Roger's t-tests and F-tests of mixed models, Satterthwaite's F-tests (pbkrtest's `vcovAdj` and `KRmodcomp`, lmerTest's `contestMD`); dream's choice of degrees of freedom and its F-tests | {api}`fixed_effect_tests`, {api}`fixed_effect_f_test`, {api}`dream_voom`, {api}`moderated_fixed_effect_tests` | pbkrtest 0.5.3, lmerTest 3.1-3, lme4 1.1-38, variancePartition 1.36.2 | exact at given estimates; close at the optimum and for dream | [`kenward_roger_and_f_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/kenward_roger_and_f_tests.py) | | voom precision weights, limma's span rule | {api}`voom`, {api}`choose_lowess_span` | limma 3.62.1, edgeR 4.4.0 | exact; close on rarefied counts | [`voom.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/voom.py) | | Sample quality weights (`arrayWeights`, `voomWithQualityWeights`) | {api}`sample_quality_weights`, {api}`voom_with_quality_weights` | limma 3.62.1, edgeR 4.4.0 | close | [`sample_quality_weights.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/sample_quality_weights.py) | | Negative binomial models of every gene (edgeR's `glmFit`), deviances | {api}`fit_negative_binomial_model`, {api}`negative_binomial_deviance`, {api}`negative_binomial_unit_deviance` | edgeR 4.4.0, mpmath 1.3.0 | close | [`negative_binomial_fits.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/negative_binomial_fits.py) | @@ -55,8 +56,7 @@ size; its sections give every check. These methods come up in the questions that studies ask ([Overview](../index.md)) and have no CytePool function yet: -- random slopes, Kenward-Roger's degrees of freedom and F-tests of several contrasts in - mixed models; +- random slopes in mixed models; - composition tests with precision weights or random effects (crumblr), and Bayesian count models of composition; - pathway and transcription-factor activity scores; diff --git a/docs/source/reference/citing.md b/docs/source/reference/citing.md index 5f16deb..01d7d5d 100644 --- a/docs/source/reference/citing.md +++ b/docs/source/reference/citing.md @@ -20,7 +20,7 @@ the credit for the methods belongs to their authors. | `quality_weights=True`, {api}`sample_quality_weights`, {api}`voom_with_quality_weights` | [@ritchie2006; @liu2015] | | `prior_weights="n_cells"` | dreamlet [@hoffman2026a] | | {api}`dream_voom`, {api}`voom_with_dream_weights`, {api}`moderated_fixed_effect_tests` | dream [@hoffman2021] | -| {api}`fit_mixed_model`, {api}`fixed_effect_tests` | lme4 [@bates2015]; Satterthwaite's tests as lmerTest computes them [@kuznetsova2017] | +| {api}`fit_mixed_model`, {api}`fixed_effect_tests`, {api}`fixed_effect_f_test` | lme4 [@bates2015]; Satterthwaite's tests as lmerTest computes them [@kuznetsova2017]; Kenward and Roger's [@kenward1997] as pbkrtest computes them [@halekoh2014] | | {api}`treat` | [@mccarthy2009] | | {api}`fit_quasi_likelihood`, {api}`quasi_likelihood_f_test` | edgeR 4 [@chen2025]; quasi-likelihood [@lund2012] | | {api}`fit_negative_binomial_model`, {api}`common_dispersion` | [@mccarthy2012] | diff --git a/docs/source/reference/glossary.md b/docs/source/reference/glossary.md index 6a0c260..ccea6d6 100644 --- a/docs/source/reference/glossary.md +++ b/docs/source/reference/glossary.md @@ -65,6 +65,10 @@ pooled per group. A bulk cube has one group. sampling error, measured by the between-study variance tau², Cochran's Q and I². [Meta-analysis](../statistics/meta_analysis.md) +**Kenward-Roger's method.** An adjustment of the standard errors of a mixed model's +estimates, and of the degrees of freedom of their tests, for the uncertainty of its +variance components in small samples. [Mixed models, REML and dream](../statistics/mixed_models.md) + **Layer.** One statistic of a cube, summed over cells: counts, detected cells, or sums of normalized and log-normalized values. [The cube and its additive layers](../statistics/cubes.md) diff --git a/docs/source/references.bib b/docs/source/references.bib index 423bdfa..238cc71 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -436,6 +436,17 @@ @book{dunn2018 doi = {10.1007/978-1-4419-0118-7}, } +@article{fai1996, + author = {Fai, Alex Hrong-Tai and Cornelius, Paul L.}, + title = {Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments}, + journal = {Journal of Statistical Computation and Simulation}, + year = {1996}, + volume = {54}, + number = {4}, + pages = {363--378}, + doi = {10.1080/00949659608811740}, +} + @article{fisher1922, author = {Fisher, R. A.}, title = {On the interpretation of χ² from contingency tables, and the calculation of P}, @@ -571,6 +582,17 @@ @article{hafner2025 doi = {10.1093/bib/bbaf397}, } +@article{halekoh2014, + author = {Halekoh, Ulrich and Højsgaard, Søren}, + title = {A Kenward-Roger approximation and parametric bootstrap methods for tests in linear mixed models: the R package pbkrtest}, + journal = {Journal of Statistical Software}, + year = {2014}, + volume = {59}, + number = {9}, + pages = {1--32}, + doi = {10.18637/jss.v059.i09}, +} + @article{hartung1999, author = {Hartung, Joachim}, title = {An alternative method for meta-analysis}, @@ -828,6 +850,17 @@ @article{kang2018 doi = {10.1038/nbt.4042}, } +@article{kenward1997, + author = {Kenward, Michael G. and Roger, James H.}, + title = {Small sample inference for fixed effects from restricted maximum likelihood}, + journal = {Biometrics}, + year = {1997}, + volume = {53}, + number = {3}, + pages = {983--997}, + doi = {10.2307/2533558}, +} + @article{kendall1945, author = {Kendall, M. G.}, title = {The treatment of ties in ranking problems}, diff --git a/docs/source/statistics/choosing_a_method.md b/docs/source/statistics/choosing_a_method.md index d9d5326..d1f42d5 100644 --- a/docs/source/statistics/choosing_a_method.md +++ b/docs/source/statistics/choosing_a_method.md @@ -17,7 +17,7 @@ finds differences that do not hold across donors [@squair2021; @zimmerman2021; | --- | --- | --- | | Conditions, with covariates or a pairing as fixed effects, such as `~ donor + condition` | {api}`limma_voom` | Log-CPM values with voom's precision weights, a weighted linear model, empirical Bayes moderation [@law2014; @ritchie2015] | | The same, with library sizes within about threefold of each other | {api}`limma_trend` | Log-CPM values, a linear model, a prior variance that follows the expression level | -| Several samples of each donor, and a variable that differs between donors, such as disease or sex | {api}`dream_voom` | voom's weights from mixed models, a random intercept per donor, Satterthwaite's t-tests, empirical Bayes [@hoffman2021] | +| Several samples of each donor, and a variable that differs between donors, such as disease or sex | {api}`dream_voom` | voom's weights from mixed models, a random intercept per donor, Kenward-Roger's or Satterthwaite's t-tests as dream chooses them, empirical Bayes [@hoffman2021] | Each runs on every group (cell type) of a cube, one at a time, and returns a {api}`DifferentialExpressionResult` with a table for every cell type and contrast diff --git a/docs/source/statistics/mixed_models.md b/docs/source/statistics/mixed_models.md index 77ca423..ccf149c 100644 --- a/docs/source/statistics/mixed_models.md +++ b/docs/source/statistics/mixed_models.md @@ -1,5 +1,5 @@ --- -description: Linear mixed models with random intercepts fitted to every gene by REML, Satterthwaite's degrees of freedom, dream's precision weights and moderated tests for repeated measures designs, and variance partitioning. +description: Linear mixed models with random intercepts fitted to every gene by REML, Satterthwaite's and Kenward-Roger's degrees of freedom, F-tests of several contrasts, dream's precision weights and moderated tests for repeated measures designs, and variance partitioning. --- # Mixed models, REML and dream @@ -61,6 +61,49 @@ is their asymptotic covariance matrix. {api}`fixed_effect_tests` gives them in derivatives where lmerTest takes numerical ones. With every component at zero, $\nu = n - p$, the degrees of freedom of the fixed-effects model. +## Kenward and Roger's method + +With few samples the variance components are estimated imprecisely, and $C$, computed +as if they were known, understates the variance of the estimates. Kenward and Roger +[@kenward1997] adjust it: with $\Sigma = \sum_i \phi_i G_i$ the covariance of the +observations, linear in the variance parameters $\phi$, +$P_i = -X'\Sigma^{-1} G_i \Sigma^{-1} X$, +$Q_{ij} = X'\Sigma^{-1} G_i \Sigma^{-1} G_j \Sigma^{-1} X$ and $W$ the inverse of +the expected information of $\phi$ under REML, + +$$ +C_A = C + 2\,C \Big(\sum_{ij} W_{ij} \big(Q_{ij} - P_i C P_j\big)\Big) C . +$$ + +`fixed_effect_tests(fit, df_method="kenward_roger")` divides each estimate by its +standard error from $C_A$, with Satterthwaite's formula for the degrees of freedom on +the expected information, as pbkrtest computes them [@halekoh2014] and lmerTest's +`ddf = "Kenward-Roger"` reports them. In balanced designs the adjustment vanishes and +the tests are those of the analysis of variance. A variance component estimated at +zero takes part, which lowers the degrees of freedom below Satterthwaite's. + +With precision weights, pbkrtest builds $\Sigma$ without lme4's prior weights while +it takes $C$ from the weighted fit. dream divides each gene's weights by their mean, +so its tests treat every observation as having its gene's average weight. +`kenward_roger_weights="average"`, the default, computes these tests; `"precision"` +uses each observation's own weight, which is Kenward and Roger's method for the +weighted model. With voom's weights the first can make standard errors smaller than +the unadjusted ones, and the second never does; the [Validation](../validation.md) +page compares them. + +## Several contrasts at once + +{api}`fixed_effect_f_test` tests that several contrasts are all zero, such as every +difference among three time points, in a {api}`FixedEffectFTest`. With Satterthwaite's +method, as lmerTest's `contestMD`, the contrasts are replaced by the eigenvectors of +their covariance matrix, whose estimates are uncorrelated; $F$ is the mean of their +squared t-statistics, and its denominator degrees of freedom $2E/(E-r)$, with +$E = \sum_m \nu_m / (\nu_m - 2)$ over their degrees of freedom, match the mean of +$F$ [@fai1996; @kuznetsova2017]. With Kenward and Roger's method, as pbkrtest's +`KRmodcomp`, the Wald statistic with $C_A$ is scaled and referred to an $F$ +distribution whose first two moments it matches [@kenward1997]. The numerator degrees +of freedom $r$ are the number of linearly independent contrasts. + ## dream dream is limma-voom with the mixed model in place of the linear model: @@ -68,18 +111,34 @@ dream is limma-voom with the mixed model in place of the linear model: 1. {api}`voom_with_dream_weights`: voom's precision weights, from the mean-variance trend of the mixed model of every gene, fitted by maximum likelihood, read at the fitted values with the predicted random effects. -2. {api}`fit_mixed_model` by REML with those weights, and {api}`fixed_effect_tests`. +2. {api}`fit_mixed_model` by REML with those weights, and {api}`fixed_effect_tests` + with dream's degrees of freedom: by default Kenward and Roger's for fewer than 20 + samples, as pbkrtest computes them with every observation at the average weight, + and Satterthwaite's otherwise. 3. {api}`moderated_fixed_effect_tests`: empirical Bayes moderation of each gene's residual variance on the approximate residual degrees of freedom of the mixed model, $\operatorname{tr}((I - H)'(I - H))$ for its hat matrix $H$; the moderated - t-statistic has Satterthwaite's degrees of freedom plus the prior's, and the result + t-statistic has the degrees of freedom of its test plus the prior's, and the result ({api}`ModeratedFixedEffectTests`) adds a z-statistic with the same p-value, so that - tests on different degrees of freedom compare. - -{api}`dream_voom` runs these steps on every cell type of a cube. As variancePartition -does, it moderates the residual variance of an observation of each gene's average -weight, which keeps the mean-variance trend that voom's weights remove; `trend=True` -lets the prior follow it. + tests on different degrees of freedom compare. Its moderated F-statistic, + $t' R^{-1} t / r$ for the moderated t-statistics and the correlation $R$ of their + estimates, tests several contrasts at once, as `topTable` does. + +{api}`dream_voom` runs these steps on every cell type of a cube, with +`df_method="adaptive"` choosing the degrees of freedom in each cell type from the +samples left in it, and `f_tests` naming sets of contrasts to test at once. As +variancePartition does, it moderates the residual variance of an observation of each +gene's average weight, which keeps the mean-variance trend that voom's weights remove; +`trend=True` lets the prior follow it. + +dream switches from Kenward and Roger's method to Satterthwaite's at 20 samples +(its help page says "<= 20"; with 20 it uses Satterthwaite's), and for a gene whose +adjusted covariance matrix is ill conditioned, which depends on the units of +continuous covariates. Its `topTable` computes the p-values of an F-test of every gene +on the degrees of freedom of the first gene; here each gene has its own, those of $r$ +uncorrelated moderated t-statistics combined as lmerTest combines those of its +F-tests, so that an F-test of one contrast is its t-test. The +[Validation](../validation.md) page reports both differences and their sizes. ## Variance partitioning @@ -131,8 +190,7 @@ whether variables are missing from the design. ## Not yet implemented -Kenward-Roger's degrees of freedom (dream's default for 20 samples or fewer), F-tests of -several contrasts at once and random slopes are planned; random slopes would also -partition the variance within subsets of the samples, such as the variance between -donors within each tissue. A gene that the fixed and random effects reproduce exactly -has no REML estimate; it is left unfitted rather than given a variance of zero. +Random slopes are planned; they would also partition the variance within subsets of +the samples, such as the variance between donors within each tissue. A gene that the +fixed and random effects reproduce exactly has no REML estimate; it is left unfitted +rather than given a variance of zero. diff --git a/docs/source/validation.md b/docs/source/validation.md index f72ec21..452ecfd 100644 --- a/docs/source/validation.md +++ b/docs/source/validation.md @@ -16,13 +16,15 @@ exits with an error when any check falls outside its tolerance. This page summarises the tables the scripts printed when they were run against the package in this revision of the repository: on **2026-10-26** for `plotting.py` and -`variance_partitioning.py` and `gene_set_files_and_ora.py`; on +`variance_partitioning.py` and `gene_set_files_and_ora.py` and +`kenward_roger_and_f_tests.py`, and for `dream.py`, which printed the same table again +after changes to the code it checks; on **2026-10-25** for the six scripts on cubes and their files (`cube.py`, `pool_cells.py`, `cube_files.py`, `cube_reads.py`, `cube_fingerprints.py` and `anndata_conversion.py`), for `linear_models.py`, `mixed_models.py`, `voom.py`, `empirical_bayes.py`, `sample_quality_weights.py`, `negative_binomial_fits.py`, `quasi_likelihood.py`, `deseq2_tests.py`, -`differential_expression.py`, `differential_expression_options.py`, `dream.py`, +`differential_expression.py`, `differential_expression_options.py`, `interferon_stimulated_pbmcs.py`, `gene_set_tests.py`, `rotation_gene_set_tests.py`, `composition_tests.py` and `meta_analysis.py`, and for `design.py`, which printed the same table again after changes to the code it checks; and on **2026-10-24** for the @@ -58,6 +60,7 @@ covers and the largest difference among them. Run a script to see every check. | [`differential_expression.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression.py) | Differential expression | limma 3.62.1, edgeR 4.4.0 | 32 | 28 pass, 4 deliberate differences | 1.8 s | 189 MiB | | [`differential_expression_options.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression_options.py) | Options for pseudobulk differential expression | limma 3.62.1, edgeR 4.4.0 | 66 | 66 pass | 5.0 s | 182 MiB | | [`dream.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/dream.py) | dream | variancePartition 1.36.2, limma 3.62.1, lme4 1.1-38, lmerTest 3.1-3 | 63 | 53 pass, 5 information, 5 deliberate differences | 30.5 s | 390 MiB | +| [`kenward_roger_and_f_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/kenward_roger_and_f_tests.py) | Kenward-Roger's tests and F-tests | R 4.4.3, variancePartition 1.36.2, limma 3.62.1, lme4 1.1-38, lmerTest 3.1-3, pbkrtest 0.5.3 | 143 | 115 pass, 16 information, 12 deliberate differences | 3 min 59 s | 390 MiB | | [`variance_partitioning.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/variance_partitioning.py) | Variance partitioning | variancePartition 1.36.2, lme4 1.1-38, ggplot2 4.0.3 | 31 | 22 pass, 6 information, 3 deliberate differences | 63.1 s | 394 MiB | | [`interferon_stimulated_pbmcs.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/interferon_stimulated_pbmcs.py) | Differential expression on published data | limma 3.62.1, edgeR 4.4.0 | 41 | 17 pass, 24 information | 15.2 s | 617 MiB | | [`composition_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/composition_tests.py) | Cell-type composition | speckle 1.6.0, limma 3.62.1 | 32 | 25 pass, 7 deliberate differences | 13.1 s | 860 MiB | @@ -65,7 +68,7 @@ covers and the largest difference among them. Run a script to see every check. | [`plotting.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/plotting.py) | Figures | limma 3.62.1, edgeR 4.4.0, DESeq2 1.46.0, metafor 5.0.1, speckle 1.6.0 | 26 | 23 pass, 2 information, 1 deliberate difference | 25.4 s | 810 MiB | | [`distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) | Distributions and special functions | mpmath 1.3.0 | 46 | 46 pass | 6 min 39 s | 221 MiB | | [`smoothing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/smoothing.py) | LOWESS and spline bases | R 4.4.3 | 90 | 87 pass, 1 information, 2 deliberate differences | 1.5 s | 175 MiB | -| **All** | | | **2,007** | **1,736 pass, none failed; 156 information, 114 deliberate differences, 1 control** | | | +| **All** | | | **2,150** | **1,851 pass, none failed; 172 information, 126 deliberate differences, 1 control** | | | Run time is the wall time of the whole script and peak memory the largest resident set of any of its processes, reference programs included, both measured with GNU @@ -1314,7 +1317,9 @@ CytePool refuses the two together as well. compares `voom_with_dream_weights`, `moderated_fixed_effect_tests` and `dream_voom`, and the maximum likelihood fits, fitted values and approximate residual degrees of freedom of `fit_mixed_model` behind them, with variancePartition's -`voomWithDreamWeights`, `dream(ddf = "Satterthwaite")`, `eBayes` and `topTable`, on +`voomWithDreamWeights`, `dream(ddf = "Satterthwaite")`, `eBayes` and `topTable` +(`dream_voom` with `df_method="satterthwaite"`; [the next section](#validation-kenward-roger) +checks its default, the other choices of degrees of freedom and the F-tests), on negative binomial counts of 60 genes from fixed seeds (20261025 to 20261029): a random donor intercept with a condition, an age and the contrast `conditionstim - age` (10 donors x 3 samples, with two hypervariable genes), moderated also with `robust = TRUE` @@ -1398,6 +1403,102 @@ Deliberate differences: model failed"); `dream_voom` reports the status `random_effects_not_estimable` and goes on with the other cell types. +(validation-kenward-roger)= +## Kenward-Roger's tests and F-tests + +[`validation/kenward_roger_and_f_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/kenward_roger_and_f_tests.py) +compares `fixed_effect_tests(df_method="kenward_roger")` with lmerTest's `summary` and +`contest1D` (`ddf = "Kenward-Roger"`) and pbkrtest's `vcovAdj`, and +`fixed_effect_f_test` with pbkrtest's `KRmodcomp` and lmerTest's `contestMD` (by +Kenward-Roger and by Satterthwaite), gene by gene, in 9 cases of 24 genes from fixed +seeds: a time course of three time points within 6 donors with an age per donor, with +the F-tests of two contrasts and of three that span two dimensions; unbalanced donors +with 1 to 5 samples; crossed donor and batch intercepts; nested intercepts +(`(1 | donor/sample)`); a design simulated without donor variance, so that 15 of the 24 +fits lie on the boundary; and observation weights like voom's, with one and with two +random intercepts, each taken both ways: every observation at the average weight, +against lme4 fits whose weights are divided by their mean, as dream fits them, and the +precision weights, against pbkrtest on the same model without weights, its values, +fixed effects and random effects multiplied by the square roots of the weights, which is +the weighted model. It then compares dream with variancePartition's `dream` with +`ddf = "Kenward-Roger"`, `"adaptive"` and `"Satterthwaite"`, `eBayes` and `topTable`, +on negative binomial counts of 40 genes in studies of 12, 18, 19, 20 and 24 samples (a +time course of three time points and an age), including `topTable`'s F-test of time, +and `dream_voom` with its defaults on two cubes: cell types left with 12 and with 24 +samples, and one of 15 samples whose age, 32 years on average and about 1 apart, makes +the adjusted covariance matrices ill conditioned. 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 and minqa 1.2.8. + +| Check family | Reference | Largest difference | Tolerance | Outcome | +| --- | --- | --- | --- | --- | +| Adjusted covariance matrix of the coefficients, relative to each gene's largest entry | pbkrtest `vcovAdj` | 1.4e-13 | 1e-10 | 9 pass | +| Kenward-Roger t-tests of every coefficient and a contrast: estimates (relative to the larger of their absolute value and standard error), standard errors and t-statistics (relative to max(\|t\|, 1)) | lmerTest `summary`, `contest1D` | 1.7e-12 | 1e-10 | 9 pass | +| Their df and p-values | lmerTest `summary`, `contest1D` | 1.8e-12 | 1e-10 | 9 pass | +| Kenward-Roger F-tests: F, numerator and denominator df, scaling and p-value, and lmerTest's F and df | pbkrtest `KRmodcomp`, lmerTest `contestMD` | 8.2e-13 | 1e-10 | 13 pass | +| Satterthwaite F-statistics | lmerTest `contestMD` | 8.2e-13 | 1e-10 | 13 pass | +| Satterthwaite F-tests: denominator df and p-values | lmerTest `contestMD` | 1.9e-7 | 1e-5 | 13 pass | +| Kenward-Roger standard errors and df at lme4's own optimum | lmerTest `summary` | 1.2e-7 | 1e-6 | 9 pass | +| dream: precision weights | `voomWithDreamWeights` | 2.1e-8 | 1e-6 | 5 pass | +| dream with each `ddf`: standard errors and df of the tests | `dream` | 8.0e-8 | 1e-6 | 9 pass | +| Moderated t-statistics, B-statistics and z-statistics (relative to max(\|x\|, 1)), total df and p-values | `eBayes`, `topTable` | 5.0e-8 | 1e-6 | 9 pass | +| Moderated F-statistics of time | `topTable` | 4.1e-8 | 1e-6 | 11 pass | +| `dream_voom`: samples used and features tested | the rules of `dream_voom` | 0 differences | 0 | 3 pass | +| `dream_voom` with its defaults: t, B and z-statistics, df and p-values of every contrast | `dream`, `eBayes`, `topTable` | 3.3e-7 | 1e-6 | 3 pass | +| Fits with a variance at zero; the two treatments of the weights; the default against Satterthwaite's method | CytePool | below | n/a | 16 information | +| Deliberate differences, listed below | variancePartition | 12 inputs | n/a | 12 deliberate differences | + +At CytePool's estimates, which the script passes to lme4 as starting values with one +evaluation allowed, as for [mixed models](#validation-mixed-models), pbkrtest computes +Kenward and Roger's quantities from dense matrices without numerical derivatives, so +they are compared to the exact tier; lmerTest's Satterthwaite F-tests differentiate the +deviance numerically, compared to 1e-5. At lme4's own optimum, with Powell's BOBYQA and +its trust region shrunk to 1e-10 as in the section on dream, and in dream's fits, which +the script runs with that optimiser, the estimates agree to the close tier. P-values are +scaled by max(p, 1e-300) max(1, \|log p\|), as for voom. + +The checks confirm what probing pbkrtest, lmerTest and variancePartition as black boxes +found, and the documentation of `cytepool.stats.mixed_models` and +`cytepool.stats.dream` lists: pbkrtest takes the inverse of the expected information +of REML, includes a variance component estimated at zero, which 15 fits of the +boundary case have, and builds the covariance matrix of the observations without lme4's +prior weights while it takes the covariance matrix of the coefficients from the +weighted fit; lmerTest's `contest1D` and `summary` use the adjusted standard errors with +pbkrtest's degrees of freedom; `dream` uses Kenward and Roger's method with the adjusted +standard errors for 19 samples and Satterthwaite's for 20 under `ddf = "adaptive"`, and +`topTable`'s F-statistic is limma's, from the moderated t-statistics and each gene's +covariance matrix of the coefficients, the adjusted one under Kenward and Roger's +method. In the ill-conditioned cell type `dream` tested 24 genes by Satterthwaite's +method because R's estimate of the condition number of their adjusted covariance +exceeded 1e6, and `dream_voom` the same 24, from the exact condition number. + +The information rows give the size of the two treatments of the weights. Against the +precision weights, every observation at the average weight gave standard errors 0.976 +times as large at the median (0.51 to 1.05) and df 0.91 times (0.26 to 3.7) with one +random intercept, and 0.977 (0.81 to 1.02) and 0.89 (0.27 to 3.4) with two; its standard +errors were below the unadjusted ones for 61% and 33% of the genes, the precision +weights' for none. And they give how dream's default changes the results of small +studies against Satterthwaite's method, with which `dream_voom` tested them before: +with 12 samples, p-values 1.12 times as large at the median (0.97 to 1.66 for 90% of the +genes), with standard errors 1.06 times as large at the median and smaller for 20% of +the genes; with 18 samples, 1.00 (0.91 to 1.07), and standard errors 1.00, smaller for +48%; in the cube's cell type with 12 samples, 0.99 (0.73 to 1.12), with as many genes +below an adjusted p-value of 0.05 (12) as before. + +Deliberate differences: + +- The p-values of `topTable`'s F-tests: dream computes those of every gene on the total + degrees of freedom of the first gene for the first coefficient listed, which the + script reproduces to 3e-14 in the 11 F-tests; CytePool gives each gene its own (see + `cytepool.stats.dream`). Its p-values were 0.14 to 11.5 times dream's. +- Which genes fall back to Satterthwaite's method: dream decides from R's estimate of + the condition number (`kappa`), CytePool from the exact one; they can differ within + about 2% of 1e6. Here they decided alike for every gene tested, and the tests agree + to 3.3e-7. + +dream also computes `F.std` with the number of coefficients where they are linearly +dependent, and CytePool with their rank, the numerator degrees of freedom of the F-test; +the F-tests here have independent coefficients. + (validation-variance-partitioning)= ## Variance partitioning @@ -1989,6 +2090,7 @@ RSCRIPT=$REF/r/bin/Rscript METAFOR_RSCRIPT=$REF/r-metafor/bin/Rscript \ SPECKLE_RSCRIPT=$REF/r-speckle/bin/Rscript python validation/plotting.py # variancePartition, in an R environment of its own with lme4 1.1-38 (below). RSCRIPT=$REF/r-dream/bin/Rscript python validation/dream.py +RSCRIPT=$REF/r-dream/bin/Rscript python validation/kenward_roger_and_f_tests.py RSCRIPT=$REF/r-dream/bin/Rscript python validation/variance_partitioning.py python validation/distributions.py REFERENCE_PYTHON=$REF/python-scanpy/bin/python python validation/cube.py diff --git a/tests/data/kenward_roger_and_f_tests/boundary.json b/tests/data/kenward_roger_and_f_tests/boundary.json new file mode 100644 index 0000000..4b050ae --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/boundary.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D02", "D02", "D02", "D03", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [41.3, 41.3, 41.3, 54.3, 54.3, 54.3, 27.6, 27.6, 27.6, 61.3, 61.3, 61.3, 73.1, 73.1, 73.1, 61.2, 61.2, 61.2]}, "formula": "~ time + age", "random": "(1 | donor)", "variant": "none"} diff --git a/tests/data/kenward_roger_and_f_tests/boundary.npz b/tests/data/kenward_roger_and_f_tests/boundary.npz new file mode 100644 index 0000000..9175c4d Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/boundary.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/crossed.json b/tests/data/kenward_roger_and_f_tests/crossed.json new file mode 100644 index 0000000..965dc89 --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/crossed.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D01", "D02", "D02", "D02", "D02", "D03", "D03", "D03", "D03", "D04", "D04", "D04", "D04", "D05", "D05", "D05", "D05", "D06", "D06", "D06", "D06", "D07", "D07", "D07", "D07", "D08", "D08", "D08", "D08"], "batch": ["B03", "B02", "B04", "B03", "B01", "B02", "B03", "B01", "B01", "B04", "B01", "B03", "B02", "B04", "B03", "B01", "B03", "B01", "B02", "B04", "B01", "B02", "B02", "B01", "B04", "B01", "B01", "B02", "B04", "B03", "B01", "B01"], "condition": ["ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim"], "x": [0.435, -0.893, 0.596, -0.539, 0.051, 1.5, 0.751, 2.569, -0.907, 0.676, -0.687, -0.186, -0.914, 1.364, -0.189, -0.129, 0.925, 0.182, -1.378, 1.289, 0.869, -0.401, -0.842, -2.832, 0.406, 0.157, -1.432, -0.614, 1.902, 0.766, 0.043, -1.843]}, "formula": "~ condition + x", "random": "(1 | donor) + (1 | batch)", "variant": "none"} diff --git a/tests/data/kenward_roger_and_f_tests/crossed.npz b/tests/data/kenward_roger_and_f_tests/crossed.npz new file mode 100644 index 0000000..5f91f55 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/crossed.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/cube.json b/tests/data/kenward_roger_and_f_tests/cube.json new file mode 100644 index 0000000..37bebf7 --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/cube.json @@ -0,0 +1 @@ +{"donor": ["D01", "D01", "D01", "D02", "D02", "D02", "D03", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06", "D07", "D07", "D07", "D08", "D08", "D08"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [0.909, 0.437, -1.713, 0.814, 0.786, 0.493, -0.126, 0.21, -0.518, 0.181, 0.612, 2.118, 1.044, -0.181, -0.04, -1.024, -0.754, -0.501, -0.781, 0.762, -0.91, 1.22, 1.076, 0.099], "sample_id": ["S00", "S01", "S02", "S03", "S04", "S05", "S06", "S07", "S08", "S09", "S10", "S11", "S12", "S13", "S14", "S15", "S16", "S17", "S18", "S19", "S20", "S21", "S22", "S23"]} diff --git a/tests/data/kenward_roger_and_f_tests/cube.npz b/tests/data/kenward_roger_and_f_tests/cube.npz new file mode 100644 index 0000000..37d2f1a Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/cube.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/nested.json b/tests/data/kenward_roger_and_f_tests/nested.json new file mode 100644 index 0000000..cd0b54f --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/nested.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D01", "D01", "D01", "D02", "D02", "D02", "D02", "D02", "D02", "D03", "D03", "D03", "D03", "D03", "D03", "D04", "D04", "D04", "D04", "D04", "D04"], "sample": ["S01", "S01", "S02", "S02", "S03", "S03", "S01", "S01", "S02", "S02", "S03", "S03", "S01", "S01", "S02", "S02", "S03", "S03", "S01", "S01", "S02", "S02", "S03", "S03"], "time": ["t0", "t0", "t1", "t1", "t2", "t2", "t0", "t0", "t1", "t1", "t2", "t2", "t0", "t0", "t1", "t1", "t2", "t2", "t0", "t0", "t1", "t1", "t2", "t2"]}, "formula": "~ time", "random": "(1 | donor/sample)", "variant": "none"} diff --git a/tests/data/kenward_roger_and_f_tests/nested.npz b/tests/data/kenward_roger_and_f_tests/nested.npz new file mode 100644 index 0000000..daed079 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/nested.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/study_n12.json b/tests/data/kenward_roger_and_f_tests/study_n12.json new file mode 100644 index 0000000..6b597c4 --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/study_n12.json @@ -0,0 +1 @@ +{"donor": ["D01", "D01", "D02", "D02", "D03", "D03", "D04", "D04", "D05", "D05", "D06", "D06"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [-1.789, -1.835, 0.981, 0.669, 0.067, 0.012, 0.301, -0.182, -0.596, 0.908, 0.531, 0.102]} diff --git a/tests/data/kenward_roger_and_f_tests/study_n12.npz b/tests/data/kenward_roger_and_f_tests/study_n12.npz new file mode 100644 index 0000000..cb9b4d4 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/study_n12.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/study_n20.json b/tests/data/kenward_roger_and_f_tests/study_n20.json new file mode 100644 index 0000000..815703f --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/study_n20.json @@ -0,0 +1 @@ +{"donor": ["D01", "D01", "D01", "D01", "D02", "D02", "D02", "D02", "D03", "D03", "D03", "D03", "D04", "D04", "D04", "D04", "D05", "D05", "D05", "D05"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1"], "age": [0.841, 1.531, -0.245, -2.037, -1.377, -0.828, 0.685, -0.466, 1.128, 0.122, -0.81, 0.12, -1.801, 0.0, 0.085, -0.197, -1.501, 0.866, -0.311, 1.209]} diff --git a/tests/data/kenward_roger_and_f_tests/study_n20.npz b/tests/data/kenward_roger_and_f_tests/study_n20.npz new file mode 100644 index 0000000..e48cfc6 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/study_n20.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/time_course.json b/tests/data/kenward_roger_and_f_tests/time_course.json new file mode 100644 index 0000000..4b050ae --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/time_course.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D02", "D02", "D02", "D03", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [41.3, 41.3, 41.3, 54.3, 54.3, 54.3, 27.6, 27.6, 27.6, 61.3, 61.3, 61.3, 73.1, 73.1, 73.1, 61.2, 61.2, 61.2]}, "formula": "~ time + age", "random": "(1 | donor)", "variant": "none"} diff --git a/tests/data/kenward_roger_and_f_tests/time_course.npz b/tests/data/kenward_roger_and_f_tests/time_course.npz new file mode 100644 index 0000000..0551991 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/time_course.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/unbalanced.json b/tests/data/kenward_roger_and_f_tests/unbalanced.json new file mode 100644 index 0000000..8b5ca51 --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/unbalanced.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D02", "D02", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06", "D06", "D07", "D07", "D07", "D07", "D08", "D08", "D08", "D08", "D08"], "condition": ["ctrl", "ctrl", "stim", "stim", "stim", "stim", "stim", "stim", "ctrl", "stim", "ctrl", "stim", "ctrl", "ctrl", "stim", "stim", "stim", "ctrl", "ctrl", "ctrl", "ctrl", "stim", "ctrl", "stim"], "x": [0.464, 0.165, 0.564, -1.615, 0.607, -0.284, -0.186, 0.647, -0.428, 0.758, 0.429, 0.864, 0.764, 1.358, 0.449, 0.039, 1.0, -0.169, 0.669, -0.002, 1.982, 0.417, -1.215, -0.6]}, "formula": "~ condition + x", "random": "(1 | donor)", "variant": "none"} diff --git a/tests/data/kenward_roger_and_f_tests/unbalanced.npz b/tests/data/kenward_roger_and_f_tests/unbalanced.npz new file mode 100644 index 0000000..93486fe Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/unbalanced.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/weights.json b/tests/data/kenward_roger_and_f_tests/weights.json new file mode 100644 index 0000000..82c2e6b --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/weights.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D02", "D02", "D02", "D03", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [41.3, 41.3, 41.3, 54.3, 54.3, 54.3, 27.6, 27.6, 27.6, 61.3, 61.3, 61.3, 73.1, 73.1, 73.1, 61.2, 61.2, 61.2]}, "formula": "~ time + age", "random": "(1 | donor)", "variant": "average"} diff --git a/tests/data/kenward_roger_and_f_tests/weights.npz b/tests/data/kenward_roger_and_f_tests/weights.npz new file mode 100644 index 0000000..c0c1267 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/weights.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/weights_crossed.json b/tests/data/kenward_roger_and_f_tests/weights_crossed.json new file mode 100644 index 0000000..f9735d4 --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/weights_crossed.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D01", "D02", "D02", "D02", "D02", "D03", "D03", "D03", "D03", "D04", "D04", "D04", "D04", "D05", "D05", "D05", "D05", "D06", "D06", "D06", "D06", "D07", "D07", "D07", "D07", "D08", "D08", "D08", "D08"], "batch": ["B03", "B02", "B04", "B03", "B01", "B02", "B03", "B01", "B01", "B04", "B01", "B03", "B02", "B04", "B03", "B01", "B03", "B01", "B02", "B04", "B01", "B02", "B02", "B01", "B04", "B01", "B01", "B02", "B04", "B03", "B01", "B01"], "condition": ["ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim"], "x": [0.435, -0.893, 0.596, -0.539, 0.051, 1.5, 0.751, 2.569, -0.907, 0.676, -0.687, -0.186, -0.914, 1.364, -0.189, -0.129, 0.925, 0.182, -1.378, 1.289, 0.869, -0.401, -0.842, -2.832, 0.406, 0.157, -1.432, -0.614, 1.902, 0.766, 0.043, -1.843]}, "formula": "~ condition + x", "random": "(1 | donor) + (1 | batch)", "variant": "average"} diff --git a/tests/data/kenward_roger_and_f_tests/weights_crossed.npz b/tests/data/kenward_roger_and_f_tests/weights_crossed.npz new file mode 100644 index 0000000..24e9c9b Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/weights_crossed.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.json b/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.json new file mode 100644 index 0000000..d8e294f --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D01", "D02", "D02", "D02", "D02", "D03", "D03", "D03", "D03", "D04", "D04", "D04", "D04", "D05", "D05", "D05", "D05", "D06", "D06", "D06", "D06", "D07", "D07", "D07", "D07", "D08", "D08", "D08", "D08"], "batch": ["B03", "B02", "B04", "B03", "B01", "B02", "B03", "B01", "B01", "B04", "B01", "B03", "B02", "B04", "B03", "B01", "B03", "B01", "B02", "B04", "B01", "B02", "B02", "B01", "B04", "B01", "B01", "B02", "B04", "B03", "B01", "B01"], "condition": ["ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim", "ctrl", "ctrl", "stim", "stim"], "x": [0.435, -0.893, 0.596, -0.539, 0.051, 1.5, 0.751, 2.569, -0.907, 0.676, -0.687, -0.186, -0.914, 1.364, -0.189, -0.129, 0.925, 0.182, -1.378, 1.289, 0.869, -0.401, -0.842, -2.832, 0.406, 0.157, -1.432, -0.614, 1.902, 0.766, 0.043, -1.843]}, "formula": "~ condition + x", "random": "(1 | donor) + (1 | batch)", "variant": "precision"} diff --git a/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.npz b/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.npz new file mode 100644 index 0000000..c0f5c02 Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/weights_crossed_precision.npz differ diff --git a/tests/data/kenward_roger_and_f_tests/weights_precision.json b/tests/data/kenward_roger_and_f_tests/weights_precision.json new file mode 100644 index 0000000..c114e1e --- /dev/null +++ b/tests/data/kenward_roger_and_f_tests/weights_precision.json @@ -0,0 +1 @@ +{"samples": {"donor": ["D01", "D01", "D01", "D02", "D02", "D02", "D03", "D03", "D03", "D04", "D04", "D04", "D05", "D05", "D05", "D06", "D06", "D06"], "time": ["t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2", "t0", "t1", "t2"], "age": [41.3, 41.3, 41.3, 54.3, 54.3, 54.3, 27.6, 27.6, 27.6, 61.3, 61.3, 61.3, 73.1, 73.1, 73.1, 61.2, 61.2, 61.2]}, "formula": "~ time + age", "random": "(1 | donor)", "variant": "precision"} diff --git a/tests/data/kenward_roger_and_f_tests/weights_precision.npz b/tests/data/kenward_roger_and_f_tests/weights_precision.npz new file mode 100644 index 0000000..204bd5a Binary files /dev/null and b/tests/data/kenward_roger_and_f_tests/weights_precision.npz differ diff --git a/tests/test_dream.py b/tests/test_dream.py index 319ba81..1979382 100644 --- a/tests/test_dream.py +++ b/tests/test_dream.py @@ -186,10 +186,18 @@ def fixture_cube() -> tuple[Cube, Design]: def test_dream_voom_matches_variance_partition_on_every_cell_type(variant: str) -> None: data, _ = fixture("cube") cube, design = fixture_cube() + # The fixture is dream(ddf = "Satterthwaite"); the default, "adaptive", would + # use Kenward and Roger's method in these groups of fewer than 20 samples. result = dream_voom( - cube, design, prior_weights="n_cells", span=0.5, trend=variant == "trend" + cube, + design, + prior_weights="n_cells", + span=0.5, + trend=variant == "trend", + df_method="satterthwaite", ) assert result.method == "dream_voom" + assert result.summary["df_method"].tolist() == ["satterthwaite"] * 2 + [""] for group in ("B", "NK"): table = result.table(group, "stim_vs_ctrl") reference = { @@ -516,13 +524,18 @@ def test_dream_voom_records_its_settings_and_the_features_not_fitted() -> None: assert { "weights_fit": "ml", "fit": "reml", - "df_method": "satterthwaite", + "df_method": "adaptive", + "kenward_roger_weights": "average", + "f_tests": None, "trend": False, }.items() <= record["settings"].items() assert "quality_weights" not in record["settings"] group = record["groups"][0] assert group["features_not_fitted"] == 1 and group["features_not_converged"] == 0 assert group["voom_span"] == 1.0 + # 12 samples: dream's adaptive choice is Kenward and Roger's method. + assert group["df_method"] == "kenward_roger" + assert group["features_satterthwaite_fallback"] == 0 and group["f_test_df1"] == {} assert pickle.loads(pickle.dumps(result)).provenance == record # The design's rows in another order give the same tables. reversed_samples = {key: value[::-1] for key, value in samples.items()} diff --git a/tests/test_kenward_roger_and_f_tests.py b/tests/test_kenward_roger_and_f_tests.py new file mode 100644 index 0000000..ae6fe9a --- /dev/null +++ b/tests/test_kenward_roger_and_f_tests.py @@ -0,0 +1,728 @@ +"""Tests for Kenward-Roger's tests and the F-tests of mixed models, and dream's use of them. + +Reference values in ``tests/data/kenward_roger_and_f_tests`` come from pbkrtest 0.5.3 and +lmerTest 3.1-3 on lme4 1.1-38, evaluated at CytePool's REML estimates, and from +variancePartition 1.36.2's ``dream``, ``eBayes`` and ``topTable``, written by +``validation/kenward_roger_and_f_tests.py``. The other tests check closed forms of +balanced designs, where Kenward and Roger's degrees of freedom are those of the analysis +of variance and the adjustment vanishes, and invariants that a wrong implementation +breaks: reparametrised designs and equivalent contrasts, relabelled donors and +reordered samples, the scale of the weights, and an F-test of one contrast, which is the +square of its t-test. +""" + +import json +import pickle +from functools import cache +from pathlib import Path +from typing import Any + +import numpy as np +import pytest + +from cytepool._jit import thread_limit +from cytepool.cube import Cube, from_bulk +from cytepool.design import Design +from cytepool.differential_expression import dream_voom +from cytepool.stats._fixed_effect_inference import combined_df +from cytepool.stats.dream import moderated_fixed_effect_tests, voom_with_dream_weights +from cytepool.stats.mixed_models import ( + FixedEffectFTest, + MixedModelFit, + fit_mixed_model, + fixed_effect_f_test, + fixed_effect_tests, +) + +DATA = Path(__file__).parent / "data" / "kenward_roger_and_f_tests" +# Closed forms at given estimates. +EXACT = 1e-10 +# variancePartition's fits converge to about 1e-8 of the optimum (see the validation +# script), and its moderated tests inherit it. +CLOSE = 1e-6 +# lmerTest differentiates the deviance numerically for Satterthwaite's F-tests +# (validation/mixed_models.py). +SATTERTHWAITE = 1e-5 +CASES = ( + "time_course", + "unbalanced", + "crossed", + "nested", + "boundary", + "weights", + "weights_precision", + "weights_crossed", + "weights_crossed_precision", +) +CONTRASTS = { + "t2_vs_t1": [0, -1, 1, 0], + "(Intercept)": [1, 0, 0, 0], + "timet1": [0, 1, 0, 0], + "timet2": [0, 0, 1, 0], + "age": [0, 0, 0, 1], +} + + +@cache +def fixture(name: str) -> tuple[dict[str, np.ndarray], dict[str, Any]]: + with np.load(DATA / f"{name}.npz", allow_pickle=False) as data: + arrays = {key: data[key] for key in data.files} + return arrays, json.loads((DATA / f"{name}.json").read_text()) + + +@cache +def fitted(name: str) -> tuple[MixedModelFit, str]: + """The fit of a case's fixture genes, and the treatment of its weights.""" + data, meta = fixture(name) + design = Design(meta["samples"], meta["formula"], random=meta["random"]) + fit = fit_mixed_model(data["y"], design, weights=data.get("weights")) + return fit, "precision" if meta["variant"] == "precision" else "average" + + +def assert_relative( + actual: Any, desired: Any, tolerance: float, floor: Any = 0.0 +) -> None: + """|actual - desired| <= tolerance * max(|desired|, floor), NaN where desired is.""" + actual, desired = np.asarray(actual, float), np.asarray(desired, float) + np.testing.assert_array_equal(np.isnan(actual), np.isnan(desired)) + scale = np.fmax(np.abs(desired), floor) + finite = ~np.isnan(desired) + excess = np.abs(actual - desired)[finite] - tolerance * scale[finite] + assert np.all(excess <= 0.0), float(np.max(excess)) + + +def p_scale(p: np.ndarray) -> np.ndarray: + """max(p, 1e-300) max(1, |log p|): p-values compared on the scale of log p.""" + bounded = np.fmax(p, 1e-300) + return bounded * np.fmax(1.0, np.abs(np.log(bounded))) + + +# Against pbkrtest and lmerTest. + + +@pytest.mark.parametrize("name", CASES) +def test_kenward_roger_t_tests_match_pbkrtest_and_lmertest(name: str) -> None: + data, _ = fixture(name) + fit, variant = fitted(name) + tests = fixed_effect_tests( + fit, df_method="kenward_roger", kenward_roger_weights=variant + ) + adjusted = data["adjusted"] + scale = np.max(np.abs(adjusted), axis=(1, 2))[:, None, None] + assert_relative(tests.covariance, adjusted, EXACT, scale) + coef = data["coef"] # estimate, standard error, df, t, p + assert_relative(tests.estimate, coef[..., 0], EXACT, coef[..., 1]) + assert_relative(tests.standard_error, coef[..., 1], EXACT) + assert_relative(tests.t, coef[..., 3], EXACT, 1.0) + assert_relative(tests.df, coef[..., 2], EXACT) + assert_relative(tests.p_value, coef[..., 4], EXACT, p_scale(coef[..., 4])) + contrast = fixed_effect_tests( + fit, + data["contrast_weights"], + df_method="kenward_roger", + kenward_roger_weights=variant, + ) + expected = data["contrast"] + assert_relative(contrast.standard_error[:, 0], expected[:, 1], EXACT) + assert_relative(contrast.df[:, 0], expected[:, 2], EXACT) + assert_relative(contrast.t[:, 0], expected[:, 3], EXACT, 1.0) + assert tests.df_method == "kenward_roger" and tests.kenward_roger_weights == variant + + +@pytest.mark.parametrize("name", CASES) +def test_f_tests_match_krmodcomp_and_contestmd(name: str) -> None: + data, _ = fixture(name) + fit, variant = fitted(name) + for key in [key.split(" ", 1)[1] for key in data if key.startswith("f_matrix ")]: + matrix = data[f"f_matrix {key}"].T + # KRmodcomp's stat, ndf, ddf, F.scaling, p; contestMD's F, NumDF, DenDF, p by + # Kenward-Roger, then by Satterthwaite. + f = data[f"f {key}"] + kr = fixed_effect_f_test( + fit, matrix, df_method="kenward_roger", kenward_roger_weights=variant + ) + assert kr.df1 == 2 and np.all(f[:, 1] == 2) # 'pairs' spans 2 dimensions + assert_relative(kr.f_statistic, f[:, 0], EXACT) + assert_relative(kr.df2, f[:, 2], EXACT) + assert_relative(kr.scaling, f[:, 3], EXACT) + assert_relative(kr.p_value, f[:, 4], EXACT, p_scale(f[:, 4])) + satterthwaite = fixed_effect_f_test(fit, matrix) + assert_relative(satterthwaite.f_statistic, f[:, 9], EXACT) + np.testing.assert_array_equal(satterthwaite.scaling, 1.0) + assert_relative(satterthwaite.df2, f[:, 11], SATTERTHWAITE) + assert_relative( + satterthwaite.p_value, f[:, 12], SATTERTHWAITE, p_scale(f[:, 12]) + ) + + +def test_the_cases_cover_weights_two_random_effects_and_the_boundary() -> None: + # Without them the comparison would leave pbkrtest's treatment of weights, the + # residual's derivatives and the boundary components untested. + assert fixture("weights")[0]["weights"].std() > 0.3 + assert fitted("weights_crossed")[0].component_names == ("donor", "batch") + assert fitted("nested")[0].component_names == ("donor", "sample:donor") + assert np.count_nonzero(fitted("boundary")[0].boundary) >= 3 + # The two treatments of the weights differ on these data. + fit, _ = fitted("weights") + average = fixed_effect_tests(fit, df_method="kenward_roger") + precision = fixed_effect_tests( + fit, df_method="kenward_roger", kenward_roger_weights="precision" + ) + assert np.max(np.abs(average.df - precision.df) / precision.df) > 0.01 + + +# Against dream. + + +def study(name: str, method: str) -> tuple[Any, Any]: + data, samples = fixture(f"study_{name}") + design = Design(samples, "~ time + age", random="(1 | donor)") + voom = voom_with_dream_weights(data["counts"], design, span=0.5) + fit = fit_mixed_model(voom.log_cpm, design, weights=voom.weights) + tests = fixed_effect_tests(fit, CONTRASTS, df_method=method) # type: ignore[arg-type] + time = fixed_effect_tests( + fit, + {key: CONTRASTS[key] for key in ("timet1", "timet2")}, + df_method=method, # type: ignore[arg-type] + ) + return moderated_fixed_effect_tests(tests), moderated_fixed_effect_tests(time) + + +@pytest.mark.parametrize( + ("name", "ddf", "method"), + [ + ("n12", "Kenward-Roger", "kenward_roger"), + ("n12", "adaptive", "kenward_roger"), + ("n12", "Satterthwaite", "satterthwaite"), + ("n20", "adaptive", "satterthwaite"), + ], +) +def test_moderated_tests_match_dream(name: str, ddf: str, method: str) -> None: + """dream's adaptive choice is Kenward and Roger's method for 12 samples and + Satterthwaite's for 20, as the tests of each method reproduce.""" + data, _ = fixture(f"study_{name}") + moderated, time = study(name, method) + reference = { + key.split(" ", 1)[1]: value + for key, value in data.items() + if key.startswith(f"{ddf} ") + } + assert_relative(moderated.tests.standard_error, reference["se"], CLOSE) + assert_relative(moderated.tests.df, reference["df"], CLOSE) + assert_relative(moderated.t, reference["t"], CLOSE, 1.0) + assert_relative(moderated.lods, reference["lods"], CLOSE, 1.0) + assert_relative(moderated.z, reference["z"], CLOSE, 1.0) + assert_relative(moderated.df_total, reference["df_total"], CLOSE) + assert_relative(moderated.p_value, reference["p"], CLOSE, p_scale(reference["p"])) + # topTable's F; its p-values use the first gene's degrees of freedom for all. + assert_relative(time.f_statistic, reference["F"], CLOSE) + assert time.f_df1 == 2 and np.ptp(time.f_df2) > 1.0 + + +@cache +def cube_and_design() -> tuple[Cube, Design]: + data, samples = fixture("cube") + cube = Cube( + {"counts": data["counts"]}, + samples=samples, + groups={"group_id": ["small", "large"]}, + features={"feature_id": [f"gene{i}" for i in range(data["counts"].shape[2])]}, + n_cells=data["n_cells"], + ) + design = Design( + cube.samples, + "~ time + age", + random="(1 | donor)", + levels={"time": ["t0", "t1", "t2"]}, + contrasts={ + "t1_vs_t0": "timet1", + "t2_vs_t0": "timet2", + "t2_vs_t1": "timet2 - timet1", + }, + ) + return cube, design + + +@cache +def cube_result() -> Any: + cube, design = cube_and_design() + return dream_voom( + cube, design, span=0.5, f_tests={"time": ["t1_vs_t0", "t2_vs_t0"]} + ) + + +def test_dream_voom_by_default_chooses_the_degrees_of_freedom_as_dream_does() -> None: + """12 samples left in "small": Kenward and Roger's method; 24 in "large": + Satterthwaite's; the tables are dream's with ddf = "adaptive".""" + data, _ = fixture("cube") + result = cube_result() + assert result.summary["df_method"].tolist() == ["kenward_roger", "satterthwaite"] + assert result.summary["n_samples_used"].tolist() == [12, 24] + for group in ("small", "large"): + for name in ("t1_vs_t0", "t2_vs_t1"): + table = result.table(group, name) + reference = { + key.split(" ")[-1]: value + for key, value in data.items() + if key.startswith(f"{group} {name} ") + } + assert_relative(table["t"], reference["t"], CLOSE, 1.0) + assert_relative(table["df"], reference["df_total"], CLOSE) + assert_relative(table["b_statistic"], reference["lods"], CLOSE, 1.0) + assert_relative(table["z_statistic"], reference["z"], CLOSE, 1.0) + assert_relative( + table["p_value"], reference["p"], CLOSE, p_scale(reference["p"]) + ) + f_table = result.table(group, "time") + assert_relative(f_table["f_statistic"], data[f"{group} time F"], CLOSE) + + +def test_the_f_tables_of_dream_voom() -> None: + result = cube_result() + table = result.table("small", "time") + assert table.columns == ( + "feature_id", + "log2_fold_change_t1_vs_t0", + "log2_fold_change_t2_vs_t0", + "average_log2_cpm", + "f_statistic", + "df", + "p_value", + "adjusted_p_value", + "standardized_f_statistic", + ) + contrast = result.table("small", "t1_vs_t0") + np.testing.assert_array_equal( + table["log2_fold_change_t1_vs_t0"], contrast["log2_fold_change"] + ) + np.testing.assert_array_equal(table["feature_id"], contrast["feature_id"]) + # F.std is the F-statistic on (2, infinity) df, half a chi-squared on 2 df, with + # the test's p-value: exp(-F.std) is that p-value. + np.testing.assert_allclose( + np.exp(-table["standardized_f_statistic"]), table["p_value"], rtol=1e-12 + ) + assert result.f_tests == {"time": ("t1_vs_t0", "t2_vs_t0")} + groups = result.provenance["groups"] + assert [group["f_test_df1"] for group in groups] == [{"time": 2}, {"time": 2}] + assert [group["features_satterthwaite_fallback"] for group in groups] == [0, 0] + assert result.provenance["settings"]["f_tests"] == { + "time": ["t1_vs_t0", "t2_vs_t0"] + } + random = result.provenance["design"] + assert random["random"] == ["(1 | donor)"] and "donor" in random["samples"] + + +# Closed forms of balanced designs. + + +def test_a_balanced_one_way_design_has_kenward_rogers_closed_forms() -> None: + """In a balanced one-way layout the mean has the variance MSB / n on q - 1 degrees + of freedom exactly, and the adjustment vanishes: G_k maps the fixed effects' + columns into their span (Kenward and Roger 1997, section 5).""" + q, m = 6, 4 + rng = np.random.default_rng(31) + groups = np.repeat([f"G{i}" for i in range(q)], m) + design = Design({"group": groups}, "~ 1", random="(1 | group)") + y = 5.0 + rng.normal(size=(q, 5))[np.repeat(np.arange(q), m)] * [ + 0.6, + 1.0, + 2.0, + 3.0, + 5.0, + ] + y = y + rng.normal(size=(q * m, 5)) + fit = fit_mixed_model(y, design) + assert not fit.boundary.any() + kr = fixed_effect_tests(fit, df_method="kenward_roger") + satterthwaite = fixed_effect_tests(fit) + np.testing.assert_allclose(kr.df[:, 0], q - 1, rtol=EXACT) + np.testing.assert_allclose(kr.covariance, satterthwaite.covariance, rtol=EXACT) + means = y.reshape(q, m, -1).mean(axis=1) + msb = m * np.sum((means - y.mean(axis=0)) ** 2, axis=0) / (q - 1) + np.testing.assert_allclose(kr.standard_error[:, 0] ** 2, msb / (q * m), rtol=EXACT) + + +def test_a_randomised_block_design_has_the_classical_degrees_of_freedom() -> None: + """Donors as random blocks with each treatment once per donor: a difference of + treatments lies within donors, on (q - 1)(t - 1) df, and a treatment's mean has + Satterthwaite's combination of the donor and residual mean squares, as Kenward + and Roger's method gives; the adjustment vanishes. The F-test of the treatments + is the analysis of variance's, MS treatments / MSE on (t - 1, (q - 1)(t - 1)).""" + q, t = 7, 3 + rng = np.random.default_rng(5) + donors = np.repeat([f"D{i}" for i in range(q)], t) + treatment = np.tile(["a", "b", "c"], q) + design = Design( + {"donor": donors, "treatment": treatment}, + "~ 0 + treatment", + random="(1 | donor)", + ) + y = ( + np.array([4.0, 5.0, 4.5])[np.tile(np.arange(t), q)][:, None] + + rng.normal(0.0, 2.0, (q, 4))[np.repeat(np.arange(q), t)] + + rng.normal(0.0, 0.5, (q * t, 4)) + ) + fit = fit_mixed_model(y, design) + assert not fit.boundary.any() + tests = fixed_effect_tests( + fit, {"a": [1.0, 0, 0], "b_vs_a": [-1.0, 1.0, 0]}, df_method="kenward_roger" + ) + cells = y.reshape(q, t, -1) + donor_means, treatment_means = cells.mean(axis=1), cells.mean(axis=0) + grand = cells.mean(axis=(0, 1)) + msd = t * np.sum((donor_means - grand) ** 2, axis=0) / (q - 1) + residual = cells - donor_means[:, None] - treatment_means[None] + grand + mse = np.sum(residual**2, axis=(0, 1)) / ((q - 1) * (t - 1)) + combined = msd + (t - 1) * mse + expected = combined**2 / ( + msd**2 / (q - 1) + ((t - 1) * mse) ** 2 / ((q - 1) * (t - 1)) + ) + np.testing.assert_allclose(tests.df[:, 0], expected, rtol=EXACT) + np.testing.assert_allclose(tests.df[:, 1], (q - 1) * (t - 1), rtol=EXACT) + np.testing.assert_allclose(tests.standard_error[:, 1] ** 2, 2 * mse / q, rtol=EXACT) + f = fixed_effect_f_test( + fit, + {"b_vs_a": [-1.0, 1.0, 0], "c_vs_a": [-1.0, 0, 1.0]}, + df_method="kenward_roger", + ) + mst = q * np.sum((treatment_means - grand) ** 2, axis=0) / (t - 1) + np.testing.assert_allclose(f.f_statistic, mst / mse, rtol=EXACT) + np.testing.assert_allclose(f.df2, (q - 1) * (t - 1), rtol=EXACT) + np.testing.assert_allclose(f.scaling, 1.0, rtol=EXACT) + satterthwaite = fixed_effect_f_test( + fit, {"b_vs_a": [-1.0, 1.0, 0], "c_vs_a": [-1.0, 0, 1.0]} + ) + np.testing.assert_allclose(satterthwaite.f_statistic, mst / mse, rtol=EXACT) + np.testing.assert_allclose(satterthwaite.df2, (q - 1) * (t - 1), rtol=EXACT) + + +def test_the_denominator_df_of_an_f_test_combine_those_of_its_t_tests() -> None: + # Values derived by hand from the rule lmerTest documents for get_Fstat_ddf. + df = np.array( + [ + [5.0, 5.0], + [4.0, 12.0], + [2.0, 30.0], + [np.inf, 6.0], + [np.inf, np.inf], + [7.0, np.nan], + ] + ) + # [4, 12]: E = 4/2 + 12/10 = 3.2, 2 E / (E - 2) = 16/3; [inf, 6]: E = 1 + 3/2. + expected = [5.0, 16.0 / 3.0, 2.0, 10.0, np.inf, np.nan] + np.testing.assert_allclose(combined_df(df), expected, rtol=1e-15) + # One keeps its own, below 2 too, as lmerTest's contestMD does (found by running it). + np.testing.assert_allclose(combined_df(np.array([[9.0], [1.5]])), [9.0, 1.5]) + + +# Invariants. + + +def test_a_reparametrised_design_changes_only_the_coefficients() -> None: + data, meta = fixture("time_course") + fit, _ = fitted("time_course") + # Means of time points instead of a reference, and age in months: x_new = x t. + samples = dict( + meta["samples"], age_months=[12.0 * a for a in meta["samples"]["age"]] + ) + other = Design(samples, "~ 0 + time + age_months", random="(1 | donor)") + refit = fit_mixed_model(data["y"], other) + old = Design(meta["samples"], meta["formula"], random=meta["random"]).matrix + t = np.linalg.lstsq(old, other.matrix, rcond=None)[0] + first = fixed_effect_tests(fit, df_method="kenward_roger") + second = fixed_effect_tests(refit, df_method="kenward_roger") + inverse = np.linalg.inv(t) + transformed = np.einsum("ij,gjk,lk->gil", inverse, first.covariance, inverse) + scale = np.max(np.abs(second.covariance), axis=(1, 2))[:, None, None] + assert_relative(second.covariance, transformed, 1e-8, scale) + # The same contrasts in both parametrisations. + one = {"t2_vs_t1": [0, -1, 1, 0], "age": [0, 0, 0, 1]} + two = {"t2_vs_t1": [0, -1, 1, 0], "age": [0, 0, 0, 12.0]} + a = fixed_effect_tests(fit, one, df_method="kenward_roger") + b = fixed_effect_tests(refit, two, df_method="kenward_roger") + for name in ("estimate", "standard_error", "t", "df", "p_value"): + np.testing.assert_allclose(getattr(b, name), getattr(a, name), rtol=1e-8) + for method in ("kenward_roger", "satterthwaite"): + f1 = fixed_effect_f_test( + fit, {"t1": [0, 1, 0, 0], "t2": [0, 0, 1, 0]}, df_method=method + ) # type: ignore[arg-type] + f2 = fixed_effect_f_test( + refit, {"t1": [-1, 1, 0, 0], "t2": [-1, 0, 1, 0]}, df_method=method + ) # type: ignore[arg-type] + np.testing.assert_allclose(f2.f_statistic, f1.f_statistic, rtol=1e-8) + np.testing.assert_allclose(f2.p_value, f1.p_value, rtol=1e-7) + + +def test_equivalent_contrasts_give_the_same_kenward_roger_f_test() -> None: + fit, _ = fitted("weights") + spans = ( + {"t1": [0, 1, 0, 0], "t2": [0, 0, 1, 0]}, + {"t1": [0, 2, 0, 0], "t2_vs_t1": [0, -1, 1, 0]}, + {"t1": [0, 1, 0, 0], "t2": [0, 0, 1, 0], "t2_vs_t1": [0, -1, 1, 0]}, + ) + tests = [ + fixed_effect_f_test(fit, span, df_method="kenward_roger") for span in spans + ] + for other in tests[1:]: + assert other.df1 == 2 + for name in ("f_statistic", "df2", "scaling", "p_value"): + np.testing.assert_allclose( + getattr(other, name), getattr(tests[0], name), rtol=1e-9, err_msg=name + ) + # Satterthwaite's F-statistic is the same Wald statistic; its df depend on how the + # hypothesis is spanned, as lmerTest's do. + plain = [fixed_effect_f_test(fit, span) for span in spans] + for other in plain[1:]: + np.testing.assert_allclose(other.f_statistic, plain[0].f_statistic, rtol=1e-9) + + +@pytest.mark.parametrize("method", ["kenward_roger", "satterthwaite"]) +def test_an_f_test_of_one_contrast_is_the_square_of_its_t_test(method: str) -> None: + fit, _ = fitted("crossed") + contrast = {"stim_minus_x": [0.0, 1.0, -1.0]} + t = fixed_effect_tests(fit, contrast, df_method=method) # type: ignore[arg-type] + f = fixed_effect_f_test(fit, contrast, df_method=method) # type: ignore[arg-type] + assert f.df1 == 1 + np.testing.assert_allclose(f.f_statistic, t.t[:, 0] ** 2, rtol=1e-12) + np.testing.assert_allclose(f.df2, t.df[:, 0], rtol=1e-12) + np.testing.assert_allclose(f.p_value, t.p_value[:, 0], rtol=1e-10) + np.testing.assert_allclose(f.scaling, 1.0, rtol=1e-12) + # So is dream's moderated F-test of one contrast. + moderated = moderated_fixed_effect_tests(t) + assert moderated.f_df1 == 1 + np.testing.assert_allclose( + moderated.f_statistic, moderated.t[:, 0] ** 2, rtol=1e-12 + ) + np.testing.assert_allclose(moderated.f_df2, moderated.df_total[:, 0], rtol=1e-12) + np.testing.assert_allclose(moderated.f_p_value, moderated.p_value[:, 0], rtol=1e-10) + + +@pytest.mark.parametrize("method", ["kenward_roger", "satterthwaite"]) +def test_an_f_test_of_one_contrast_keeps_degrees_of_freedom_below_two( + method: str, +) -> None: + """A difference between donors, two against one, has about one degree of freedom. + Its F-test is still its t-test, as lmerTest's contestMD gives it (found by running + it), not the 2 that the rule for several contrasts gives when one is at most 2.""" + rng = np.random.default_rng(3) + donors = np.repeat(["A", "B", "C"], 3) + group = np.repeat(["x", "x", "y"], 3) + design = Design({"donor": donors, "group": group}, "~ group", random="(1 | donor)") + y = ( + rng.normal(size=(9, 6)) + + rng.normal(0.0, 2.0, (3, 6))[np.repeat(np.arange(3), 3)] + + (group == "y")[:, np.newaxis] + ) + fit = fit_mixed_model(y, design) + contrast = {"y_vs_x": [0.0, 1.0]} + t = fixed_effect_tests(fit, contrast, df_method=method) # type: ignore[arg-type] + f = fixed_effect_f_test(fit, contrast, df_method=method) # type: ignore[arg-type] + assert np.count_nonzero(t.df[:, 0] < 2.0) >= 3 + np.testing.assert_allclose(f.df2, t.df[:, 0], rtol=1e-12) + np.testing.assert_allclose(f.p_value, t.p_value[:, 0], rtol=1e-10) + + +def test_relabelled_donors_and_reordered_samples_change_nothing() -> None: + data, meta = fixture("weights_crossed") + fit, _ = fitted("weights_crossed") + order = np.random.default_rng(2).permutation(data["y"].shape[0]) + samples = {key: np.asarray(value)[order] for key, value in meta["samples"].items()} + samples["donor"] = np.array(["donor " + str(d)[::-1] for d in samples["donor"]]) + samples["batch"] = np.array(["x" + str(b) for b in samples["batch"]]) + design = Design(samples, meta["formula"], random=meta["random"]) + refit = fit_mixed_model(data["y"][order], design, weights=data["weights"][order]) + for variant in ("average", "precision"): + first = fixed_effect_tests( + fit, + df_method="kenward_roger", + kenward_roger_weights=variant, # type: ignore[arg-type] + ) + second = fixed_effect_tests( + refit, + df_method="kenward_roger", + kenward_roger_weights=variant, # type: ignore[arg-type] + ) + for name in ("standard_error", "df", "p_value"): + # Sums taken in another order differ in rounding. + np.testing.assert_allclose( + getattr(second, name), getattr(first, name), rtol=1e-8, err_msg=name + ) + + +def test_the_scale_of_the_weights_changes_nothing_and_equal_weights_are_no_weights() -> ( + None +): + data, meta = fixture("weights") + design = Design(meta["samples"], meta["formula"], random=meta["random"]) + fit = fit_mixed_model(data["y"], design, weights=data["weights"]) + larger = fit_mixed_model(data["y"], design, weights=7.0 * data["weights"]) + for variant in ("average", "precision"): + a = fixed_effect_tests( + fit, df_method="kenward_roger", kenward_roger_weights=variant + ) # type: ignore[arg-type] + b = fixed_effect_tests( + larger, df_method="kenward_roger", kenward_roger_weights=variant + ) # type: ignore[arg-type] + np.testing.assert_allclose(b.df, a.df, rtol=1e-8) + np.testing.assert_allclose(b.standard_error, a.standard_error, rtol=1e-8) + plain = fit_mixed_model(data["y"], design) + equal = fit_mixed_model(data["y"], design, weights=np.full(data["y"].shape[0], 3.0)) + reference = fixed_effect_tests(plain, df_method="kenward_roger") + for variant in ("average", "precision"): + tests = fixed_effect_tests( + equal, df_method="kenward_roger", kenward_roger_weights=variant + ) # type: ignore[arg-type] + np.testing.assert_allclose(tests.df, reference.df, rtol=1e-8) + np.testing.assert_allclose( + tests.standard_error, reference.standard_error, rtol=1e-8 + ) + + +def test_the_number_of_threads_changes_nothing() -> None: + # The features are computed in runs on Python threads, each by the same + # operations. + fit, _ = fitted("weights_crossed") + with thread_limit(1): + serial = fixed_effect_tests(fit, df_method="kenward_roger") + with thread_limit(3): + threaded = fixed_effect_tests(fit, df_method="kenward_roger") + for name in ("standard_error", "df", "covariance"): + np.testing.assert_array_equal(getattr(threaded, name), getattr(serial, name)) + + +# dream_voom: its choices, its records and its errors. + + +def ill_conditioned_cube() -> tuple[Cube, Design]: + """15 samples whose age, 32 years on average and about 1 apart, leaves the + condition numbers of the adjusted covariance matrices near 1e6.""" + rng = np.random.default_rng(41) + donor = np.repeat(np.arange(5), 3) + samples = { + "sample_id": [f"S{i}" for i in range(15)], + "donor": [f"D{d}" for d in donor], + "time": ["t0", "t1", "t2"] * 5, + "age": (32.0 + np.round(rng.normal(size=15), 3)).tolist(), + } + log_mean = rng.normal(3.0, 1.0, 40) + rng.normal(0.0, 0.4, (5, 40))[donor] + counts = rng.negative_binomial( + 10, 10 / (10 + np.exp(log_mean) * rng.uniform(0.5, 2.0, (15, 1))) + ) + cube = from_bulk(counts, samples=samples) + design = Design( + cube.samples, "~ time + age", random="(1 | donor)", contrasts={"t1": "timet1"} + ) + return cube, design + + +def test_dream_voom_falls_back_to_satterthwaite_where_dream_does() -> None: + """Features whose adjusted covariance has a condition number above 1e6 take + Satterthwaite's tests; the others Kenward and Roger's.""" + cube, design = ill_conditioned_cube() + kr = dream_voom(cube, design, df_method="kenward_roger", expression_filter=False) + satterthwaite = dream_voom( + cube, design, df_method="satterthwaite", expression_filter=False + ) + group = kr.provenance["groups"][0] + first, second = kr.table("bulk", "t1"), satterthwaite.table("bulk", "t1") + same = np.isclose(first["df"], second["df"], rtol=1e-12) + assert 0 < group["features_satterthwaite_fallback"] < len(same) + assert np.count_nonzero(same) == group["features_satterthwaite_fallback"] + np.testing.assert_allclose( + np.asarray(first["t"])[same], np.asarray(second["t"])[same], rtol=1e-10 + ) + # dream_voom with its default is Kenward and Roger's for 15 samples. + default = dream_voom(cube, design, expression_filter=False) + np.testing.assert_array_equal(default.table("bulk", "t1")["df"], first["df"]) + assert default.summary["df_method"][0] == "kenward_roger" + + +def test_an_f_test_with_a_contrast_a_group_cannot_estimate_is_not_tested() -> None: + cube, design = cube_and_design() + cells = np.asarray(cube.n_cells).copy() + cells[np.asarray(cube.samples["time"]) == "t2", 0] = 5 # no t2 left in "small" + trimmed = Cube( + {"counts": np.asarray(cube.layer("counts"))}, + samples=cube.samples.to_dict(), + groups={"group_id": ["small", "large"]}, + features=cube.features.to_dict(), + n_cells=cells, + ) + result = dream_voom( + trimmed, + design, + contrasts={"t1_vs_t0": "timet1", "t2_vs_t0": "timet2"}, + f_tests={"time": ["t1_vs_t0", "t2_vs_t0"]}, + span=0.5, + ) + reason = result.provenance["groups"][0]["not_estimable"]["time"] + assert "cannot estimate its contrast 't2_vs_t0'" in reason + assert result.summary["not_estimable"][0] == "t2_vs_t0, time" + with pytest.raises(KeyError, match="cannot estimate its contrast"): + result.table("small", "time") + assert ("large", "time") in result.tables + + +def test_invalid_arguments_are_refused_with_what_to_do() -> None: + cube, design = cube_and_design() + with pytest.raises(ValueError, match="df_method must be one of 'adaptive'"): + dream_voom(cube, design, df_method="lme4") # type: ignore[arg-type] + with pytest.raises(TypeError, match="f_tests must be a mapping"): + dream_voom(cube, design, f_tests=["t1_vs_t0"]) # type: ignore[arg-type] + with pytest.raises(TypeError, match=r"f_tests\['time'\] must be a list"): + dream_voom(cube, design, f_tests={"time": "t1_vs_t0"}) + with pytest.raises(ValueError, match="must name at least one contrast"): + dream_voom(cube, design, f_tests={"time": []}) + with pytest.raises(ValueError, match="'t1_vs_t0' more than once"): + dream_voom(cube, design, f_tests={"time": ["t1_vs_t0", "t1_vs_t0"]}) + with pytest.raises(ValueError, match="names 't3_vs_t0', which is not among"): + dream_voom(cube, design, f_tests={"time": ["t1_vs_t0", "t3_vs_t0"]}) + with pytest.raises(ValueError, match="has the name of a contrast"): + dream_voom(cube, design, f_tests={"t1_vs_t0": ["t1_vs_t0", "t2_vs_t0"]}) + fit, _ = fitted("time_course") + ml = fit_mixed_model( + fixture("time_course")[0]["y"], + Design( + fixture("time_course")[1]["samples"], "~ time + age", random="(1 | donor)" + ), + method="ml", + ) + with pytest.raises(ValueError, match="defined for REML fits"): + fixed_effect_tests(ml, df_method="kenward_roger") + with pytest.raises(ValueError, match="defined for REML fits"): + fixed_effect_f_test(ml, np.eye(4)[:, 1:3], df_method="kenward_roger") + with pytest.raises(ValueError, match="the contrasts are all zero"): + fixed_effect_f_test(fit, np.zeros((4, 2))) + with pytest.raises(TypeError, match="fit must be a MixedModelFit"): + fixed_effect_f_test(fit.design, np.eye(4)) # type: ignore[arg-type] + + +def test_features_without_a_fit_have_nan_kenward_roger_tests() -> None: + data, meta = fixture("time_course") + design = Design(meta["samples"], meta["formula"], random=meta["random"]) + y = np.column_stack([data["y"][:, 0], np.full(design.n_samples, 4.0)]) + fit = fit_mixed_model(y, design) + tests = fixed_effect_tests(fit, df_method="kenward_roger") + f = fixed_effect_f_test(fit, np.eye(4)[:, 1:3], df_method="kenward_roger") + assert np.isfinite(tests.df[0]).all() and np.isnan(tests.df[1]).all() + assert np.isnan(tests.standard_error[1]).all() and np.isnan(f.p_value[1]) + assert np.isfinite(f.p_value[0]) + + +def test_the_f_test_is_a_read_only_record() -> None: + fit, _ = fitted("unbalanced") + f = fixed_effect_f_test( + fit, {"stim": [0, 1, 0], "x": [0, 0, 1]}, df_method="kenward_roger" + ) + assert isinstance(f, FixedEffectFTest) and f.fit is fit + assert f.contrast_names == ("stim", "x") and f.kenward_roger_weights == "average" + assert repr(f) == ( + "FixedEffectFTest(n_features=6, contrast_names=('stim', 'x'), df1=2, " + "df_method='kenward_roger')" + ) + for array in (f.f_statistic, f.df2, f.p_value, f.scaling, f.contrasts): + with pytest.raises(ValueError, match="read-only"): + array[0] = 0 + copy = pickle.loads(pickle.dumps(f)) + np.testing.assert_array_equal(copy.p_value, f.p_value) + assert fixed_effect_f_test(fit, {"stim": [0, 1, 0]}).kenward_roger_weights is None diff --git a/tests/test_masked_inputs.py b/tests/test_masked_inputs.py index 2343b49..2d11ba6 100644 --- a/tests/test_masked_inputs.py +++ b/tests/test_masked_inputs.py @@ -274,14 +274,17 @@ def _weighted_by_cells(context: Context) -> dict[str, Any]: def _dream(context: Context) -> dict[str, Any]: """As _differential_expression, with a random donor intercept and no quality - weights, which dream_voom does not estimate.""" + weights, which dream_voom does not estimate, and an F-test.""" arguments = _differential_expression(context) del arguments["quality_weights"] return arguments | { "design": Design( SAMPLE_TABLE, "~ condition", random="(1 | donor)", levels=LEVELS ), + "f_tests": {"both": ["stim"]}, "span": 0.6, + "df_method": "kenward_roger", + "kenward_roger_weights": "precision", "n_threads": 1, } @@ -804,7 +807,19 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "names": ["u", "v"], }, {"contrasts": None, "names": None}, - options=("df_method",), + options=("df_method", "kenward_roger_weights"), + objects=("fit",), + ), + "cytepool.stats.mixed_models.fixed_effect_f_test": Entry( + lambda c: { + "fit": fit_mixed_model(**_mixed_model(c)), + "contrasts": np.eye(2), + "names": ["u", "v"], + "df_method": "kenward_roger", + "kenward_roger_weights": "precision", + }, + {"contrasts": None, "names": None}, + options=("df_method", "kenward_roger_weights"), objects=("fit",), ), # dream: voom weights from mixed models, and their moderated tests. @@ -1495,6 +1510,7 @@ def _forest_arguments(context: Context) -> dict[str, Any]: _dream, { "contrasts": None, + "f_tests": None, "groups": None, "min_cells": None, "expression_filter": None, @@ -1504,12 +1520,14 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "min_counts", "normalization", "span", + "df_method", + "kenward_roger_weights", "robust", "trend", "adjust_across_groups", ), objects=("cube", "design"), - keys={"contrasts": "stim", "expression_filter": "min_count"}, + keys={"contrasts": "stim", "f_tests": "both", "expression_filter": "min_count"}, ), "cytepool.differential_expression.dream_voom[prior_weights]": Entry( _dream_by_cells, {}, options=("prior_weights",) @@ -2001,6 +2019,7 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "cytepool.stats.deseq2.DESeq2Dispersions", "cytepool.stats.deseq2.DESeq2Fit", "cytepool.stats.deseq2.DESeq2Results", + "cytepool.stats.mixed_models.FixedEffectFTest", "cytepool.stats.mixed_models.FixedEffectTests", "cytepool.stats.mixed_models.MixedModelFit", "cytepool.meta_analysis.MetaAnalysisFit", diff --git a/tests/test_mixed_models.py b/tests/test_mixed_models.py index 144541c..6e21231 100644 --- a/tests/test_mixed_models.py +++ b/tests/test_mixed_models.py @@ -794,8 +794,12 @@ def test_invalid_tests_are_refused() -> None: fixed_effect_tests( fit_linear_model(reference("one_intercept")["y"], fit.design) ) - with pytest.raises(ValueError, match="Kenward-Roger's degrees of freedom are not"): - fixed_effect_tests(fit, df_method="kenward_roger") # type: ignore[arg-type] + with pytest.raises( + ValueError, match="df_method must be one of 'satterthwaite', 'kenward_roger'" + ): + fixed_effect_tests(fit, df_method="containment") # type: ignore[arg-type] + with pytest.raises(TypeError, match="kenward_roger_weights must be one of"): + fixed_effect_tests(fit, kenward_roger_weights=1) # type: ignore[arg-type] with pytest.raises(ValueError, match="names are given only with contrasts"): fixed_effect_tests(fit, names=["a"]) with pytest.raises(ValueError, match=r"must weigh exactly the fit's coefficients"): diff --git a/tests/test_private_names.py b/tests/test_private_names.py index 5867847..31fbb6f 100644 --- a/tests/test_private_names.py +++ b/tests/test_private_names.py @@ -63,6 +63,10 @@ ("cytepool.stats.distributions", "_f_tail_quantile"): _F_DISTRIBUTION, ("cytepool.stats.distributions", "_f_tails"): _F_DISTRIBUTION, ("cytepool.stats.dream", "_dream_weights"): _VOOM_RESULT, + ("cytepool.stats.dream", "_moderated_tests"): ( + "builds the public ModeratedFixedEffectTests of cytepool.stats.dream from the " + "parts of tests that dream_voom chooses feature by feature" + ), ("cytepool.stats.mixed_models", "_fit"): ( "builds the public MixedModelFit of cytepool.stats.mixed_models" ), diff --git a/validation/dream.py b/validation/dream.py index b658f71..a87065d 100644 --- a/validation/dream.py +++ b/validation/dream.py @@ -3,8 +3,10 @@ Compares ``cytepool.stats.dream`` (``voom_with_dream_weights``, ``moderated_fixed_effect_tests``), the mixed models behind them (``fit_mixed_model`` by maximum likelihood and REML, with fitted values and approximate residual degrees -of freedom) and ``dream_voom`` with variancePartition's ``voomWithDreamWeights``, -``dream(ddf = "Satterthwaite")``, ``eBayes`` and ``topTable``, on counts simulated from +of freedom) and ``dream_voom(df_method="satterthwaite")`` with variancePartition's +``voomWithDreamWeights``, ``dream(ddf = "Satterthwaite")``, ``eBayes`` and ``topTable`` +(``validation/kenward_roger_and_f_tests.py`` checks the other choices of degrees of +freedom and the F-tests), on counts simulated from fixed seeds: a random donor intercept with a condition, a covariate and a contrast, and ``eBayes`` with ``robust = TRUE`` and ``trend = TRUE``; crossed random intercepts with sample prior weights; unbalanced donors, some with a single sample; and @@ -612,8 +614,15 @@ def compare_cube( results: Results, ) -> None: for variant in ("default", "trend"): + # dream(ddf = "Satterthwaite"), which the R script runs; dream_voom's default + # is dream's adaptive choice (validation/kenward_roger_and_f_tests.py). result = dream_voom( - cube, design, prior_weights="n_cells", span=0.5, trend=variant == "trend" + cube, + design, + prior_weights="n_cells", + span=0.5, + trend=variant == "trend", + df_method="satterthwaite", ) for group, given in inputs.items(): ref, label = ( diff --git a/validation/kenward_roger_and_f_tests.py b/validation/kenward_roger_and_f_tests.py new file mode 100644 index 0000000..97a5f12 --- /dev/null +++ b/validation/kenward_roger_and_f_tests.py @@ -0,0 +1,1052 @@ +"""Check Kenward-Roger's tests and the F-tests of mixed models, and dream's use of them. + +Compares, gene by gene, on data simulated from fixed seeds: + +- ``fixed_effect_tests(df_method="kenward_roger")`` with lmerTest's ``summary`` and + ``contest1D`` (``ddf = "Kenward-Roger"``) and pbkrtest's ``vcovAdj``; +- ``fixed_effect_f_test`` with pbkrtest's ``KRmodcomp`` and lmerTest's ``contestMD`` + (``ddf = "Kenward-Roger"`` and ``"Satterthwaite"``), for two contrasts and for three + that span two dimensions; +- with observation weights, ``kenward_roger_weights="average"`` with pbkrtest on lme4 + fits whose weights are divided by their mean, as dream fits them, and + ``"precision"`` with pbkrtest on the same model without weights, the values, fixed + effects and random effects multiplied by the square roots of the weights (y* = W^1/2 + y, X* = W^1/2 X, ``(0 + sw | donor)``): the weighted model itself, whose covariance + pbkrtest then builds from its matrices; +- variancePartition's ``dream`` with ``ddf = "Kenward-Roger"``, ``"adaptive"`` and + ``"Satterthwaite"`` after ``voomWithDreamWeights``, with ``eBayes`` and ``topTable``, + including its F-tests of several coefficients, against + ``voom_with_dream_weights``, ``fixed_effect_tests`` and + ``moderated_fixed_effect_tests``, on studies of 12 to 24 samples, the adaptive + choice at 19 and 20 samples, and genes whose adjusted covariance is ill conditioned; +- ``dream_voom`` with its default, ``df_method="adaptive"``, and an F-test, on a cube + whose cell types have 12 and 24 samples, against ``dream`` on each cell type's + samples, features and library sizes. + +Designs: one random intercept with a three-level time course and a donor covariate; +unbalanced donors, some with one sample; crossed and nested random intercepts; a design +simulated without donor variance, so that most fits lie on the boundary; and +observation weights like voom's with one and two random intercepts. + +At given estimates: lme4 and lmerTest are evaluated once at CytePool's REML estimates +(passed as lme4's starting values with one evaluation allowed), where Kenward-Roger's +adjusted covariance, degrees of freedom and F-tests are closed forms, compared to the +exact tier (1e-10); lmerTest's Satterthwaite F-tests differentiate the deviance +numerically and are compared to 1e-5, as in ``validation/mixed_models.py``. At the +optimum, lme4 fits each gene from its own start with Powell's BOBYQA (minqa) and a trust +region shrunk to 1e-10, and the tests are compared to the close tier (1e-6); dream's +moderated tests too, the fits made with that optimiser. + +Deliberate differences, printed as such: dream's p-values of an F-test use the total +degrees of freedom of the first coefficient of the first gene for every gene, and its +``F.std`` the number of coefficients rather than their rank; CytePool gives each gene +its own degrees of freedom. And dream decides from R's estimate of the condition number +(``kappa``) which genes fall back to Satterthwaite's method, CytePool from the exact +one: the two can decide differently within about 2% of the bound, 1e6. + +The script also reports how the default, dream's adaptive choice, changes the results +of studies with fewer than 20 samples, which ``dream_voom`` tested with +Satterthwaite's method before it. + +Requirements: + Python 3.13 with numpy and numba, run from the repository root so that + ``cytepool`` imports from the working tree. + R 4.4 with variancePartition 1.36.2, pbkrtest, lmerTest, lme4 1.1-x, minqa and + jsonlite. Set ``RSCRIPT`` to the Rscript executable when it is not on PATH. + +Usage: + PYTHONPATH=. python validation/kenward_roger_and_f_tests.py + PYTHONPATH=. python validation/kenward_roger_and_f_tests.py --write-fixtures +""" + +import argparse +import json +import os +import subprocess +import tempfile +from pathlib import Path +from typing import Any + +import numpy as np + +from cytepool.cube import Cube +from cytepool.design import Design +from cytepool.differential_expression import dream_voom +from cytepool.normalization import filter_by_expression, library_sizes, tmm_factors +from cytepool.stats.distributions import f_sf +from cytepool.stats.dream import moderated_fixed_effect_tests, voom_with_dream_weights +from cytepool.stats.mixed_models import ( + fit_mixed_model, + fixed_effect_f_test, + fixed_effect_tests, +) + +FIXTURES = ( + Path(__file__).resolve().parents[1] / "tests" / "data" / "kenward_roger_and_f_tests" +) +EXACT, CLOSE = 1e-10, 1e-6 +SATTERTHWAITE = 1e-5 # lmerTest's numerical derivatives (validation/mixed_models.py) +N_GENES, N_FIXTURE = 24, 6 +RELATIVE = "max rel. difference" +SCALED = "max diff. / max(abs, 1)" +P_VALUE = "max diff. / (p max(1, abs log p))" + +R_SCRIPT = r""" +suppressPackageStartupMessages({library(variancePartition); library(lmerTest); library(pbkrtest); library(jsonlite)}) +args <- commandArgs(trailingOnly = TRUE) +input <- read_json(args[[1]], simplifyVector = TRUE) +quiet <- function(expr) suppressWarnings(suppressMessages(expr)) +given <- lmerControl(optimizer = "nloptwrap", calc.derivs = FALSE, check.conv.singular = "ignore", + restart_edge = FALSE, boundary.tol = 0, optCtrl = list(maxeval = 1)) +tight_bobyqa <- function(par, fn, lower, upper, control) { + res <- minqa::bobyqa(par, fn, lower = lower, upper = upper, + control = list(rhobeg = 0.2, rhoend = 1e-10, maxfun = 1e5)) + list(par = res$par, fval = res$fval, conv = 0, message = res$msg) +} +tight <- lmerControl(optimizer = "tight_bobyqa", calc.derivs = FALSE, check.conv.singular = "ignore") +groups_of <- function(fit) sub("\\..*$", "", names(getME(fit, "theta"))) +out <- list(versions = as.list(c(R = as.character(getRversion()), sapply( + c("variancePartition", "limma", "lme4", "lmerTest", "pbkrtest", "minqa"), + function(p) as.character(packageVersion(p)))))) +for (name in names(input$fixed)) { + case <- input$fixed[[name]] + d <- as.data.frame(case$samples, stringsAsFactors = FALSE) + G <- nrow(case$y); p <- ncol(case$x) + res <- list(adjusted = list(), coef = list(), contrast = list(), f = list(), optimum = list()) + for (g in seq_len(G)) { + w <- if (is.null(case$w)) rep(1, nrow(d)) else case$w[g, ] + if (case$variant == "precision") { # the weighted model, its rows scaled + d$y <- sqrt(w) * case$y[g, ]; d$X <- sqrt(w) * case$x; d$sw <- sqrt(w); d$ww <- 1 + } else { + d$y <- case$y[g, ]; d$X <- case$x; d$sw <- 1; d$ww <- w / mean(w) + } + form <- as.formula(paste("y ~ 0 + X +", case$random)) + probe <- quiet(lmer(form, d, weights = ww, REML = TRUE, control = given)) + theta <- case$theta[g, match(groups_of(probe), case$components)] + fit <- quiet(lmerTest::lmer(form, d, weights = ww, REML = TRUE, start = list(theta = theta), control = given)) + res$adjusted[[g]] <- as.matrix(vcovAdj(fit)) + res$coef[[g]] <- summary(fit, ddf = "Kenward-Roger")$coefficients[, c("Estimate", "Std. Error", "df", "t value", "Pr(>|t|)")] + res$contrast[[g]] <- unlist(contest1D(fit, case$contrast, ddf = "Kenward-Roger")[1, ]) + res$f[[g]] <- lapply(case$f, function(L) { + L <- matrix(L, ncol = p) + c(unlist(KRmodcomp(fit, L)$test[1, c("stat", "ndf", "ddf", "F.scaling", "p.value")]), + unlist(contestMD(fit, L, ddf = "Kenward-Roger")[1, c("F value", "NumDF", "DenDF", "Pr(>F)")]), + unlist(contestMD(fit, L, ddf = "Satterthwaite")[1, c("F value", "NumDF", "DenDF", "Pr(>F)")])) + }) + best <- quiet(lmerTest::lmer(form, d, weights = ww, REML = TRUE, control = tight)) + res$optimum[[g]] <- summary(best, ddf = "Kenward-Roger")$coefficients[, c("Std. Error", "df")] + } + out$fixed[[name]] <- res +} +for (name in names(input$dream)) { + case <- input$dream[[name]] + counts <- case$counts; data <- as.data.frame(case$samples, stringsAsFactors = FALSE) + colnames(counts) <- rownames(data) <- paste0("s", seq_len(ncol(counts))) + rownames(counts) <- paste0("g", seq_len(nrow(counts))) + form <- as.formula(case$formula) + v <- quiet(voomWithDreamWeights(counts, form, data, lib.size = case$lib_size, span = 0.5, control = tight)) + L <- if (length(case$contrasts)) quiet(makeContrastsDream(form, data, contrasts = unlist(case$contrasts))) else NULL + r <- list(weights = v$weights) + for (ddf in case$ddf) { + fit <- quiet(dream(v, form, data, L = L, ddf = ddf, control = tight)) + e <- eBayes(fit) + s <- list(names = colnames(fit$coefficients), df = fit$df.residual, se = fit$stdev.unscaled * fit$sigma, + t = e$t, df_total = e$df.total, p = e$p.value, lods = e$lods, + z = sapply(colnames(fit$coefficients), function(k) topTable(e, coef = k, number = Inf, sort.by = "none")$z.std)) + for (f in names(case$f_tests)) { + top <- topTable(e, coef = case$f_tests[[f]], number = Inf, sort.by = "none") + s$f[[f]] <- list(F = top$F, p = top$P.Value) + } + r[[ddf]] <- s + } + if (isTRUE(case$condition)) { # the condition numbers that decide dream's fallback + r$kappa <- sapply(seq_len(nrow(v$E)), function(g) { + d <- data; d$yy <- v$E[g, ]; d$ww <- v$weights[g, ] / mean(v$weights[g, ]) + m <- quiet(lmerTest::lmer(update(form, yy ~ .), d, weights = ww, control = tight)) + kappa(as.matrix(vcovAdj(m))) + }) + } + out$dream[[name]] <- r +} +write_json(out, args[[2]], digits = NA, matrix = "rowmajor", auto_unbox = TRUE) +""" + + +class Results: + def __init__(self) -> None: + self.rows: list[tuple[str, ...]] = [] + self.failures = 0 + + def check( + self, check: str, reference: str, metric: str, value: float, tolerance: float + ) -> None: + passed = bool(value <= tolerance) + self.failures += not passed + result = "pass" if passed else "FAIL" + self.rows.append( + (check, reference, metric, f"{value:.1e}", f"{tolerance:.0e}", result) + ) + + def other( + self, + check: str, + reference: str, + metric: str, + value: str, + expected: bool | None = None, + ) -> None: + """An information row, or a deliberate difference that must be as expected.""" + self.failures += expected is False + result = ( + "info" + if expected is None + else "deliberate difference" + if expected + else "FAIL" + ) + self.rows.append((check, reference, metric, value, "n/a", result)) + + def print(self) -> None: + print("| check | reference | metric | value | tolerance | result |") + print("| --- | --- | --- | --- | --- | --- |") + for row in self.rows: + print("| " + " | ".join(row) + " |") + + +def relative(actual: Any, desired: Any, floor: Any = 0.0) -> float: + """Largest difference over max(|desired|, floor); infinite where NaN differs.""" + actual, desired = np.asarray(actual, float), np.asarray(desired, float) + if actual.shape != desired.shape or not np.array_equal( + np.isnan(actual), np.isnan(desired) + ): + return np.inf + finite = np.isfinite(desired) + scale = np.broadcast_to(np.fmax(np.abs(desired), floor), desired.shape)[finite] + difference = np.abs(actual[finite] - desired[finite]) + with np.errstate(invalid="ignore", divide="ignore"): + ratios = np.where(difference == 0.0, 0.0, difference / scale) + return float(np.max(ratios, initial=0.0)) + + +def p_value_difference(actual: Any, desired: Any) -> float: + """Difference relative to max(p, 1e-300) max(1, |log p|), the scale of log p.""" + p = np.fmax(np.asarray(desired, float), 1e-300) + return relative(actual, desired, p * np.fmax(1.0, np.abs(np.log(p)))) + + +def labels(prefix: str, codes: Any) -> list[str]: + return [f"{prefix}{code + 1:02d}" for code in np.asarray(codes)] + + +# The tests of the fixed effects: pbkrtest and lmerTest at CytePool's estimates. + + +def simulate(rng: Any, design: Design, n_genes: int, random_sd: float) -> np.ndarray: + """Values with gene-wise effects, residual SDs and random intercepts, whose SD is + up to ``random_sd`` times the residual one and zero for a fifth of the genes.""" + x = design.matrix + effects = rng.normal(0.0, 1.0, (x.shape[1], n_genes)) * ( + rng.random((x.shape[1], n_genes)) < 0.6 + ) + effects[0] = rng.uniform(2.0, 10.0, n_genes) + sd = 0.3 * np.sqrt(4.0 / rng.chisquare(4.0, n_genes)) + y = x @ effects + for effect in design.random: + joined = [ + "\t".join(map(str, row)) + for row in zip(*(design.samples[g] for g in effect.grouping)) + ] + codes = np.unique(joined, return_inverse=True)[1] + scale = rng.uniform(0.0, random_sd, n_genes) * (rng.random(n_genes) > 0.2) + y = y + rng.normal(0.0, 1.0, (codes.max() + 1, n_genes))[codes] * scale * sd + return y + rng.normal(0.0, 1.0, y.shape) * sd + + +def fixed_cases(seed: int, n_genes: int) -> dict[str, dict[str, Any]]: + """The designs, values, weights, contrasts and F-tests of the fixed effects' tests.""" + rng = np.random.default_rng(seed) + donor = np.repeat(np.arange(6), 3) + time_course = { + "donor": labels("D", donor), + "time": ["t0", "t1", "t2"] * 6, + "age": np.round(rng.uniform(20.0, 80.0, 6), 1)[donor].tolist(), + } + sizes = [1, 2, 2, 3, 3, 4, 4, 5] + unbalanced = { + "donor": labels("D", np.repeat(np.arange(8), sizes)), + "condition": rng.choice(["ctrl", "stim"], sum(sizes)).tolist(), + "x": np.round(rng.normal(0.0, 1.0, sum(sizes)), 3).tolist(), + } + crossed = { + "donor": labels("D", np.repeat(np.arange(8), 4)), + "batch": labels("B", rng.integers(0, 4, 32)), + "condition": ["ctrl", "ctrl", "stim", "stim"] * 8, + "x": np.round(rng.normal(0.0, 1.0, 32), 3).tolist(), + } + nested = { + "donor": labels("D", np.repeat(np.arange(4), 6)), + "sample": labels("S", np.repeat(np.arange(3), 2).tolist() * 4), + "time": np.repeat(["t0", "t1", "t2"], 2).tolist() * 4, + } + # F-tests: two contrasts, and three that span two dimensions. + time_f = { + "time": [[0, 1, 0, 0], [0, 0, 1, 0]], + "pairs": [[0, 1, 0, 0], [0, 0, 1, 0], [0, -1, 1, 0]], + } + two_f = {"both": [[0, 1, 0], [0, 0, 1]]} + specifications = [ # samples, formula, random, F-tests, a contrast, donor SD, weights + ( + "time_course", + time_course, + "~ time + age", + "(1 | donor)", + time_f, + [0, -1, 1, 0], + 2.0, + False, + ), + ( + "unbalanced", + unbalanced, + "~ condition + x", + "(1 | donor)", + two_f, + [0, 1, -1], + 2.0, + False, + ), + ( + "crossed", + crossed, + "~ condition + x", + "(1 | donor) + (1 | batch)", + two_f, + [0, 1, -1], + 2.0, + False, + ), + ( + "nested", + nested, + "~ time", + "(1 | donor/sample)", + {"time": [[0, 1, 0], [0, 0, 1]]}, + [0, -1, 1], + 2.0, + False, + ), + ( + "boundary", + time_course, + "~ time + age", + "(1 | donor)", + time_f, + [0, -1, 1, 0], + 0.0, + False, + ), + ( + "weights", + time_course, + "~ time + age", + "(1 | donor)", + time_f, + [0, -1, 1, 0], + 2.0, + True, + ), + ( + "weights_crossed", + crossed, + "~ condition + x", + "(1 | donor) + (1 | batch)", + two_f, + [0, 1, -1], + 2.0, + True, + ), + ] + cases = {} + for name, samples, formula, random, f, contrast, sd, weighted in specifications: + y = simulate(rng, Design(samples, formula, random=random), n_genes, sd) + w = None + if weighted: # positive weights that rise with expression, like voom's + w = np.exp(rng.normal(0.0, 0.5, y.shape)) / ( + 0.2 + 4.0 / (1.0 + np.exp(y - 5.0)) + ) + y = y + rng.normal(0.0, 0.3, y.shape) / np.sqrt(w) + case = { + "samples": samples, + "formula": formula, + "random": random, + "f": f, + "contrast": contrast, + "y": y, + "w": w, + "variant": "average" if weighted else "none", + } + cases[name] = case + if weighted: + cases[f"{name}_precision"] = case | {"variant": "precision"} + return cases + + +def run_fixed(case: dict[str, Any]) -> tuple[dict[str, Any], dict[str, Any]]: + """CytePool's fits and tests of a case, and what R needs to evaluate lme4 there.""" + design = Design(case["samples"], case["formula"], random=case["random"]) + fit = fit_mixed_model(case["y"], design, weights=case["w"]) + options = { + "df_method": "kenward_roger", + "kenward_roger_weights": "precision" + if case["variant"] == "precision" + else "average", + } + ours = { + "fit": fit, + "kr": fixed_effect_tests(fit, **options), + "contrast": fixed_effect_tests( + fit, np.array(case["contrast"], float), **options + ), + "f": { + key: ( + fixed_effect_f_test(fit, np.array(L, float).T, **options), + fixed_effect_f_test(fit, np.array(L, float).T), + ) + for key, L in case["f"].items() + }, + } + ratios = fit.variance_components / fit.residual_variance[:, None] + if case["variant"] == "average": + # lme4 fits the weights divided by their mean: its residual variance is that of + # an observation of average weight. + ratios = ratios * fit.mean_weights[:, None] + to_r = { + "samples": case["samples"], + "random": case["random"].replace("(1 |", "(0 + sw |"), + "x": design.matrix, + "y": case["y"].T, + "w": None if case["w"] is None else case["w"].T, + "theta": np.sqrt(ratios), + "components": list(fit.component_names), + "contrast": case["contrast"], + "f": case["f"], + "variant": case["variant"], + } + return ours, to_r + + +def compare_fixed( + name: str, ours: dict[str, Any], ref: dict[str, Any], results: Results +) -> None: + kr, contrast, fit = ours["kr"], ours["contrast"], ours["fit"] + adjusted = np.array(ref["adjusted"], float) + coef = np.array( + ref["coef"], float + ) # genes x coefficients x (estimate, se, df, t, p) + one = np.array(ref["contrast"], float) + scale = np.max(np.abs(adjusted), axis=(1, 2))[:, None, None] + # The tests of every coefficient use the adjusted covariance of the coefficients. + results.check( + f"{name}: adjusted covariance of the coefficients", + "pbkrtest vcovAdj", + "max diff. / max abs, per gene", + relative(kr.covariance, adjusted, scale), + EXACT, + ) + results.check( + f"{name}: Kenward-Roger t-tests of the coefficients and a contrast: estimates, standard " + "errors, t", + "lmerTest summary, contest1D", + f"{SCALED} (estimates: of max(abs, se))", + max( + relative(kr.estimate, coef[..., 0], coef[..., 1]), + relative(kr.standard_error, coef[..., 1]), + relative(kr.t, coef[..., 3], 1.0), + relative(contrast.estimate[:, 0], one[:, 0], one[:, 1]), + relative(contrast.standard_error[:, 0], one[:, 1]), + relative(contrast.t[:, 0], one[:, 3], 1.0), + ), + EXACT, + ) + results.check( + f"{name}: Kenward-Roger df and p-values", + "lmerTest summary, contest1D", + f"df: {RELATIVE}; p: {P_VALUE}", + max( + relative(kr.df, coef[..., 2]), + p_value_difference(kr.p_value, coef[..., 4]), + relative(contrast.df[:, 0], one[:, 2]), + p_value_difference(contrast.p_value[:, 0], one[:, 4]), + ), + EXACT, + ) + for key, (krf, satf) in ours["f"].items(): + # KRmodcomp's stat, ndf, ddf, F.scaling, p; contestMD's F, NumDF, DenDF, p by + # Kenward-Roger, then by Satterthwaite. + f = np.array([gene[key] for gene in ref["f"]], float) + numerator = float(np.any(f[:, [1, 6, 10]] != krf.df1)) * np.inf + results.check( + f"{name}, F-test '{key}': Kenward-Roger F, df, scaling, p-value; lmerTest's the same", + "pbkrtest KRmodcomp, lmerTest contestMD", + f"{RELATIVE}; p: {P_VALUE}", + max( + relative(krf.f_statistic, f[:, 0]), + relative(krf.df2, f[:, 2]), + numerator, + relative(krf.scaling, f[:, 3]), + p_value_difference(krf.p_value, f[:, 4]), + relative(f[:, 5], f[:, 0]), + relative(f[:, 7], f[:, 2]), + ), + EXACT, + ) + results.check( + f"{name}, F-test '{key}': Satterthwaite F", + "lmerTest contestMD", + RELATIVE, + relative(satf.f_statistic, f[:, 9]), + EXACT, + ) + results.check( + f"{name}, F-test '{key}': Satterthwaite denominator df and p-value", + "lmerTest contestMD", + f"df: {RELATIVE}; p: {P_VALUE}", + max( + relative(satf.df2, f[:, 11]), p_value_difference(satf.p_value, f[:, 12]) + ), + SATTERTHWAITE, + ) + optimum = np.array(ref["optimum"], float) + results.check( + f"{name}: Kenward-Roger standard errors and df at lme4's optimum", + "lmerTest summary", + RELATIVE, + max( + relative(kr.standard_error, optimum[..., 0]), + relative(kr.df, optimum[..., 1]), + ), + CLOSE, + ) + results.other( + f"{name}: fits with a variance at zero", + "CytePool", + "genes", + f"{int(np.count_nonzero(fit.boundary.any(axis=1)))} of {fit.boundary.shape[0]}", + ) + if kr.kenward_roger_weights == "average" and fit.weights is not None: + # How the two treatments of the weights differ, and from the unadjusted + # standard errors. + precision = fixed_effect_tests( + fit, df_method="kenward_roger", kenward_roger_weights="precision" + ) + unadjusted = fixed_effect_tests(fit).standard_error + se = kr.standard_error / precision.standard_error + df = kr.df / precision.df + results.other( + f"{name}: every observation at the average weight against the precision " + "weights", + "CytePool", + "standard errors, df: ratios, median (range); standard errors below the " + "unadjusted ones: average, precision", + f"{np.median(se):.3f} ({se.min():.3f} to {se.max():.3f}), {np.median(df):.3f} " + f"({df.min():.3f} to {df.max():.3f}); " + f"{np.mean(kr.standard_error < unadjusted):.0%}, " + f"{np.mean(precision.standard_error < unadjusted * (1 - 1e-12)):.0%}", + ) + + +# dream: studies, and cubes for dream_voom. + + +def counts_of( + seed: int, per_donor: list[int], n_genes: int, offset: float = 0.0 +) -> tuple[np.ndarray, dict[str, Any]]: + """Negative binomial counts with a donor effect, three time points and an age.""" + rng = np.random.default_rng(seed) + donor = np.repeat(np.arange(len(per_donor)), per_donor) + n = donor.size + time = np.array(["t0", "t1", "t2"])[np.arange(n) % 3] + age = np.round(rng.normal(size=n), 3) + sd = rng.uniform(0.0, 0.6, n_genes) * (rng.random(n_genes) > 0.2) + log_mean = ( + rng.normal(3.0, 1.5, n_genes) + + (rng.normal(size=(len(per_donor), n_genes)) * sd)[donor] + ) + for level in ("t1", "t2"): + log_mean += np.outer(time == level, rng.normal(0.0, 0.5, n_genes)) + log_mean += np.outer(age, rng.normal(0.0, 0.2, n_genes)) + depth = rng.uniform(5e5, 2e6, n) + counts = rng.negative_binomial( + 10, 10 / (10 + np.exp(log_mean) * (depth / depth.mean())[:, None]) + ) + return counts, { + "donor": labels("D", donor), + "time": time.tolist(), + "age": (age + offset).tolist(), + } + + +FORMULA, FIXED, RANDOM = "~ time + age + (1 | donor)", "~ time + age", "(1 | donor)" +CONTRASTS = { + "t2_vs_t1": [0, -1, 1, 0], + "(Intercept)": [1, 0, 0, 0], + "timet1": [0, 1, 0, 0], + "timet2": [0, 0, 1, 0], + "age": [0, 0, 0, 1], +} +# Studies: samples per donor, and the values of dream's ddf to run. +STUDIES = { + "n12": ([2] * 6, ["Kenward-Roger", "adaptive", "Satterthwaite"]), + "n18": ([3] * 6, ["Kenward-Roger", "adaptive", "Satterthwaite"]), + "n19": ([4, 4, 4, 4, 3], ["adaptive"]), + "n20": ([4] * 5, ["adaptive"]), + "n24": ([3] * 8, ["Kenward-Roger"]), +} + + +def study_cases() -> dict[str, dict[str, Any]]: + cases = {} + for seed, (name, (per_donor, ddf)) in enumerate(STUDIES.items(), start=20261101): + counts, samples = counts_of(seed, per_donor, 40) + cases[name] = { + "counts": counts, + "samples": samples, + "formula": FORMULA, + "ddf": ddf, + "contrasts": {"t2_vs_t1": "timet2 - timet1"}, + "f_tests": {"time": ["timet1", "timet2"]}, + "lib_size": library_sizes(counts.astype(float)), + } + return cases + + +def ours_study(case: dict[str, Any], method: str) -> tuple[Any, Any, Any]: + """CytePool's dream of a study: voom's weights, the moderated tests of every + coefficient and contrast, and the moderated F-test of time.""" + design = Design(case["samples"], FIXED, random=RANDOM) + voom = voom_with_dream_weights(case["counts"], design, span=0.5) + fit = fit_mixed_model(voom.log_cpm, design, weights=voom.weights) + tests = fixed_effect_tests(fit, CONTRASTS, df_method=method) + time = fixed_effect_tests( + fit, {k: CONTRASTS[k] for k in ("timet1", "timet2")}, df_method=method + ) + return voom, moderated_fixed_effect_tests(tests), moderated_fixed_effect_tests(time) + + +def compare_f( + label: str, + f_value: Any, + p_value: Any, + dream_f: dict[str, Any], + first_df: float, + results: Results, +) -> None: + """The moderated F-test against topTable's: the statistic, and the p-values that + dream computes on one gene's degrees of freedom.""" + results.check( + f"{label}: moderated F-statistic", + "topTable F", + RELATIVE, + relative(f_value, dream_f["F"]), + CLOSE, + ) + dream_p = np.asarray(f_sf(np.asarray(dream_f["F"], float), 2.0, first_df)) + ratio = np.asarray(p_value) / np.asarray(dream_f["p"], float) + results.other( + f"{label}: p-values", + "topTable P.Value", + "dream's on the first gene's df.total, for every gene: rel. difference; CytePool's over dream's", + f"{relative(dream_p, dream_f['p']):.0e}; {np.nanmin(ratio):.3g} to {np.nanmax(ratio):.3g}", + relative(dream_p, dream_f["p"]) < EXACT, + ) + + +def compare_study( + name: str, case: dict[str, Any], ref: dict[str, Any], results: Results +) -> list[str]: + n = len(case["samples"]["donor"]) + voom, *_ = ours_study(case, "satterthwaite") + results.check( + f"{name}: precision weights", + "voomWithDreamWeights", + RELATIVE, + relative(voom.weights, np.asarray(ref["weights"]).T), + CLOSE, + ) + found = {} + for ddf in case["ddf"]: + r = ref[ddf] + method = ( + "kenward_roger" + if ddf == "Kenward-Roger" or (ddf == "adaptive" and n < 20) + else "satterthwaite" + ) + _, moderated, time = ours_study(case, method) + label = f"{name}, ddf = {ddf}" + ( + f" (CytePool: {method})" if ddf == "adaptive" else "" + ) + columns = [r["names"].index(key) for key in CONTRASTS] + pick = { + key: np.asarray(r[key], float)[:, columns] + for key in ("se", "df", "t", "lods", "df_total", "p", "z") + } + results.check( + f"{label}: standard errors and df of the tests", + "dream", + RELATIVE, + max( + relative(moderated.tests.standard_error, pick["se"]), + relative(moderated.tests.df, pick["df"]), + ), + CLOSE, + ) + results.check( + f"{label}: moderated t, B, z.std; df.total, p-values", + "eBayes, topTable", + f"t, B, z: {SCALED}; df.total: {RELATIVE}; p: {P_VALUE}", + max( + relative(moderated.t, pick["t"], 1.0), + relative(moderated.lods, pick["lods"], 1.0), + relative(moderated.z, pick["z"], 1.0), + relative(moderated.df_total, pick["df_total"]), + p_value_difference(moderated.p_value, pick["p"]), + ), + CLOSE, + ) + first = float(np.asarray(r["df_total"], float)[0, r["names"].index("timet1")]) + compare_f( + f"{label}, F-test of time", + time.f_statistic, + time.f_p_value, + r["f"]["time"], + first, + results, + ) + found[ddf] = moderated + if "Satterthwaite" not in found: + return [] + kr, satterthwaite = found["Kenward-Roger"], found["Satterthwaite"] + ratio = kr.p_value[:, 0] / satterthwaite.p_value[:, 0] + se = kr.tests.standard_error[:, 0] / satterthwaite.tests.standard_error[:, 0] + return [ + f"study {name}, contrast t2_vs_t1", + f"{np.nanmedian(ratio):.3g}; " + f"{np.nanquantile(ratio, 0.05):.3g} to {np.nanquantile(ratio, 0.95):.3g}; " + f"{np.nanmedian(se):.3g}, {np.mean(se < 1.0):.0%}", + ] + + +def cube_cases() -> list[dict[str, Any]]: + """Cubes for dream_voom: cell types with 12 and 24 samples left, and one of 15 + samples whose age, 32 years on average and 1 apart, gives condition numbers near + 1e6. And, for R, each cell type's samples, features and library sizes.""" + contrasts = { + "t1_vs_t0": "timet1", + "t2_vs_t0": "timet2", + "t2_vs_t1": "timet2 - timet1", + } + cases = [] + for name, per_donor, offset, seeds in ( + ("time_course", [3] * 8, 0.0, (20261111, 20261112)), + ("ill_conditioned", [3] * 5, 32.0, (20261113,)), + ): + layers = [counts_of(seed, per_donor, 60, offset)[0] for seed in seeds] + samples = counts_of(seeds[0], per_donor, 60, offset)[1] + n = len(samples["donor"]) + cells = np.random.default_rng(20261114).integers(20, 300, size=(n, len(seeds))) + if len(seeds) == 2: + cells[n // 2 :, 0] = 5 # "small": half the donors have too few cells + samples["sample_id"] = [f"S{i:02d}" for i in range(n)] + groups = ["small", "large"] if len(seeds) == 2 else ["cells"] + cube = Cube( + {"counts": np.stack(layers)}, + samples=samples, + groups={"group_id": groups}, + features={"feature_id": [f"gene{i}" for i in range(60)]}, + n_cells=cells, + ) + design = Design( + cube.samples, + FIXED, + random=RANDOM, + levels={"time": ["t0", "t1", "t2"]}, + contrasts=contrasts, + ) + inputs = {} + for g, group in enumerate(groups): + values = np.asarray(cube.layer("counts", g), float) + used = cells[:, g] >= 10 + keep = filter_by_expression(values[used], design=design.matrix[used]) + y = values[used][:, keep] + inputs[group] = { + "used": used, + "keep": keep, + "counts": y.T, + "lib_size": library_sizes(y) * tmm_factors(y), + "samples": { + k: np.asarray(samples[k])[used].tolist() + for k in ("donor", "time", "age") + }, + "formula": FORMULA, + "ddf": ["adaptive"], + "condition": offset > 0, + "contrasts": contrasts, + "f_tests": {"time": ["t1_vs_t0", "t2_vs_t0"]}, + } + cases.append({"name": name, "cube": cube, "design": design, "inputs": inputs}) + return cases + + +def compare_cube( + case: dict[str, Any], refs: dict[str, Any], results: Results +) -> list[list[str]]: + """dream_voom's default against dream on every cell type, and the change from + Satterthwaite's method that the default makes.""" + changes = [] + # span = 0.5, voomWithDreamWeights' span. + result = dream_voom( + case["cube"], + case["design"], + span=0.5, + f_tests={"time": ["t1_vs_t0", "t2_vs_t0"]}, + ) + before = dream_voom( + case["cube"], case["design"], span=0.5, df_method="satterthwaite" + ) + for entry in result.provenance["groups"]: + group = entry["group_id"] + given, ref = case["inputs"][group], refs[f"{case['name']}_{group}"] + r = ref["adaptive"] + label = f"cube {case['name']}, {group} ({len(entry['samples'])} samples: {entry['df_method']})" + tested = np.asarray(result.table(group, "t1_vs_t0")["feature_id"]) + expected = np.array([f"gene{i}" for i in np.flatnonzero(given["keep"])]) + samples = np.asarray(case["design"].samples["sample_id"])[given["used"]] + differences = np.count_nonzero(np.asarray(entry["samples"]) != samples) + differences += ( + np.count_nonzero(tested != expected) + if tested.size == expected.size + else expected.size + ) + results.check( + f"{label}: samples used, features tested", + "the rules of dream_voom", + "differences", + float(differences), + 0.0, + ) + alike = np.ones(tested.size, bool) + if given["condition"]: + # dream decides from R's estimate of the condition number, CytePool from + # the exact one: a gene near the bound can be decided differently. + k = r["names"].index("t1_vs_t0") + df = np.asarray(result.table(group, "t1_vs_t0")["df"]) + alike = np.abs(df - np.asarray(r["df_total"], float)[:, k]) <= 1e-4 * df + kappa = np.asarray(ref["kappa"], float) + near = np.abs(np.log(kappa / 1e6)) < np.log(1.02) + results.other( + f"{label}: genes whose Kenward-Roger tests fall back to Satterthwaite's", + "dream, R's kappa of vcovAdj", + "dream; CytePool; decided differently; all of those within 2% of 1e6", + f"{int(np.count_nonzero(kappa > 1e6))}; {entry['features_satterthwaite_fallback']}; " + f"{int(np.count_nonzero(~alike))}; {bool(np.all(near[~alike]))}", + bool(np.all(near[~alike])) + and entry["features_satterthwaite_fallback"] > 0, + ) + worst = 0.0 + for name in case["design"].contrasts.names: + table, k = result.table(group, name), r["names"].index(name) + for column, key, floor in ( + ("t", "t", 1.0), + ("b_statistic", "lods", 1.0), + ("z_statistic", "z", 1.0), + ("df", "df_total", 0.0), + ): + worst = max( + worst, + relative( + np.asarray(table[column])[alike], + np.asarray(r[key], float)[alike, k], + floor, + ), + ) + worst = max( + worst, + p_value_difference( + np.asarray(table["p_value"])[alike], + np.asarray(r["p"], float)[alike, k], + ), + ) + alike_label = " (genes decided alike)" if given["condition"] else "" + results.check( + f"{label}: t, B, z; df; p-values of every contrast{alike_label}", + "dream, eBayes, topTable", + f"t, B, z: {SCALED}; df: {RELATIVE}; p: {P_VALUE}", + worst, + CLOSE, + ) + if not given["condition"]: + f_table = result.table(group, "time") + first = float( + np.asarray(r["df_total"], float)[0, r["names"].index("t1_vs_t0")] + ) + compare_f( + f"{label}, F-test of time", + f_table["f_statistic"], + f_table["p_value"], + r["f"]["time"], + first, + results, + ) + old, new = before.table(group, "t1_vs_t0"), result.table(group, "t1_vs_t0") + ratio = np.asarray(new["p_value"]) / np.asarray(old["p_value"]) + significant = [ + int(np.nansum(np.asarray(t["adjusted_p_value"]) < 0.05)) for t in (old, new) + ] + changes.append( + [ + f"cube {case['name']}, {group}, contrast t1_vs_t0", + f"{np.nanmedian(ratio):.3g}; {np.nanquantile(ratio, 0.05):.3g} to " + f"{np.nanquantile(ratio, 0.95):.3g}; {significant[0]} and {significant[1]} genes", + ] + ) + return changes + + +def run_r(payload: dict[str, Any]) -> dict[str, Any]: + def plain(value: Any) -> Any: + if isinstance(value, np.ndarray): + return value.tolist() + if isinstance(value, dict): + return { + key: plain(item) + for key, item in value.items() + if key not in ("used", "keep") + } + return value + + with tempfile.TemporaryDirectory() as directory: + root = Path(directory) + (root / "tests.R").write_text(R_SCRIPT) + (root / "input.json").write_text(json.dumps(plain(payload))) + command = [ + os.environ.get("RSCRIPT", "Rscript"), + "--vanilla", + str(root / "tests.R"), + str(root / "input.json"), + str(root / "out.json"), + ] + subprocess.run(command, check=True) + reference: dict[str, Any] = json.loads((root / "out.json").read_text()) + return reference + + +def write_fixtures( + cases: dict[str, Any], refs: dict[str, Any], studies: dict[str, Any], cube: Cube +) -> None: + FIXTURES.mkdir(parents=True, exist_ok=True) + genes = slice(0, N_FIXTURE) + for name, case in cases.items(): + ref = refs["fixed"][name] + arrays: dict[str, Any] = { + "y": case["y"][:, genes], + "adjusted": np.array(ref["adjusted"], float)[genes], + "coef": np.array(ref["coef"], float)[genes], + "contrast": np.array(ref["contrast"], float)[genes], + "contrast_weights": np.array(case["contrast"], float), + } + if case["w"] is not None: + arrays["weights"] = case["w"][:, genes] + for key, L in case["f"].items(): + arrays[f"f {key}"] = np.array([gene[key] for gene in ref["f"]], float)[ + genes + ] + arrays[f"f_matrix {key}"] = np.array(L, float) + np.savez_compressed(FIXTURES / f"{name}.npz", **arrays) + meta = {key: case[key] for key in ("samples", "formula", "random", "variant")} + (FIXTURES / f"{name}.json").write_text(json.dumps(meta) + "\n") + for name in ("n12", "n20"): + case, ref = studies[name], refs["dream"][name] + arrays = {"counts": case["counts"]} + for ddf in case["ddf"]: + r = ref[ddf] + columns = [r["names"].index(key) for key in CONTRASTS] + for key in ("se", "df", "t", "lods", "df_total", "p", "z"): + arrays[f"{ddf} {key}"] = np.asarray(r[key], float)[:, columns] + arrays[f"{ddf} F"] = np.asarray(r["f"]["time"]["F"], float) + np.savez_compressed(FIXTURES / f"study_{name}.npz", **arrays) + (FIXTURES / f"study_{name}.json").write_text(json.dumps(case["samples"]) + "\n") + arrays = {"counts": np.asarray(cube.layer("counts")), "n_cells": cube.n_cells} + for group in ("small", "large"): + r = refs["dream"][f"time_course_{group}"]["adaptive"] + for name in ("t1_vs_t0", "t2_vs_t1"): + k = r["names"].index(name) + for key in ("t", "df_total", "p", "lods", "z"): + arrays[f"{group} {name} {key}"] = np.asarray(r[key], float)[:, k] + arrays[f"{group} time F"] = np.asarray(r["f"]["time"]["F"], float) + np.savez_compressed(FIXTURES / "cube.npz", **arrays) + (FIXTURES / "cube.json").write_text( + json.dumps({k: cube.samples[k].tolist() for k in cube.samples.columns}) + "\n" + ) + print(f"wrote fixtures to {FIXTURES}") + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + parser.add_argument( + "--write-fixtures", action="store_true", help="write test fixtures" + ) + arguments = parser.parse_args() + cases = fixed_cases(20261031, N_GENES) + fixed = {name: run_fixed(case) for name, case in cases.items()} + studies = study_cases() + cubes = cube_cases() + dream_inputs = { + name: case | {"counts": case["counts"].T} for name, case in studies.items() + } + for cube_case in cubes: + dream_inputs |= { + f"{cube_case['name']}_{g}": given + for g, given in cube_case["inputs"].items() + } + refs = run_r( + { + "fixed": {name: to_r for name, (_, to_r) in fixed.items()}, + "dream": dream_inputs, + } + ) + print( + "Reference: " + + ", ".join(f"{name} {version}" for name, version in refs["versions"].items()) + ) + results = Results() + for name, (ours, _) in fixed.items(): + compare_fixed(name, ours, refs["fixed"][name], results) + changes = [] + for name, case in studies.items(): + changes += [compare_study(name, case, refs["dream"][name], results)] + for cube_case in cubes: + changes += compare_cube(cube_case, refs["dream"], results) + for change in changes: + if change: + results.other( + f"{change[0]}: the default against Satterthwaite's method", + "CytePool", + "p-values over Satterthwaite's: median; 5% to 95%; genes with adjusted p below 0.05 " + "before and now" + if change[0].startswith("cube") + else "p-values over Satterthwaite's: median; 5% to 95%; standard errors over " + "Satterthwaite's: median, share below 1", + change[1], + ) + results.print() + print(f"\n{len(results.rows)} checks, {results.failures} failed") + if arguments.write_fixtures: + if results.failures: + print("not writing fixtures: some checks failed") + else: + write_fixtures(cases, refs, studies, cubes[0]["cube"]) + return 1 if results.failures else 0 + + +if __name__ == "__main__": + raise SystemExit(main())