From b582f1b70ad2eb2d2c1de6a411d6a2b7be3c662f Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 25 Sep 2026 20:31:09 +0000 Subject: [PATCH] fix(analysis): stop returning placeholder statistics from run_analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit run_analysis derived its p-values from 0.01 + (hash(taxon_id) % 100)/1000.0, padj from pvalue * 1.5 and log2FoldChange from (hash % 20)/10 - 1: every one of those is a function of the feature's NAME, not of the counts, and they were returned to callers as results. The placeholder is deleted rather than deprecated. src/analysis/estimation.jl runs the declared model for real (MASS::glm.nb with a required offset for nb_glm, stats::lm for CLR/ILR, stats::glm binomial for logistic on a 0/1 response), records R and MASS versions, offset hash and result hashes, excludes failed fits from the BH family while counting them, and returns an explicit not_run with a reason — never an empty table that reads as 'nothing was significant' — when no design or no runtime is available. Conditions published in docs/statistics/method-conditions/parametric-fits.md (catalogue item 2). Evidence in test/unit/test_estimation.jl: written-in known effects, a direct R reference compared coefficient by coefficient, BH against R's p.adjust, negative controls for every refusal, and a source-level guard that the placeholder does not come back. Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- CHANGELOG.md | 29 + docs/statistics/method-catalogue-v1.md | 6 + .../method-conditions/parametric-fits.md | 186 ++++ src/MetaManifold.jl | 2 + src/analysis/Execution.jl | 152 ++-- src/analysis/estimation.jl | 849 ++++++++++++++++++ test/runtests.jl | 1 + test/unit/test_estimation.jl | 435 +++++++++ test/unit/test_execution.jl | 15 +- 9 files changed, 1611 insertions(+), 64 deletions(-) create mode 100644 docs/statistics/method-conditions/parametric-fits.md create mode 100644 src/analysis/estimation.jl create mode 100644 test/unit/test_estimation.jl diff --git a/CHANGELOG.md b/CHANGELOG.md index 14ae9c74..8b9b6971 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,35 @@ types, tests, infrastructure, and alignment. ## [Unreleased] +### Fixed — the analysis path no longer returns placeholder statistics (2026-09-25) + +- **`run_analysis` computed nothing and returned numbers anyway.** It derived + `pvalue` from `0.01 + (hash(taxon_id) % 100) / 1000.0`, `padj` from + `pvalue * 1.5` and `log2FoldChange` from `(hash % 20) / 10 - 1`: every one of + those is a deterministic function of the feature's *name*, carrying no + information about the counts, and they were returned to callers as results. The + placeholder is deleted, not deprecated, and `test/unit/test_estimation.jl` + reads `Execution.jl` as source and asserts it does not come back. +- **Real estimation** (`src/analysis/estimation.jl`, catalogue item 2): per-feature + `MASS::glm.nb` for `nb_glm` (with a required offset), `stats::lm` for + `clr_lm`/`ilr_lm`, and `stats::glm(family = binomial)` for `logistic` on a 0/1 + response — all through the shared R runtime lock, all with R and MASS versions, + offset hash and result-table hashes in the provenance. +- **Unsuccessful states are states.** A feature whose fit fails gets + `status = "failed"`, a note naming the cause, and `null` for every statistic; it + is excluded from the BH family and counted. A run that cannot happen at all — + no design, R unreachable, R busy past the timeout — returns `status = + "not_run"` with a reason and an empty result set, never an empty table that + reads as "nothing was significant". +- **Refusals instead of substitutions**: unsupported formula syntax + (interactions, transformations, random effects, nesting), a `glmGamPoi`/ + `local`/`mean`/`pooled` dispersion method that has no implementation here, an + offset on a non-count model, a count model without one, and a binomial fit on + proportions all raise with the reason attached. +- BH is implemented in Julia and compared against R's `p.adjust(method = "BH")` + in the tests; conditions are published in + `docs/statistics/method-conditions/parametric-fits.md`. + ### Added — Type-system engineering series (2026-09) - **Epistemic claims with receipts**: added `Standpoint`, `TaxonWarrant`, diff --git a/docs/statistics/method-catalogue-v1.md b/docs/statistics/method-catalogue-v1.md index c9d002a5..352195e0 100644 --- a/docs/statistics/method-catalogue-v1.md +++ b/docs/statistics/method-catalogue-v1.md @@ -43,6 +43,12 @@ means something specific. maximum likelihood where the model justifies it, with identifiability, boundary-estimate and convergence checks that produce *unsuccessful states* rather than plausible parameters. + **Implemented 2026-09-25** in `src/analysis/estimation.jl`, reached through + `Execution.run_analysis`; conditions published in + [`method-conditions/parametric-fits.md`](method-conditions/parametric-fits.md). + It replaces a placeholder that derived p-values from `hash(taxon_id)`. What it does + *not* yet have is issue #1's independent statistical review, which is why the + catalogue's other items remain gated. 3. **Nonparametric tests** — permutation/bootstrap, with resampling units and exchangeability stated, seeds recorded, and Monte Carlo limits reported alongside the estimate. diff --git a/docs/statistics/method-conditions/parametric-fits.md b/docs/statistics/method-conditions/parametric-fits.md new file mode 100644 index 00000000..7db29c01 --- /dev/null +++ b/docs/statistics/method-conditions/parametric-fits.md @@ -0,0 +1,186 @@ + +# Method conditions — parametric fits + +**Catalogue item 2** of [`method-catalogue-v1.md`](../method-catalogue-v1.md), approved +2026-09-22. **Published before implementation**, as the catalogue requires: the +implementation in `src/analysis/estimation.jl` is held to this document, so it is written to +be checked against rather than admired. + +**Status: implemented** in `src/analysis/estimation.jl`, called by +`Execution.run_analysis`. Evidence, as delivered, is in `test/unit/test_estimation.jl`: +known answers written into the data rather than read back out of the fit; an independent +reference (the same data fitted directly in R, compared coefficient by coefficient, and BH +compared against R's `p.adjust`); negative controls for every refusal path; and a +source-level guard that the placeholder statistics this replaced do not come back. + +**What this document does not claim.** There has been no independent statistical review of +this layer. Issue #1 requires one before its acceptance criteria are met, and that is still +outstanding. Nothing below is a claim that the models are *appropriate* for a given dataset; +it is a statement of what is computed, what is refused, and what is argued for the user. + +## What ran before, and why this document exists + +`Execution.run_analysis` used to return, per feature, a p-value computed as +`0.01 + (hash(taxon_id) % 100) / 1000.0`, an adjusted p-value of `p * 1.5`, a log2 fold +change of `(hash % 20) / 10 - 1` and a `baseMean` of `100 + hash % 1000`, all derived from +the *name* of the feature. Those numbers are deterministic, look plausible in a table, and +carry no information about the counts. They were labelled a stub in the source and handed to +the caller as results. + +A stub that is returned as a result is a wrong answer. This layer replaces it: where a fit +cannot be run or fails, the output says so and carries no number at all. + +## What this method is + +Per-feature estimation of a declared model, one row per feature, with the design named by the +configuration's formula. Four methods, and for each one the response it accepts: + +| Method | Response | Model | Fit in | +| --- | --- | --- | --- | +| `nb_glm` | raw counts (non-negative) | negative binomial GLM, log link | R `MASS::glm.nb` | +| `clr_lm` | CLR-transformed values | Gaussian linear model | R `stats::lm` | +| `ilr_lm` | ILR balances | Gaussian linear model | R `stats::lm` | +| `logistic` | 0/1 presence | binomial GLM, logit link | R `stats::glm(family = binomial)` | + +Every fit is deterministic: there is no resampling, no permutation and no random start. The +provenance says so, and records the seed the caller passed as *not a parameter of any number +here* rather than implying it was used. + +### Response types accepted, and what happens to zeros + +- `nb_glm` requires counts. Zeros are handled before the fit, by + `prepare_analysis_table`: either left alone (`pseudocount` is not applied to the raw counts + in this path), or replaced by the declared zero policy, whose parameters are recorded in + the config. A response that is constant, or has fewer than three finite values, is not + fitted at all — see *Unsuccessful states* below. +- `clr_lm` and `ilr_lm` receive a table that has already been transformed. A zero reaching + the transform is what the zero policy is for; `log(0)` is not silently conventionalised + anywhere in this layer. +- `logistic` requires a 0/1 response and **refuses anything else**. A binomial fit on + proportions needs the number of trials per sample, and this layer will not invent weights: + the prepared table must come from the `presence_absence` transform. + +### Study designs supported + +Additive main effects over metadata columns: `~ a`, `~ a + b`, `~ a + b + c`, where every +term is a plain column name. The first term is the **primary contrast** — the coefficient +reported per feature and the coefficient the BH family is built from. For a factor, R's +default treatment coding applies: the coefficient is the contrast between the first two +levels, and the levels are recorded in the provenance so the sign is interpretable. + +Refused, by name, with a reason (`UnsupportedFormula`): + +- interactions (`*`, `:`), nesting (`/`), polynomial terms (`^`), transformations (`I()`, + `log()`, `poly()`), random effects (`(1 | batch)`), functions of several arguments, and + `-` in any position; +- a formula with no `~`, an empty right-hand side, a duplicated term, or a dangling `+`; +- `$`, backticks, `;`, backslashes and newlines anywhere in the formula. + +**Pairing, blocking and repeated measures are not supported.** A paired design cannot be +expressed in this grammar; entering the subject as a covariate fits a fixed effect per +subject, which is not the same analysis and would be reported as though it were. Until a +paired model is implemented and validated, a paired study gets a refusal, not a p-value. + +### Overdispersion and depth + +- `nb_glm` estimates one dispersion parameter per feature by maximum likelihood + (`MASS::glm.nb`), that is `dispersion_method = "parametric"`. The other names accepted by + `AnalysisConfig` — `local`, `mean`, `pooled`, `glmGamPoi` — are **refused by name** at the + point of fitting rather than aliased to it (issue #21). Dispersion sets every standard + error in the fit; a run that asked for one estimator and received another would report + different intervals and different multiple-testing outcomes under a label saying otherwise. +- **An offset is required for `nb_glm` and refused for every other method.** The offset is the + log-scale factor computed in `prepare_analysis_table`: `log(library size)` for + `none`/`TSS`/`CSS`/`RSS`, `log(size factor)` for `size_factors`. Raw counts fitted without + one would attribute sequencing depth to biology; a Gaussian fit on CLR units has nothing to + offset, and accepting the parameter would imply it had been used. +- Dispersion at a boundary is reported: theta ≥ 1e7 gives `status = "boundary"` with a note + that the negative binomial model has collapsed onto Poisson; theta ≤ 1e-8 gives + `status = "boundary"` with a note that the standard error is not trustworthy. +- `logistic` separation (R's "fitted probabilities numerically 0 or 1") gives + `status = "boundary"` and a note saying the Wald test is not a valid p-value. It is never + presented as a clean result. + +### Uncertainty, effect sizes and multiple testing + +- **Effect size, stated in units.** `estimate_scale` names what the number is in: + `log_counts` (so `log2FoldChange = estimate / log(2)`, computed only for `nb_glm`), + `clr_difference`, `ilr_balance`, `log_odds` (with `odds_ratio` for `logistic`). A CLR + difference is not a fold change and is not labelled as one. +- **Intervals are Wald**, 95%: `estimate ± 1.96·SE` for the z-based fits, the t quantile for + the Gaussian fits. Profile-likelihood intervals are not computed — they cost orders of + magnitude more per feature, and the limitation is stated here rather than left implicit. +- **BH is mandatory** (`config.correction`). Adjusted p-values are computed by + `Estimation.bh_adjust`, a step-up implementation with the monotonicity constraint enforced, + tested against R's `p.adjust(method = "BH")`. Disabling BH requires the acknowledgement + token and produces a DANGER banner; when it is disabled, `padj` is `null` rather than equal + to `pvalue`, because absent and uncorrected are different statements. + +### Unsuccessful states, and the shape of a failure + +Three outcomes, all recorded in the result and in `diagnostics.checks["estimation"]`: + +1. The fit ran. Each feature row carries `status` ∈ {`ok`, `boundary`} and its statistics. +2. Some features failed. Those rows carry `status = "failed"`, a `note` naming the cause, and + `null` for every statistic — never a zero, never a NaN, never a substituted estimate. They + are excluded from the BH family, the exclusion count is reported as `n_failed`, and the + family is described as "features with a testable fit" so that a large exclusion rate is + visible rather than flattering. +3. Nothing ran. `status = "not_run"`, `results` empty, and a reason: no metadata (no design), + a formula term that is not in the metadata, a primary term with fewer than two distinct + values, the R runtime busy past the timeout (`RBusyError`), R unreachable, or any other + failure to execute. **An empty result is a statement about the run**, and it is never + presented as "no features were significant". + +### Provenance written for every run + +Method, engine, formula, terms, primary term and its contrast and levels, `estimate_scale`, +dispersion method, offset kind and the SHA-256 of the offset vector, correction and alpha, +R version, MASS version, the SHA-256 of the per-feature table and of the coefficient table as +they crossed the boundary, the statement that the fits are deterministic, and the interval +method. + +### Computational limits + +Fits run inside the shared R runtime lock (`RRuntime`), which serialises them against the +pipeline and every other analysis. The wait is bounded by `Execution.ESTIMATION_R_WAIT_SECONDS` +(default 10 s), and exceeding it is an explicit `not_run`, not a partial result. Per-feature +cost is one R model fit; no per-feature memory bound has been measured, and none is claimed +here. A feature count in the tens of thousands will be slow, and slow is not the same as +failed: the run either completes or reports why it did not. + +## Evidence, as delivered + +In `test/unit/test_estimation.jl`: + +- **Known answers** — a table with a written-in four-fold increase, a written-in no-effect + feature and a written-in ten-fold decrease; the recovered log2 fold changes are checked + against the arithmetic of the means, not against a recording of a previous run. +- **Independent reference** — the same counts fitted directly in R (`MASS::glm.nb` with the + same offset) and compared coefficient by coefficient, standard error by standard error and + p-value by p-value; and BH compared against R's `p.adjust`. +- **Negative controls** — every refusal above is asserted to throw, and every refusal's + message is asserted to name what it refused. A constant feature is asserted to return + `status = "failed"` with no numbers. A missing design, a wrong-length column and a + single-level grouping are asserted to return `not_run`, not a fit. +- **Regression guard on the placeholder** — `Execution.jl` is read as source and asserted not + to contain `hash(taxon_id)`, and `estimation.jl` is asserted not to contain `hash(` or + `rand(`. + +## Explicitly not claimed + +- No Poisson, negative binomial or binomial model is *endorsed* by this document for + sequencing data; the catalogue names them candidates requiring review, and this + implementation does not change that. +- Wald intervals and z-based p-values are asymptotic. Small samples per group are exactly + where they are weakest, and `min_samples_per_group` below 3 already raises a DANGER flag in + the configuration layer. +- The BH family here is one test per feature on the primary contrast. Testing several + contrasts per feature needs the contrast derivation and family accounting that issue #4 + covers; until then, other coefficients are reported but not adjusted as a family of their + own. +- No equivalence testing, no power analysis, no multiplicity across several analyses of the + same data. +- ELIGIBLE-FOR-REVIEW: independent statistical review, as issue #1 requires, has not happened. diff --git a/src/MetaManifold.jl b/src/MetaManifold.jl index c724aba7..4f609440 100644 --- a/src/MetaManifold.jl +++ b/src/MetaManifold.jl @@ -35,6 +35,8 @@ include("analysis/diversity.jl") include("analysis/analysis.jl") include("analysis/AnalysisConfig.jl") include("analysis/clade_cumulus.jl") +# Estimation fits what AnalysisConfig declares and Execution runs it, so it sits between them. +include("analysis/estimation.jl") include("analysis/Execution.jl") end diff --git a/src/analysis/Execution.jl b/src/analysis/Execution.jl index 058c1c9c..0cf6336d 100644 --- a/src/analysis/Execution.jl +++ b/src/analysis/Execution.jl @@ -1,13 +1,13 @@ # SPDX-License-Identifier: AGPL-3.0-only # SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) """ - Execution — execution harness with abstract AnalysisAdapter and concrete stubs (RAdapter, JuliaAdapter) + Execution — execution harness with abstract AnalysisAdapter and concrete adapters (RAdapter, JuliaAdapter) Implements Milestone 4: -- abstract AnalysisAdapter and concrete stubs RAdapter, JuliaAdapter +- abstract AnalysisAdapter and concrete adapters RAdapter, JuliaAdapter - prepare_analysis_table: applies declared transform, offset, zero_policy, epistemic filtering with Advanced section validation -- run_analysis: stub that records full manifest and diagnostics, hard-stop with DANGER banner on failures +- run_analysis: runs the declared model through `Estimation` and records what it did; hard-stops with a DANGER banner on failures, and reports an unsuccessful state rather than a number when a fit cannot be run - self-diagnostics and safe self-healing - Tests for prepare_analysis_table (CLR with pseudocount=1, epsilon=1e-6, drop vs impute) @@ -28,7 +28,7 @@ Design: - Epistemic filtering: avec_fibre true, epistemic_status present_in_every_admissible_world, min_prevalence, min_abundance, max_features - Self-diagnostics: NaN/Inf, zero variance, all-zero samples/taxa, library size outliers, batch confounding, prevalence/abundance, etc. - Safe self-healing: heal NaN/Inf with epsilon, drop all-zero samples/taxa with warning, record healing in diagnostics, never silent -- run_analysis: stub records full manifest and diagnostics, hard-stop with DANGER banner on failures +- run_analysis: runs the declared model (src/analysis/estimation.jl) and records full manifest and diagnostics; a fit that cannot be run is an unsuccessful state with a reason, never a placeholder number """ module Execution @@ -43,6 +43,7 @@ using Statistics # Use AnalysisConfig from parent module import ..AnalysisConfig using ..Epistemic +using ..Estimation using ..Provenance: probe_metamanifold, probe_host export AnalysisAdapter, RAdapter, JuliaAdapter, @@ -87,7 +88,7 @@ const IMPUTE_POLICY_STRINGS = Dict{String,ImputePolicy}( ) # -------------------------------------------------------------------------- -# Abstract adapter and concrete stubs +# Abstract adapter and concrete adapters # -------------------------------------------------------------------------- """ @@ -99,12 +100,12 @@ Concrete implementations: Every adapter must implement: - prepare_analysis_table (via Execution module, not adapter-specific, but adapter may provide custom transform) -- run_analysis (adapter-specific stub that records manifest and diagnostics) +- run_analysis (adapter-specific; runs the model and records manifest and diagnostics) """ abstract type AnalysisAdapter end """ - RAdapter — stub for R-based execution (DESeq2, edgeR, vegan, metagenomeSeq, etc.) + RAdapter — R-based execution (DESeq2, edgeR, vegan, metagenomeSeq, etc.) Fields: - r_binary: path to R binary (from config/tools.yml) @@ -140,7 +141,7 @@ struct RAdapter <: AnalysisAdapter end """ - JuliaAdapter — stub for pure Julia execution (GLM, MultivariateStats, etc.) + JuliaAdapter — pure Julia execution (GLM, MultivariateStats, etc.) Fields: - method: analysis method (must match AnalysisConfig) @@ -314,7 +315,7 @@ struct ExecutionManifest end """ - ExecutionResult — stub result from run_analysis, records full manifest and diagnostics + ExecutionResult — result of run_analysis: the fit's output, full manifest and diagnostics Fields: - id: UUID4 @@ -322,7 +323,7 @@ Fields: - config_id: UUID4 - config_hash: SHA256 - method: AnalysisMethod -- results: OrderedDict{String,Any} — taxon -> stats (mock for stub) +- results: OrderedDict{String,Any} — feature -> statistics from the fitted model. Empty when nothing could be fitted, never filled with placeholders. - diagnostics: ExecutionDiagnostics - manifest: ExecutionManifest - provenance: OrderedDict @@ -1369,7 +1370,15 @@ function prepare_analysis_table( end # -------------------------------------------------------------------------- -# run_analysis — stub that records full manifest and diagnostics, hard-stop with DANGER banner on failures +# run_analysis — runs the declared model, records full manifest and diagnostics, hard-stop with +# DANGER banner on failures, and reports an unsuccessful state (never a placeholder number) when +# a fit cannot be run. +# +# `Estimation.estimate_models` is the estimator. The block that used to stand here derived its +# p-values from `hash(taxon_id)`, which is a function of the feature name and not of the data: +# it has been removed rather than deprecated, and test/unit/test_estimation.jl asserts at the +# source level that it does not come back. +const ESTIMATION_R_WAIT_SECONDS = Ref(10.0) # -------------------------------------------------------------------------- function is_dangerous_execution(diagnostics::ExecutionDiagnostics) @@ -1414,20 +1423,36 @@ function log_danger_banner_execution(diagnostics::ExecutionDiagnostics, config:: end """ - run_analysis(adapter, config, prepared_table, sample_metadata; diagnostics, manifest) -> ExecutionResult - -Stub that records full manifest and diagnostics, hard-stop with DANGER banner on failures. - -- adapter: AnalysisAdapter (RAdapter or JuliaAdapter stub) -- config: AnalysisConfig.AnalysisConfig -- prepared_table: Matrix{Float64} from prepare_analysis_table -- sample_metadata: OrderedDict or nothing -- diagnostics: ExecutionDiagnostics from prepare_analysis_table -- manifest: ExecutionManifest from prepare_analysis_table - -Returns ExecutionResult stub with mock results, full manifest, diagnostics, provenance, hash chain. - -Hard-stop with DANGER banner on failures: if diagnostics has errors (should have been hard-stopped earlier in prepare_analysis_table) or if adapter method mismatches config method, or if prepared table has NaN/Inf after healing, throws ArgumentError with DANGER banner. + run_analysis(adapter, config, prepared_table; sample_metadata, diagnostics, manifest, + sample_ids, taxa_ids) -> ExecutionResult + +Run the model named by `config` and return its per-feature estimates with the full manifest, +diagnostics, provenance and hash chain. + +- adapter: AnalysisAdapter (RAdapter or JuliaAdapter). Its `method` must equal the config's. +- config: AnalysisConfig.AnalysisConfig — the declared method, formula and correction. +- prepared_table: Matrix{Float64} from `prepare_analysis_table`. +- sample_metadata: per-sample metadata for the design named by `config.formula`, as + column => per-sample vector. Without it no test is run and the result says so; the method + catalogue is explicit that a missing design means a descriptive summary, not a guess. +- diagnostics, manifest: from `prepare_analysis_table`. +- sample_ids, taxa_ids: labels, in the prepared table's order. + +Three outcomes, all explicit: + +- the fit ran: `results` holds one entry per feature, with `pvalue`/`padj` from the model and + `status` in {"ok", "boundary", "failed"} per feature; +- some features failed: those rows carry `status = "failed"`, a `note`, and `nothing` for + every statistic. They are excluded from the BH family and counted in + `diagnostics.checks["estimation"]["n_failed"]`; +- nothing ran (no design, R unavailable or busy): `results` is empty, the reason is in + `provenance["estimation"]["reason"]` and in a warning on the diagnostics. An empty table of + results is a statement about the run, and this function never dresses one up as a statement + about the data. + +Hard-stop with DANGER banner on failures: if diagnostics has errors, if the adapter method +mismatches the config method, if the prepared table still has NaN/Inf after healing, or if the +requested dispersion method has no implementation (refusing to substitute another one). """ function run_analysis( adapter::AnalysisAdapter, @@ -1489,37 +1514,43 @@ function run_analysis( # Log DANGER banner if dangerous log_danger_banner_execution(diagnostics, config) - # Mock results: for each taxon, create mock stats (p-value, log2FoldChange, etc.) - # Real implementation would call R via RCall or Julia via GLM, etc. - results = OrderedDict{String,Any}() - for (i, taxon_id) in enumerate(taxa_ids) - # Mock p-value and log2FoldChange - # For deterministic testing, use hash of taxon_id - h = hash(taxon_id) - p_val = 0.01 + (h % 100) / 1000.0 # 0.01..0.11 - log2fc = (h % 20) / 10.0 - 1.0 # -1..1 - results[taxon_id] = OrderedDict{String,Any}( - "pvalue" => p_val, - "padj" => p_val * 1.5, # mock BH - "log2FoldChange" => log2fc, - "baseMean" => 100.0 + (h % 1000), - "method" => config_method_str, - "adapter" => string(typeof(adapter)) - ) - end + # ---------------------------------------------------------------------- + # Estimation — every number below comes from a fit that ran + # + # This replaces a block that computed `p_val = 0.01 + (hash(taxon_id) % 100)/1000.0` + # and `padj = p_val * 1.5` and returned them as results. Those numbers were a + # function of the feature NAME: reproducible, plausible-looking and meaningless. + # See src/analysis/estimation.jl for what runs instead and + # docs/statistics/method-conditions/parametric-fits.md for what it supports. + # ---------------------------------------------------------------------- + outcome = Estimation.estimate_models( + config, prepared_table; + sample_metadata = sample_metadata, + offset = manifest.offset, + taxa_ids = taxa_ids, + sample_ids = sample_ids, + seed = adapter.seed, + r_wait_seconds = ESTIMATION_R_WAIT_SECONDS[] + ) - # If no taxa_ids provided, use generic - if isempty(results) - for i in 1:size(prepared_table, 1) - taxon_id = isempty(taxa_ids) ? "taxon_$i" : taxa_ids[i] - results[taxon_id] = OrderedDict{String,Any}( - "pvalue" => 0.05, - "padj" => 0.1, - "log2FoldChange" => 0.0, - "method" => config_method_str - ) - end + results = outcome.results + + # The estimation diagnostics travel with the execution diagnostics, and a run that did + # not happen is a warning on the record rather than an empty table of results. + checks = OrderedDict{String,Any}(diagnostics.checks) + checks["estimation"] = outcome.diagnostics + warnings = copy(diagnostics.warnings) + if outcome.status == :not_run + push!(warnings, "Estimation not run: $(outcome.reason)") end + exec_diagnostics = ExecutionDiagnostics( + warnings = warnings, + errors = diagnostics.errors, + healings = diagnostics.healings, + checks = checks, + is_dangerous = diagnostics.is_dangerous, + banner = diagnostics.banner + ) # Create ExecutionResult with full manifest and diagnostics exec_result = ExecutionResult( @@ -1528,7 +1559,7 @@ function run_analysis( config_hash=config.hash, method=config.method, results=results, - diagnostics=diagnostics, + diagnostics=exec_diagnostics, manifest=manifest, provenance=OrderedDict{String,Any}( "config_id" => config.id, @@ -1537,16 +1568,17 @@ function run_analysis( "manifest_hash" => manifest.hash, "adapter_type" => string(typeof(adapter)), "prepared_table_hash" => manifest.prepared_table_hash, + "estimation" => outcome.provenance, "diagnostics" => OrderedDict( - "warnings" => diagnostics.warnings, - "healings" => diagnostics.healings, - "checks" => diagnostics.checks, - "is_dangerous" => diagnostics.is_dangerous + "warnings" => exec_diagnostics.warnings, + "healings" => exec_diagnostics.healings, + "checks" => exec_diagnostics.checks, + "is_dangerous" => exec_diagnostics.is_dangerous ) ) ) - @info "run_analysis stub completed" result_id=exec_result.id config_id=config.id method=config_method_str taxa_count=length(results) dangerous=diagnostics.is_dangerous + @info "run_analysis finished" result_id=exec_result.id config_id=config.id method=config_method_str features=length(results) estimation_status=string(outcome.status) dangerous=exec_diagnostics.is_dangerous return exec_result end diff --git a/src/analysis/estimation.jl b/src/analysis/estimation.jl new file mode 100644 index 00000000..cebafcfe --- /dev/null +++ b/src/analysis/estimation.jl @@ -0,0 +1,849 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +# +# Catalogue item 2 of docs/statistics/method-catalogue-v1.md — parametric fits with +# explicit constraint, convergence and identifiability reporting. +# +# Conditions were published before this implementation, as the catalogue requires: +# docs/statistics/method-conditions/parametric-fits.md +# +# This module exists because the layer it replaces did not compute anything. The block +# it replaces in Execution.run_analysis derived its p-values from +# `p_val = 0.01 + (h % 100) / 1000.0` where `h = hash(taxon_id)` — a deterministic +# function of the feature NAME, carrying no information about the counts. The source +# called it a stub; the caller received it as a result. A stub that is returned as a +# result is a wrong answer, so it is gone rather than deprecated. +# +# The rule this module follows: a fit that cannot be performed produces NO number and +# a reason. There is no warn-and-carry-on path in this file that ends in a statistic. +# +# R is the estimator, deliberately. The pinned environment already carries R (RCall +# cannot build without it), renv.lock pins MASS, and the pipeline's other statistical +# work — NMDS, PERMANOVA, alpha significance — already runs there under one shared +# lock. Writing a negative binomial GLM by hand here would add a second, unvalidated +# implementation of a fitted model; docs/statistics/method-conditions/parametric-fits.md +# records why that trade was refused, and what this layer does and does not claim. + +module Estimation + +using Dates, Logging, OrderedCollections, SHA +using CSV +import JSON3 +using RCall +using ..RRuntime: with_r_lock, RBusyError +import ..AnalysisConfig + +export EstimationOutcome, estimate_models, bh_adjust, supported_terms, primary_term, + UnsupportedFormula, ESTIMATE_SCALE, REFUSED_DISPERSION, SUPPORTED_DISPERSION + +# --------------------------------------------------------------------------- +# What this module will not do +# --------------------------------------------------------------------------- + +""" + UnsupportedFormula(formula, detail) + +Thrown when a formula asks for a model this layer does not implement. Refusing is the +point: a term silently dropped from a design changes every estimate in the table, and a +term silently *added* would do the same. +""" +struct UnsupportedFormula <: Exception + formula :: String + detail :: String +end + +Base.showerror(io::IO, e::UnsupportedFormula) = + print(io, "Unsupported formula '", e.formula, "': ", e.detail) + +""" + REFUSED_DISPERSION + +Dispersion methods that are accepted as *configuration strings* — they are named in the +v1 catalogue and in issue #21 — but are not implemented here. Each is refused by name, +with a reason, rather than aliased to the one estimator this release implements. + +Why refusal rather than a fallback: dispersion is the parameter that sets every standard +error in a negative binomial fit. A run that asked for `glmGamPoi` and silently received +per-taxon maximum-likelihood theta would report different confidence intervals and +different multiple-testing outcomes under a label that says otherwise. See +`docs/statistics/method-conditions/parametric-fits.md`. +""" +const REFUSED_DISPERSION = Dict{String,String}( + "glmGamPoi" => "not implemented: the pinned R library (renv.lock) has no glmGamPoi, and the quasi-likelihood estimator that issue #21 asks for has not been written or validated here. Refusing the alias that v1 shipped.", + "local" => "not implemented: the trended (local) dispersion fit is not implemented. Refusing to substitute the per-taxon estimator.", + "mean" => "not implemented: the mean-dispersion fit is not implemented. Refusing to substitute the per-taxon estimator.", + "pooled" => "not implemented: the pooled single-theta fit is not implemented in this release. Refusing to substitute the per-taxon estimator.", +) + +const SUPPORTED_DISPERSION = ("parametric",) + +""" + ESTIMATE_SCALE + +What the `estimate` column is *in*, per method. A number whose units are unstated is not +a result, and "log2 fold change" is the wrong name for a difference in CLR units. +""" +const ESTIMATE_SCALE = Dict{AnalysisConfig.AnalysisMethod,String}( + AnalysisConfig.NB_GLM => "log_counts: the log-link coefficient, so log2FoldChange = estimate / log(2)", + AnalysisConfig.CLR_LM => "clr_difference: a difference in centred-log-ratio units, not a fold change", + AnalysisConfig.ILR_LM => "ilr_balance: a difference in isometric-log-ratio balance units, not a fold change", + AnalysisConfig.LOGISTIC => "log_odds: the logistic coefficient, so odds_ratio = exp(estimate)", +) + +# --------------------------------------------------------------------------- +# Formula: the supported subset, published and small +# --------------------------------------------------------------------------- + +const TERM_PATTERN = r"^[A-Za-z][A-Za-z0-9_.]*$" + +# Constructs with meaning that this layer does not implement. Refusing each by name is +# cheaper than explaining later why the design was not the one that was written down. +const REFUSED_CONSTRUCTS = Dict{String,String}( + "*" => "interaction terms ('*') are not implemented; write main effects only", + ":" => "interaction terms (':') are not implemented; write main effects only", + "|" => "random effects ('|') are not implemented", + "(" => "transformed terms such as I(x) or log(x) are not implemented; precompute the column", + ")" => "transformed terms such as I(x) or log(x) are not implemented; precompute the column", + "," => "functions of several arguments are not implemented", + "^" => "polynomial terms ('^') are not implemented; precompute the column", + "/" => "nested terms ('/') are not implemented; write main effects only", + "-" => "'-' both removes a term and appears inside column names; it is refused so that neither reading happens by accident", + "[" => "indexing is not implemented", + "]" => "indexing is not implemented", +) + +""" + supported_terms(formula) -> Vector{String} + +Parse the right-hand side of an R-style formula into the additive main effects this layer +supports, or throw `UnsupportedFormula`. + +Supported: `~ a`, `~ a + b`, `~ a + b + c`, where every term is a plain column name +(`^[A-Za-z][A-Za-z0-9_.]*`). Everything else — interactions, transformations, random +effects, nesting, polynomial terms, functions — is refused with a reason naming the +construct, because the design matrix is the one part of a fit that must never be quietly +not what was asked for. +""" +function supported_terms(formula::AbstractString)::Vector{String} + text = strip(String(formula)) + occursin("~", text) || + throw(UnsupportedFormula(text, "an R-style formula needs a '~', e.g. '~ group'")) + parts = split(text, "~"; limit = 2) + rhs = strip(parts[2]) + isempty(rhs) && + throw(UnsupportedFormula(text, "the right-hand side is empty; name at least one metadata column")) + + for (token, why) in REFUSED_CONSTRUCTS + occursin(token, rhs) && throw(UnsupportedFormula(text, why)) + end + + # Defence in depth rather than duplication: AnalysisConfig already refuses these in a + # formula, and this is the last checkpoint before the text is pasted into an R call, so it + # refuses them too. A backtick or a newline would otherwise let a formula name something + # that is not a column at all. + for (token, why) in [ + ("\$" => "a '\$' cannot appear in a formula; column names are plain identifiers"), + ("`" => "backticks are not permitted in a formula; column names are plain identifiers"), + (";" => "';' would run more than one statement; a formula names a model, not a script"), + ("\n" => "a newline cannot appear in a formula"), + ("\\" => "a backslash cannot appear in a formula"), + ] + occursin(token, rhs) && throw(UnsupportedFormula(text, why)) + end + + terms = String[] + for raw in split(rhs, "+") + term = strip(raw) + isempty(term) && + throw(UnsupportedFormula(text, "there is a '+' with no term beside it")) + occursin(TERM_PATTERN, term) || + throw(UnsupportedFormula(text, "term '$term' is not a plain column name (letters, digits, '_', '.'); see context_help('formula')")) + push!(terms, term) + end + + length(unique(terms)) == length(terms) || + throw(UnsupportedFormula(text, "the same term appears more than once")) + + return terms +end + +""" + primary_term(formula) -> String + +The first term on the right-hand side. It is the coefficient the per-feature table reports, +and the coefficient the Benjamini-Hochberg family is built from; every other coefficient is +reported beside it in `all_coefficients`, and the contrast it stands for is recorded in the +provenance. This is a stated convention, not a default that varies with the data: the first +term is the contrast of interest because the caller wrote it first. +""" +primary_term(formula::AbstractString) = first(supported_terms(formula)) + +# --------------------------------------------------------------------------- +# Multiple testing: BH, implemented here and checked against R's p.adjust +# --------------------------------------------------------------------------- + +""" + bh_adjust(pvalues) -> Vector{Float64} + +Benjamini-Hochberg step-up adjusted p-values, with the monotonicity constraint enforced. A +p-value that is not finite and in `[0, 1]` is a *failed fit*, not a number to adjust: it +raises `ArgumentError` so the caller has to surface it. + +`test/unit/test_estimation.jl` compares this against R's `p.adjust(method="BH")` rather than +against a second implementation of the same reasoning. +""" +function bh_adjust(pvalues::AbstractVector{<:Real})::Vector{Float64} + n = length(pvalues) + n == 0 && return Float64[] + for (i, p) in enumerate(pvalues) + if !(isfinite(p) && 0 <= p <= 1) + throw(ArgumentError("bh_adjust received a p-value that is not a probability at position $i ($p). A fit that did not produce a testable p-value is an unsuccessful state and must be reported as one, not adjusted.")) + end + end + + order = sortperm(collect(pvalues)) + sorted = Float64[Float64(pvalues[i]) for i in order] + adjusted = Vector{Float64}(undef, n) + running = Inf + for i in n:-1:1 + running = min(running, sorted[i] * n / i) + adjusted[i] = min(running, 1.0) + end + + out = Vector{Float64}(undef, n) + for (rank, idx) in enumerate(order) + out[idx] = adjusted[rank] + end + return out +end + +# --------------------------------------------------------------------------- +# Outcome +# --------------------------------------------------------------------------- + +""" + EstimationOutcome + +`status` is one of: + +- `:ok` — every feature in the table produced a testable fit +- `:partial` — some features did not; each one says so in its own row, and the count is in + `diagnostics["n_failed"]` +- `:not_run` — no fit was attempted, because the design or the runtime was not available. + `reason` says which, and `results` is empty. + +There is no fourth status meaning "something plausible was produced". +""" +struct EstimationOutcome + status :: Symbol + reason :: String + method :: String + results :: OrderedDict{String,Any} + diagnostics :: OrderedDict{String,Any} + provenance :: OrderedDict{String,Any} +end + +function _not_run(method::String, reason::String; err = nothing)::EstimationOutcome + @warn "Analysis not run — no statistics were produced" method reason + diagnostics = OrderedDict{String,Any}("status" => "not_run", "reason" => reason) + isnothing(err) || (diagnostics["exception_type"] = string(typeof(err))) + results = OrderedDict{String,Any}() + provenance = OrderedDict{String,Any}("status" => "not_run", "reason" => reason, + "method" => method) + return EstimationOutcome(:not_run, reason, method, results, diagnostics, provenance) +end + +# `nothing` rather than NaN: a JSON null is visibly absent, while a NaN that reaches a +# browser table renders as "#N/A" or as 0 depending on the reader. Absent must look absent. +_num(x)::Union{Float64,Nothing} = (ismissing(x) || !isfinite(x)) ? nothing : Float64(x) +_str(x)::String = ismissing(x) ? "" : String(x) +_bool(x)::Bool = x === true || (x isa AbstractString && uppercase(String(x)) == "TRUE") + +# --------------------------------------------------------------------------- +# Design +# --------------------------------------------------------------------------- + +""" + _column_vectors(sample_metadata) -> Dict{String,Vector} + +Accept the shapes callers actually have to hand: a dictionary of column name to per-sample +vector, or a vector of per-sample dictionaries. Anything else is refused rather than guessed +at, and a column with a missing value is refused too — a fit cannot be run on a design that +is partly absent without deciding what the absent entries were, and this layer does not make +that decision for the caller. +""" +function _column_vectors(sample_metadata)::Dict{String,Vector} + cols = Dict{String,Vector}() + + if sample_metadata isa AbstractDict + for (k, v) in sample_metadata + name = String(k) + v isa AbstractVector || + throw(ArgumentError("sample_metadata['$name'] is a $(typeof(v)); this layer needs one value per sample, as a vector, e.g. \"group\" => [\"a\", \"a\", \"b\", \"b\"].")) + values = collect(v) + if any(x -> x === missing || x === nothing || (x isa AbstractString && isempty(x)), values) + throw(ArgumentError("sample_metadata['$name'] has missing value(s). A design with holes in it is not the design the formula describes; refusing rather than dropping samples silently.")) + end + cols[name] = values + end + return cols + end + + if sample_metadata isa AbstractVector + for row in sample_metadata + row isa AbstractDict || + throw(ArgumentError("sample_metadata rows must be dictionaries keyed by column name; got $(typeof(row)).")) + for (k, v) in row + push!(get!(cols, String(k), Vector{Any}()), v) + end + end + return cols + end + + throw(ArgumentError("sample_metadata must be a dictionary of column => per-sample vector, or a vector of per-sample dictionaries; got $(typeof(sample_metadata)).")) +end + +# RCall needs a concrete vector type to cross the boundary; a Vector{Any} of strings is not +# one. Non-numeric columns travel as character and are made factors in R. +function _r_vector(values::Vector)::Union{Vector{String},Vector{Float64}} + if all(v -> v isa Real && !(v isa Bool), values) + return Float64[Float64(v) for v in values] + end + return String[string(v) for v in values] +end + +# --------------------------------------------------------------------------- +# Fitting +# --------------------------------------------------------------------------- + +""" + estimate_models(config, prepared; sample_metadata, offset, taxa_ids, sample_ids, + seed, r_wait_seconds) -> EstimationOutcome + +Fit the model named by `config.method` to every row of `prepared`, with the design named by +`config.formula`, and return per-feature estimates. + +What is supported, and what happens otherwise: + +- `formula` — additive main effects over metadata columns. Anything else is refused with + `UnsupportedFormula`. +- `nb_glm` — `MASS::glm.nb` per feature, with an offset. An offset is **required**: comparing + counts across samples without one compares library sizes. +- `clr_lm`, `ilr_lm` — `stats::lm` per feature on the transformed table. An offset is refused + here, because an offset has no meaning for a Gaussian fit on CLR units. +- `logistic` — `stats::glm(family = binomial)` on a 0/1 response. A non-binary response is + refused: a binomial fit on proportions needs the number of trials, and this layer will not + invent weights. +- `advanced.dispersion_method` — `parametric` only. Everything else is refused by name; see + `REFUSED_DISPERSION`. + +Unsuccessful states are reported, not smoothed over: a feature whose fit fails gets +`status = "failed"`, a `note` saying why, and `nothing` for every statistic. It is excluded +from the BH family, and the number of exclusions is reported — an exclusion rate high enough +to matter is a finding about the data, not a detail. +""" +function estimate_models(config::AnalysisConfig.AnalysisConfig, + prepared::Matrix{Float64}; + sample_metadata = nothing, + offset::Union{Vector{Float64},Nothing} = nothing, + taxa_ids::Vector{String} = String[], + sample_ids::Vector{String} = String[], + seed::Integer = 0, + r_wait_seconds::Real = 10.0)::EstimationOutcome + + method = config.method + method_str = AnalysisConfig.METHOD_TO_STRING[method] + n_features, n_samples = size(prepared) + + # -- the design ------------------------------------------------------ + terms = supported_terms(config.formula) # throws on unsupported syntax + prim = first(terms) + + if isnothing(sample_metadata) + return _not_run(method_str, + "no sample metadata was supplied, so the design named by '$(config.formula)' cannot be built. The method catalogue is explicit that a missing design means a descriptive summary: this layer will not guess a grouping. Supplying metadata (column => value per sample) is what turns this into a test.") + end + + cols = _column_vectors(sample_metadata) + missing_cols = [t for t in terms if !haskey(cols, t)] + if !isempty(missing_cols) + return _not_run(method_str, + "the formula names $(join(missing_cols, ", ")) but the supplied metadata has only $(join(sort(collect(keys(cols))), ", ")). Refusing rather than dropping a term: a design that is not the one written down is a different study.") + end + for t in terms + length(cols[t]) == n_samples || + return _not_run(method_str, + "metadata column '$t' has $(length(cols[t])) values for $n_samples samples. Refusing.") + end + + primary_levels = unique(cols[prim]) + if length(primary_levels) < 2 + return _not_run(method_str, + "the primary term '$prim' has $(length(primary_levels)) distinct value(s) across the samples, so there is no contrast to estimate. This is a property of the data, not a failure of the fit.") + end + + # -- method preconditions (configuration errors: hard refusals) ------- + if method == AnalysisConfig.NB_GLM + dispersion = lowercase(strip(config.advanced.dispersion_method)) + haskey(REFUSED_DISPERSION, dispersion) && + throw(ArgumentError("dispersion_method '$dispersion' is $(REFUSED_DISPERSION[dispersion]) See issue #21 and docs/statistics/method-conditions/parametric-fits.md.")) + dispersion in SUPPORTED_DISPERSION || + throw(ArgumentError("dispersion_method '$dispersion' is not one of $(join(SUPPORTED_DISPERSION, ", ")). Refusing.")) + + isnothing(offset) && + throw(ArgumentError("nb_glm needs an offset: raw counts are not comparable across samples, and a negative binomial fit without library-size or size-factor offsets would attribute sequencing depth to biology. prepare_analysis_table supplies log library size for 'none'/'TSS'/'CSS'/'RSS' and log size factors for 'size_factors'. Refusing.")) + length(offset) == n_samples || + throw(ArgumentError("the offset has $(length(offset)) values for $n_samples samples. Refusing.")) + all(isfinite, offset) || + throw(ArgumentError("the offset contains non-finite values; it is a log scale factor, so a non-finite entry means the scaling failed upstream. Refusing.")) + else + isnothing(offset) || + throw(ArgumentError("an offset is only defined for count models (nb_glm); method '$method_str' is not one. Refusing to accept and then ignore a parameter that means something to a different method.")) + end + + if method == AnalysisConfig.LOGISTIC + binary = all(v -> v == 0.0 || v == 1.0, prepared) + binary || + throw(ArgumentError("logistic regression needs a 0/1 response and this prepared table holds values that are not 0 or 1. Use normalization 'presence_absence': a binomial fit on proportions needs the number of trials per sample, and this layer will not supply invented weights. Refusing.")) + end + + # -- row labels ------------------------------------------------------ + row_labels, label_note = _row_labels(method, n_features, taxa_ids) + + # -- the fit, in R --------------------------------------------------- + dir = mktempdir() + fits_path = joinpath(dir, "fits.csv") + coefs_path = joinpath(dir, "coefficients.csv") + versions = Ref{Dict{String,String}}(Dict{String,String}()) + fit_rows = Ref{Vector}(Any[]) + coef_rows = Ref{Vector}(Any[]) + fits_hash = Ref{String}("") + coefs_hash = Ref{String}("") + + try + with_r_lock(; timeout = r_wait_seconds) do + RCall.reval("est_cols <- list()") + for (i, t) in enumerate(terms) + RCall.globalEnv[Symbol("est_c$i")] = _r_vector(cols[t]) + RCall.reval("est_cols[[\"$t\"]] <- est_c$i") + end + if isnothing(offset) + RCall.reval("est_offset <- NULL") + else + RCall.globalEnv[:est_offset] = collect(Float64, offset) + end + RCall.globalEnv[:est_counts] = prepared + RCall.globalEnv[:est_n] = n_features + RCall.globalEnv[:est_ns] = n_samples + RCall.globalEnv[:est_rhs] = join(terms, " + ") + RCall.globalEnv[:est_primary] = prim + RCall.globalEnv[:est_method] = method_str + RCall.globalEnv[:est_dispersion] = lowercase(strip(config.advanced.dispersion_method)) + RCall.globalEnv[:est_fits_path] = fits_path + RCall.globalEnv[:est_coefs_path] = coefs_path + + RCall.reval(R_ESTIMATION_SETUP) + RCall.reval(R_ESTIMATION_FIT) + + versions[] = Dict( + "r" => RCall.rcopy(String, RCall.reval("est_r_version")), + "mass" => RCall.rcopy(String, RCall.reval("est_mass_version")), + "primary_levels" => RCall.rcopy(String, RCall.reval("est_primary_levels_text")), + ) + fit_rows[] = collect(CSV.File(fits_path)) + coef_rows[] = collect(CSV.File(coefs_path)) + RCall.reval("rm(list = ls(pattern = \"^est_\")); gc()") + return nothing + end + + fits_hash[] = bytes2hex(sha256(read(fits_path))) + coefs_hash[] = bytes2hex(sha256(read(coefs_path))) + + return _assemble(config, method, method_str, terms, prim, row_labels, label_note, + fit_rows[], coef_rows[], offset, versions[], fits_hash[], coefs_hash[], + seed, n_features, n_samples) + catch err + # A runtime that cannot produce a fit produces no number. RBusyError is the + # documented case (a pipeline stage holds the R lock); anything else is reported + # with its own exception type so it can be diagnosed rather than guessed at. + reason = err isa RBusyError ? + "the embedded R runtime was busy for $(r_wait_seconds)s — another task held it. No statistics were produced; retry when the pipeline run has finished." : + "the fit could not be run: $(sprint(showerror, err))" + return _not_run(method_str, reason; err = err) + finally + rm(dir; recursive = true, force = true) + end +end + +""" + _row_labels(method, n_features, taxa_ids) -> (Vector{String}, String) + +Label rows by feature id when the ids line up with the rows. For ILR the prepared table has +one row per balance, which is not one per taxon, so those rows are named as balances and the +mismatch is recorded rather than papered over with a taxon name that does not describe them. +""" +function _row_labels(method::AnalysisConfig.AnalysisMethod, n_features::Int, + taxa_ids::Vector{String}) + if length(taxa_ids) == n_features && !isempty(taxa_ids) + return (copy(taxa_ids), "row labels are the feature ids supplied by the caller") + end + if method == AnalysisConfig.ILR_LM + return (["balance_$i" for i in 1:n_features], + "row labels are balances: an ILR table has one row per balance, not one per taxon, so no taxon id is attached to these numbers") + end + return (["feature_$i" for i in 1:n_features], + "no feature ids were supplied; rows are numbered in the order of the prepared table") +end + +# --------------------------------------------------------------------------- +# Assembly: R's output -> one outcome, with the BH family stated +# --------------------------------------------------------------------------- + +function _assemble(config, method, method_str, terms, prim, row_labels, label_note, + fit_rows::Vector, coef_rows::Vector, offset, versions, + fits_hash, coefs_hash, seed, n_features, n_samples) + by_feature = Dict{Int,Vector{Any}}() + for row in coef_rows + push!(get!(by_feature, Int(row.taxon_index), Any[]), row) + end + + statuses = [_str(row.status) for row in fit_rows] + n_failed = count(==("failed"), statuses) + n_boundary = count(==("boundary"), statuses) + + # The BH family is the set of features whose fit produced a p-value. Exclusions are + # counted and reported; they are never replaced by 1.0, which would look like a negative + # result rather than like a fit that did not happen. + testable = Int[] + raw_p = Float64[] + for (i, row) in enumerate(fit_rows) + statuses[i] == "failed" && continue + p = _num(row.pvalue) + isnothing(p) && continue + push!(testable, i) + push!(raw_p, p) + end + + bh_disabled = config.correction.allow_no_correction + adjusted = (bh_disabled || isempty(raw_p)) ? Float64[] : bh_adjust(raw_p) + + results = OrderedDict{String,Any}() + for (i, row) in enumerate(fit_rows) + status = statuses[i] + p_raw = _num(row.pvalue) + padj = nothing + if !bh_disabled && !isnothing(p_raw) && !isempty(adjusted) + pos = findfirst(==(i), testable) + isnothing(pos) || (padj = _num(adjusted[pos])) + end + + estimate = _num(row.estimate) + theta = _num(row.dispersion_theta) + + all_coefs = OrderedDict{String,Any}() + for c in get(by_feature, i, Any[]) + all_coefs[_str(c.term)] = OrderedDict{String,Any}( + "estimate" => _num(c.estimate), + "standard_error" => _num(c.standard_error), + "statistic" => _num(c.statistic), + "pvalue" => _num(c.pvalue), + "is_primary" => _bool(c.is_primary), + ) + end + + entry = OrderedDict{String,Any}( + "status" => status, + "note" => _str(row.note), + "estimate" => estimate, + "estimate_scale" => ESTIMATE_SCALE[method], + "standard_error" => _num(row.standard_error), + "statistic" => _num(row.statistic), + "statistic_name" => _str(row.statistic_name), + "pvalue" => p_raw, + "padj" => padj, + "ci_low" => _num(row.ci_low), + "ci_high" => _num(row.ci_high), + "ci_method" => "Wald, 95%", + "residual_df" => _num(row.residual_df), + "n_observations" => _num(row.n_observations), + "primary_term" => prim, + "all_coefficients" => all_coefs, + ) + if method == AnalysisConfig.NB_GLM + entry["log2FoldChange"] = isnothing(estimate) ? nothing : estimate / log(2) + entry["dispersion_theta"] = theta + elseif method == AnalysisConfig.LOGISTIC + entry["odds_ratio"] = isnothing(estimate) ? nothing : exp(estimate) + end + results[row_labels[i]] = entry + end + + overall = n_failed == 0 ? :ok : :partial + reason = n_failed == 0 ? "" : + "$n_failed of $n_features features did not produce a testable fit; each one says why in its own row. They are excluded from the multiple-testing family, and the exclusions are counted here rather than replaced by a p-value of 1." + + offset_kind = if isnothing(offset) + "none" + elseif method == AnalysisConfig.NB_GLM && lowercase(strip(config.normalization.method)) == "size_factors" + "log(size factor) per sample (median-of-ratios, from the prepared manifest)" + elseif method == AnalysisConfig.NB_GLM + "log(library size) per sample (from the prepared manifest)" + else + "none" + end + + diagnostics = OrderedDict{String,Any}( + "status" => string(overall), + "reason" => reason, + "method" => method_str, + "n_features" => n_features, + "n_samples" => n_samples, + "n_tested" => length(testable), + "n_failed" => n_failed, + "n_boundary" => n_boundary, + "family" => bh_disabled ? "BH explicitly disabled by the acknowledging caller (DANGER)" : + "BH (Benjamini-Hochberg) over the $(length(testable)) features with a testable fit", + "alpha" => config.correction.alpha, + "terms" => terms, + "primary_term" => prim, + "offset" => offset_kind, + "row_labels" => label_note, + "failure_note" => n_failed == 0 ? "no failures" : + "features whose fit failed are excluded from the BH family, which changes what the family is: if a large share failed, the adjusted p-values describe the features that fitted, and whoever reports them has to say so.", + ) + + provenance = OrderedDict{String,Any}( + "status" => string(overall), + "method" => method_str, + "engine" => "R through RCall, serialised on the shared runtime lock", + "formula" => config.formula, + "terms" => terms, + "primary_term" => prim, + "primary_contrast" => "coefficient of the first term; for a factor this is the contrast between the first two levels of $(get(versions, "primary_levels", "")) (R's default treatment coding)", + "estimate_scale" => ESTIMATE_SCALE[method], + "dispersion_method" => method == AnalysisConfig.NB_GLM ? lowercase(strip(config.advanced.dispersion_method)) : "not applicable to $method_str", + "offset" => offset_kind, + "offset_sha256" => isnothing(offset) ? nothing : bytes2hex(sha256(join(string.(offset), ","))), + "correction" => bh_disabled ? "NONE (acknowledged override; DANGER)" : "BH", + "alpha" => config.correction.alpha, + "r_version" => get(versions, "r", "unknown"), + "mass_version" => get(versions, "mass", "unknown"), + "fits_csv_sha256" => fits_hash, + "coefficients_csv_sha256" => coefs_hash, + "randomness" => "none — these fits are deterministic (IRLS / least squares; no permutations or resampling). seed=$seed was supplied by the caller and is not a parameter of any number here.", + "interval_method" => "Wald (estimate ± 1.96·SE for z-based fits, t quantile for the Gaussian fits). Profile intervals are not computed: they cost orders of magnitude more per feature and would not change which features are reported.", + ) + + return EstimationOutcome(overall, reason, method_str, results, diagnostics, provenance) +end + +# --------------------------------------------------------------------------- +# The R that does the work +# +# Held to the repository's boundary rule: values cross through CSV rather than as R objects, +# so an NA arrives in Julia as a missing value and never as a zero or a NaN that a downstream +# reader would render as a number. R list elements are reached with [[ ]] rather than with $, +# which keeps the embeddable-R lint (config/ci/lint_source.jl, check 2) honest instead of +# having to exempt this file. +# --------------------------------------------------------------------------- + +const R_ESTIMATION_SETUP = """ + suppressPackageStartupMessages(library(MASS)) + options(stringsAsFactors = FALSE) + + est_meta <- as.data.frame(est_cols, stringsAsFactors = FALSE, check.names = FALSE) + if (nrow(est_meta) != est_ns) { + stop("metadata has ", nrow(est_meta), " rows for ", est_ns, " samples") + } + for (nm in names(est_meta)) { + if (is.character(est_meta[[nm]])) est_meta[[nm]] <- factor(est_meta[[nm]]) + } + if (is.factor(est_meta[[est_primary]])) { + est_primary_levels_text <- paste(levels(est_meta[[est_primary]]), collapse = " then ") + } else { + est_primary_levels_text <- paste0("a numeric column: ", est_primary) + } + + est_r_version <- paste(R.version[["major"]], R.version[["minor"]], sep = ".") + est_mass_version <- tryCatch(as.character(utils::packageVersion("MASS")), + error = function(e) "unavailable") + + est_form <- if (est_method == "nb_glm") { + as.formula(paste("y ~", est_rhs, "+ offset(off)")) + } else { + as.formula(paste("y ~", est_rhs)) + } +""" + +const R_ESTIMATION_FIT = """ + est_status <- character(est_n) + est_note <- character(est_n) + est_estimate <- rep(NA_real_, est_n) + est_se <- rep(NA_real_, est_n) + est_stat <- rep(NA_real_, est_n) + est_stat_name <- rep(NA_character_, est_n) + est_p <- rep(NA_real_, est_n) + est_df <- rep(NA_real_, est_n) + est_theta <- rep(NA_real_, est_n) + est_nobs <- rep(NA_real_, est_n) + est_ci_low <- rep(NA_real_, est_n) + est_ci_high <- rep(NA_real_, est_n) + est_meta[["y"]] <- rep(0, est_ns) + if (!is.null(est_offset)) est_meta[["off"]] <- as.numeric(est_offset) + + est_one_fit <- function() { + if (est_method == "nb_glm") { + if (est_dispersion != "parametric") { + stop("dispersion method '", est_dispersion, "' has no implementation here") + } + return(MASS::glm.nb(est_form, data = est_meta, control = glm.control(maxit = 100))) + } else if (est_method == "logistic") { + return(stats::glm(est_form, data = est_meta, family = stats::binomial())) + } + return(stats::lm(est_form, data = est_meta)) + } + + est_coef_rows <- list() + est_k <- 0 + + for (i in seq_len(est_n)) { + y <- as.numeric(est_counts[i, ]) + est_nobs[i] <- sum(is.finite(y)) + if (sum(is.finite(y)) < 3 || length(unique(y[is.finite(y)])) < 2) { + est_status[i] <- "failed" + est_note[i] <- "the response is constant or has fewer than three finite values: there is nothing to estimate" + next + } + est_meta[["y"]] <- y + warns <- character(0) + fit <- tryCatch( + withCallingHandlers(est_one_fit(), warning = function(w) { + warns <<- c(warns, conditionMessage(w)) + invokeRestart("muffleWarning") + }), + error = function(e) e + ) + if (inherits(fit, "error")) { + est_status[i] <- "failed" + est_note[i] <- conditionMessage(fit) + next + } + cf <- tryCatch(summary(fit)[["coefficients"]], error = function(e) NULL) + if (is.null(cf) || is.null(dim(cf))) { + est_status[i] <- "failed" + est_note[i] <- "the fit produced no coefficient table" + next + } + coef_name <- est_primary + if (!(coef_name %in% rownames(cf))) { + cand <- setdiff(rownames(cf), "(Intercept)") + cand <- cand[startsWith(cand, est_primary)] + if (length(cand) == 0) { + est_status[i] <- "failed" + est_note[i] <- paste("the fit has no coefficient for the primary term", est_primary) + next + } + coef_name <- cand[1] + } + + status <- "ok" + notes <- character(0) + if (est_method == "nb_glm") { + if (!isTRUE(fit[["converged"]])) { + status <- "failed" + notes <- c(notes, "the fit did not converge") + } else if (fit[["theta"]] >= 1e7) { + status <- "boundary" + notes <- c(notes, "theta >= 1e7: dispersion is indistinguishable from Poisson at this precision, so the negative binomial model adds nothing here") + } else if (fit[["theta"]] <= 1e-8) { + status <- "boundary" + notes <- c(notes, "theta <= 1e-8: the dispersion estimate is at its boundary and the standard error of this coefficient is not trustworthy") + } + est_theta[i] <- fit[["theta"]] + } + if (est_method == "logistic" && any(grepl("numerically 0 or 1", warns))) { + status <- "boundary" + notes <- c(notes, "complete or quasi-complete separation: fitted probabilities reached 0 or 1, so this Wald test is not a valid p-value") + } + if (any(grepl("iteration limit|alternation limit|did not converge|NaNs produced", warns))) { + status <- "failed" + notes <- c(notes, paste(unique(warns), collapse = "; ")) + } + + est_stat_name[i] <- if (est_method == "logistic" || est_method == "nb_glm") "z" else "t" + est_df[i] <- fit[["df.residual"]] + est_note[i] <- paste(c(notes, unique(warns)), collapse = "; ") + + if (status == "failed") { + est_status[i] <- "failed" + next + } + + est_estimate[i] <- cf[coef_name, 1] + est_se[i] <- cf[coef_name, 2] + est_stat[i] <- cf[coef_name, 3] + est_p[i] <- cf[coef_name, 4] + + q <- if (est_method == "logistic" || est_method == "nb_glm") stats::qnorm(0.975) else stats::qt(0.975, fit[["df.residual"]]) + est_ci_low[i] <- cf[coef_name, 1] - q * cf[coef_name, 2] + est_ci_high[i] <- cf[coef_name, 1] + q * cf[coef_name, 2] + + if (!is.finite(est_p[i]) || est_p[i] < 0 || est_p[i] > 1) { + est_status[i] <- "failed" + est_note[i] <- paste(est_note[i], "the fit produced a p-value that is not a probability; refused rather than reported") + est_estimate[i] <- NA_real_ + est_se[i] <- NA_real_ + est_stat[i] <- NA_real_ + est_p[i] <- NA_real_ + est_ci_low[i] <- NA_real_ + est_ci_high[i] <- NA_real_ + next + } + est_status[i] <- status + + for (term in setdiff(rownames(cf), "(Intercept)")) { + est_k <- est_k + 1 + est_coef_rows[[est_k]] <- data.frame( + taxon_index = i, + term = term, + estimate = cf[term, 1], + standard_error = cf[term, 2], + statistic = cf[term, 3], + pvalue = cf[term, 4], + is_primary = identical(term, coef_name), + stringsAsFactors = FALSE + ) + } + } + + est_fits <- data.frame( + taxon_index = seq_len(est_n), + status = est_status, + note = est_note, + estimate = est_estimate, + standard_error = est_se, + statistic = est_stat, + statistic_name = est_stat_name, + pvalue = est_p, + residual_df = est_df, + dispersion_theta = est_theta, + n_observations = est_nobs, + ci_low = est_ci_low, + ci_high = est_ci_high, + stringsAsFactors = FALSE + ) + utils::write.csv(est_fits, est_fits_path, row.names = FALSE, na = "NA") + + est_coefs <- if (est_k > 0) { + do.call(rbind, est_coef_rows) + } else { + data.frame( + taxon_index = integer(0), term = character(0), estimate = numeric(0), + standard_error = numeric(0), statistic = numeric(0), pvalue = numeric(0), + is_primary = logical(0), stringsAsFactors = FALSE + ) + } + utils::write.csv(est_coefs, est_coefs_path, row.names = FALSE, na = "NA") +""" + +end # module Estimation diff --git a/test/runtests.jl b/test/runtests.jl index 5305f1f4..a1e4f1ea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,7 @@ using MetaManifold: AnalysisConfig include("unit/test_migrate_composition.jl") include("unit/test_analysis_config.jl") include("unit/test_execution.jl") + include("unit/test_estimation.jl") ## Integration tests (opt-in) if RUN_INTEGRATION diff --git a/test/unit/test_estimation.jl b/test/unit/test_estimation.jl new file mode 100644 index 00000000..acc71d39 --- /dev/null +++ b/test/unit/test_estimation.jl @@ -0,0 +1,435 @@ +# SPDX-License-Identifier: AGPL-3.0-only +# Evidence for catalogue item 2 of docs/statistics/method-catalogue-v1.md, held to the +# conditions published before the implementation existed in +# docs/statistics/method-conditions/parametric-fits.md: +# +# * known answers, where the true effect is written into the data rather than read back +# out of the fit; +# * an independent reference — the same data fitted directly in R, compared coefficient by +# coefficient, and BH compared against R's p.adjust; +# * negative controls — refusals that must throw, and unsuccessful states that must come +# back as unsuccessful rather than as numbers; +# * a source-level guard that the placeholder statistics this layer replaced do not return. + +@testset "Estimation — parametric fits (catalogue item 2)" begin + + using OrderedCollections + using RCall + using MetaManifold.Estimation + using MetaManifold.RRuntime: with_r_lock + + # ------------------------------------------------------------------ + # Synthetic data with a written-in effect + # + # Deterministic on purpose: a fixed table cannot drift under a changed RNG between Julia + # versions, and the expected log2 fold changes below are arithmetic on the means rather + # than a recording of what the fit happened to say. + # ------------------------------------------------------------------ + est_groups = vcat(fill("A", 12), fill("B", 12)) + + # taxon 1: mean 10 vs 40 (four-fold up) -> log2FC ~ +2 + # taxon 2: mean 30 vs 30 (no effect) -> log2FC ~ 0 + # taxon 3: mean 50 vs 5 (ten-fold down) -> log2FC ~ -3.3 + counts_effect = [ + 8.0 12 9 11 10 13 7 10 12 9 11 10 38 42 40 36 45 41 39 44 37 43 40 41; + 30.0 28 32 31 29 33 30 27 31 29 32 30 31 29 30 32 28 31 30 32 29 31 30 30; + 50.0 52 48 51 49 53 50 47 51 49 52 50 5 6 4 5 7 5 6 4 5 6 5 5 + ] + taxa_effect = ["t_up", "t_flat", "t_down"] + offsets_effect = log.(vec(sum(counts_effect, dims = 1))) + + meta_effect = OrderedDict{String,Any}("group" => est_groups) + + function nb_config(; dispersion = "parametric", normalization = "size_factors", + min_samples = 2) + norm = AnalysisConfig.NormalizationConfig(method = normalization) + adv = AnalysisConfig.AdvancedConfig(min_samples_per_group = min_samples, + dispersion_method = dispersion) + AnalysisConfig.AnalysisConfig( + method = "nb_glm", + formula = "~ group", + metadata_columns = ["group"], + normalization = norm, + advanced = adv, + created_by = "test_estimation", + ) + end + + pv(entry) = entry["pvalue"] === nothing ? NaN : Float64(entry["pvalue"]) + pa(entry) = entry["padj"] === nothing ? NaN : Float64(entry["padj"]) + rel_err(a, b) = abs(a - b) / max(abs(b), 1e-300) + + @testset "Benjamini-Hochberg, against answers and against R" begin + # Hand answer: BH with n = 4 on p = (0.01, 0.02, 0.03, 0.04). + # raw rank-adjusted: 0.04, 0.04, 0.04, 0.04 -> all 0.04. + @test Estimation.bh_adjust([0.01, 0.02, 0.03, 0.04]) == [0.04, 0.04, 0.04, 0.04] + + # Hand answer with a monotonicity violation to enforce: p = (0.001, 0.5, 0.5) gives + # 0.003, 0.75, 0.5 before the step-up constraint, and 0.003, 0.5, 0.5 after it, + # because an adjusted p-value can never be smaller than one with a smaller raw p. + adjusted = Estimation.bh_adjust([0.001, 0.5, 0.5]) + @test adjusted[1] ≈ 0.003 + @test adjusted[2] ≈ 0.5 + @test adjusted[3] ≈ 0.5 + @test issorted(adjusted) # monotone in the raw order + + # A p-value that is not a probability is a failed fit, not something to adjust. + @test_throws ArgumentError Estimation.bh_adjust([0.01, NaN, 0.2]) + @test_throws ArgumentError Estimation.bh_adjust([0.01, 1.4, 0.2]) + @test_throws ArgumentError Estimation.bh_adjust([0.01, -0.2, 0.2]) + @test Estimation.bh_adjust(Float64[]) == Float64[] + + # Independent reference: R's own p.adjust. + r_ok_bh = try + with_r_lock(; timeout = 10.0) do + RCall.reval("t_bh_p <- c(0.001, 0.5, 0.5, 0.02, 0.7, 0.04, 0.0009, 0.31)") + RCall.rcopy(RCall.reval("p.adjust(t_bh_p, method = \"BH\")")) + end + catch + nothing + end + if r_ok_bh === nothing + @info "SKIPPED (by name): BH comparison against R p.adjust — R is not reachable here." + @test_skip Estimation.bh_adjust([0.001, 0.5, 0.5, 0.02, 0.7, 0.04, 0.0009, 0.31]) + else + mine = Estimation.bh_adjust([0.001, 0.5, 0.5, 0.02, 0.7, 0.04, 0.0009, 0.31]) + @test all(rel_err(mine[i], r_ok_bh[i]) < 1e-12 for i in eachindex(mine)) + end + end + + @testset "the formula subset is small and refusal is by name" begin + @test Estimation.supported_terms("~ group") == ["group"] + @test Estimation.supported_terms("~ group + batch") == ["group", "batch"] + @test Estimation.primary_term("~ batch + group") == "batch" + @test Estimation.supported_terms("disease ~ group + age") == ["group", "age"] + + for (formula, why) in [ + ("~ group * batch" => "interaction"), + ("~ group:batch" => "interaction"), + ("~ (1|batch)" => "random effect"), + ("~ I(age^2)" => "transform"), + ("~ log(abundance)" => "transform"), + ("~ poly(age, 2)" => "transform"), + ("~ group - batch" => "minus"), + ("~ group + group" => "duplicate"), + ("~ group +" => "dangling plus"), + ("group" => "missing tilde"), + (raw"~ group$bad" => "dollar"), + (raw"~ `group`" => "backtick"), + ] + # The message has to say which construct was refused: a refusal the caller cannot + # act on is only slightly better than a silent change of design. + @test_throws Estimation.UnsupportedFormula Estimation.supported_terms(formula) + end + err = try + Estimation.supported_terms("~ group * batch") + catch e + e + end + @test err isa Estimation.UnsupportedFormula + @test occursin("interaction", sprint(showerror, err)) + end + + @testset "a design that cannot be built is an unsuccessful state, not a fit" begin + # No metadata at all. + no_design = Estimation.estimate_models(nb_config(), counts_effect; + offset = offsets_effect, taxa_ids = taxa_effect) + @test no_design.status == :not_run + @test isempty(no_design.results) + @test occursin("descriptive summary", no_design.reason) + + # Metadata missing the term the formula names. + wrong_col = Estimation.estimate_models(nb_config(), counts_effect; + sample_metadata = OrderedDict{String,Any}("treatment" => est_groups), + offset = offsets_effect, taxa_ids = taxa_effect) + @test wrong_col.status == :not_run + @test occursin("group", wrong_col.reason) + + # Metadata of the wrong length. + short = Estimation.estimate_models(nb_config(), counts_effect; + sample_metadata = OrderedDict{String,Any}("group" => ["A", "B"]), + offset = offsets_effect, taxa_ids = taxa_effect) + @test short.status == :not_run + @test occursin("values for", short.reason) + + # A grouping column with one level has no contrast to estimate. + one_level = Estimation.estimate_models(nb_config(), counts_effect; + sample_metadata = OrderedDict{String,Any}("group" => fill("A", 24)), + offset = offsets_effect, taxa_ids = taxa_effect) + @test one_level.status == :not_run + @test occursin("distinct value", one_level.reason) + + # A hole in the design is refused rather than papered over. + meta_with_hole = OrderedDict{String,Any}( + "group" => vcat([x for x in fill("A", 12)], [missing for _ in 1:12])) + @test_throws ArgumentError Estimation.estimate_models( + nb_config(), counts_effect; sample_metadata = meta_with_hole, + offset = offsets_effect, taxa_ids = taxa_effect) + end + + @testset "unsupported requests are refused, never silently substituted" begin + # glmGamPoi is named in the v1 catalogue and in issue #21; it is not implemented, and + # the refusal has to say so rather than quietly return per-taxon theta. + err = try + Estimation.estimate_models(nb_config(dispersion = "glmGamPoi"), counts_effect; + sample_metadata = meta_effect, offset = offsets_effect, + taxa_ids = taxa_effect) + nothing + catch e + e + end + @test err isa ArgumentError + @test occursin("glmGamPoi", err.msg) + @test occursin("#21", err.msg) + + # Counting without an offset would compare library sizes instead of groups. + @test_throws ArgumentError Estimation.estimate_models( + nb_config(), counts_effect; sample_metadata = meta_effect, taxa_ids = taxa_effect) + + # An offset means nothing to a Gaussian fit on CLR units, so it is refused rather than + # accepted and ignored. + clr_norm = AnalysisConfig.NormalizationConfig(method = "clr") + clr_config = AnalysisConfig.AnalysisConfig( + method = "clr_lm", formula = "~ group", metadata_columns = ["group"], + normalization = clr_norm, created_by = "test_estimation_clr") + @test_throws ArgumentError Estimation.estimate_models( + clr_config, counts_effect; sample_metadata = meta_effect, + offset = offsets_effect, taxa_ids = taxa_effect) + + # A binomial fit on proportions has no number of trials unless one is invented. + logistic_norm = AnalysisConfig.NormalizationConfig(method = "relative") + logistic_config = AnalysisConfig.AnalysisConfig( + method = "logistic", formula = "~ group", outcome_column = "group", + metadata_columns = ["group"], normalization = logistic_norm, + created_by = "test_estimation_logistic") + @test_throws ArgumentError Estimation.estimate_models( + logistic_config, counts_effect; sample_metadata = meta_effect, taxa_ids = taxa_effect) + end + + # ------------------------------------------------------------------ + # Fits that need R. Skipped loudly and by name where R is absent, because a silent skip + # would turn a provisioning defect into a green run. + # ------------------------------------------------------------------ + r_ok = try + with_r_lock(; timeout = 10.0) do + RCall.reval("suppressPackageStartupMessages(library(MASS))") + true + end + catch + false + end + + if r_ok + @testset "a written-in effect is recovered, and matches a direct R fit" begin + outcome = Estimation.estimate_models(nb_config(), counts_effect; + sample_metadata = meta_effect, + offset = offsets_effect, + taxa_ids = taxa_effect) + + @test outcome.status == :ok + @test length(outcome.results) == 3 + @test outcome.diagnostics["n_tested"] == 3 + @test outcome.diagnostics["n_failed"] == 0 + @test outcome.provenance["status"] == "ok" + @test outcome.provenance["correction"] == "BH" + + up = outcome.results["t_up"] + flat = outcome.results["t_flat"] + down = outcome.results["t_down"] + + # Known answers: four-fold up, no effect, ten-fold down. + @test up["estimate"] ≈ log(4) atol = 0.5 + @test up["log2FoldChange"] ≈ 2.0 atol = 0.8 + @test down["estimate"] < 0 + @test abs(flat["estimate"]) < 0.6 + @test pv(up) < 0.01 + @test up["status"] == "ok" + + # Adjusted p-values are never smaller than raw ones, and never above 1. + for entry in values(outcome.results) + @test pa(entry) >= pv(entry) - 1e-12 + @test pa(entry) <= 1.0 + @test entry["estimate_scale"] isa String + @test occursin("log_counts", entry["estimate_scale"]) + end + + # The BH family is the three features here, computed by the same function the + # module used; comparing against it checks the wiring, not the arithmetic (which + # the p.adjust comparison above covers). + raw = [pv(outcome.results["t_up"]), pv(outcome.results["t_flat"]), pv(outcome.results["t_down"])] + expected = Estimation.bh_adjust(raw) + @test rel_err(pa(outcome.results["t_up"]), expected[1]) < 1e-12 + @test rel_err(pa(outcome.results["t_flat"]), expected[2]) < 1e-12 + @test rel_err(pa(outcome.results["t_down"]), expected[3]) < 1e-12 + + # Independent reference: the same fit, written directly in R. If the module wired + # the offset, the response or the contrast wrongly, these diverge. + reference = with_r_lock(; timeout = 30.0) do + RCall.globalEnv[:t_counts] = counts_effect + RCall.globalEnv[:t_group] = est_groups + RCall.globalEnv[:t_offset] = offsets_effect + RCall.reval(""" + t_group_f <- factor(t_group) + t_fit <- MASS::glm.nb(t_counts[1, ] ~ t_group_f + offset(t_offset)) + t_cf <- summary(t_fit)[["coefficients"]] + t_est <- t_cf[2, 1] + t_se <- t_cf[2, 2] + t_p <- t_cf[2, 4] + t_theta <- t_fit[["theta"]] + """) + ( + est = RCall.rcopy(RCall.reval("t_est")), + se = RCall.rcopy(RCall.reval("t_se")), + p = RCall.rcopy(RCall.reval("t_p")), + theta = RCall.rcopy(RCall.reval("t_theta")), + ) + end + @test rel_err(up["estimate"], reference.est) < 1e-6 + @test rel_err(up["standard_error"], reference.se) < 1e-6 + @test rel_err(pv(up), reference.p) < 1e-6 + @test rel_err(up["dispersion_theta"], reference.theta) < 1e-6 + + # The provenance has to be usable by someone else: which R, which MASS, which + # family, which offset, and the hashes of the tables the numbers came from. + for key in ("r_version", "mass_version", "family", "offset", "fits_csv_sha256", + "coefficients_csv_sha256", "primary_contrast", "randomness") + @test haskey(outcome.provenance, key) + end + @test length(outcome.provenance["fits_csv_sha256"]) == 64 + end + + @testset "a feature that cannot be fitted is reported, not invented" begin + counts_with_constant = vcat(counts_effect, permutedims(fill(7.0, 24))) + outcome = Estimation.estimate_models(nb_config(), counts_with_constant; + sample_metadata = meta_effect, + offset = offsets_effect, + taxa_ids = vcat(taxa_effect, ["t_constant"])) + @test outcome.status == :partial + @test outcome.diagnostics["n_failed"] == 1 + @test outcome.diagnostics["n_tested"] == 3 + + constant = outcome.results["t_constant"] + @test constant["status"] == "failed" + @test constant["pvalue"] === nothing + @test constant["padj"] === nothing + @test constant["estimate"] === nothing + @test occursin("constant", constant["note"]) + @test occursin("excluded", outcome.diagnostics["failure_note"]) + end + + @testset "logistic: presence/absence, with the design stated" begin + presence = zeros(3, 20) + presence[1, 1:18] .= 1.0; presence[1, 19:20] .= 0.0 # 18/20 vs 0/20 + presence[2, 1:10] .= 1.0; presence[2, 11:20] .= 1.0 # present everywhere + presence[3, 1:10] .= 0.0; presence[3, 11:20] .= 0.0 # absent everywhere + groups = vcat(fill("A", 10), fill("B", 10)) + + norm = AnalysisConfig.NormalizationConfig(method = "presence_absence") + config = AnalysisConfig.AnalysisConfig( + method = "logistic", formula = "~ group", outcome_column = "group", + metadata_columns = ["group"], normalization = norm, + advanced = AnalysisConfig.AdvancedConfig(min_samples_per_group = 2), + created_by = "test_estimation_logistic") + + outcome = Estimation.estimate_models(config, presence; + sample_metadata = OrderedDict{String,Any}("group" => groups), + taxa_ids = ["present_mostly_A", "present_all", "absent_all"]) + @test outcome.status == :partial # the two constant features are excluded + @test outcome.diagnostics["n_tested"] == 1 + + # The constant features say so and carry no number. + @test outcome.results["present_all"]["status"] == "failed" + @test outcome.results["absent_all"]["pvalue"] === nothing + + separator = outcome.results["present_mostly_A"] + @test separator["status"] in ("ok", "boundary") + @test separator["odds_ratio"] ≈ exp(separator["estimate"]) + # Separation is not silently forgiven: if R says it happened, the row says so. + if separator["status"] == "boundary" + @test occursin("separation", separator["note"]) + end + end + + @testset "Gaussian fit on CLR units records that its estimate is not a fold change" begin + clr = similar(counts_effect) + for j in 1:size(counts_effect, 2) + column = log.(counts_effect[:, j]) + clr[:, j] = column .- mean(column) + end + clr_norm = AnalysisConfig.NormalizationConfig(method = "clr") + config = AnalysisConfig.AnalysisConfig( + method = "clr_lm", formula = "~ group", metadata_columns = ["group"], + normalization = clr_norm, + advanced = AnalysisConfig.AdvancedConfig(min_samples_per_group = 2), + created_by = "test_estimation_clr_fit") + + outcome = Estimation.estimate_models(config, clr; + sample_metadata = meta_effect, + taxa_ids = taxa_effect) + @test outcome.status == :ok + @test outcome.results["t_up"]["estimate"] > 0 + @test outcome.results["t_down"]["estimate"] < 0 + @test occursin("clr_difference", outcome.results["t_up"]["estimate_scale"]) + @test !haskey(outcome.results["t_up"], "log2FoldChange") + @test outcome.results["t_up"]["statistic_name"] == "t" + end + + @testset "end to end: prepare_analysis_table -> run_analysis -> a number from the data" begin + config = nb_config() + sample_ids = ["s$i" for i in 1:24] + + (prepared, diagnostics, manifest, sids, tids) = Execution.prepare_analysis_table( + config, counts_effect; + sample_metadata = meta_effect, sample_ids = sample_ids, + taxa_ids = taxa_effect, drop_policy = "drop") + + # The offset the manifest carries is what the fit consumes; without it the GLM + # would compare library sizes. + @test !isnothing(manifest.offset) + @test length(manifest.offset) == 24 + + result = Execution.run_analysis( + Execution.RAdapter(method = "nb_glm"), config, prepared; + sample_metadata = meta_effect, diagnostics = diagnostics, manifest = manifest, + sample_ids = sids, taxa_ids = tids) + + @test result isa Execution.ExecutionResult + @test length(result.results) == 3 + @test result.provenance["estimation"]["status"] == "ok" + @test result.diagnostics.checks["estimation"]["n_tested"] == 3 + @test result.results["t_up"]["pvalue"] < 0.01 + @test result.results["t_up"]["padj"] >= result.results["t_up"]["pvalue"] - 1e-12 + @test result.results["t_up"]["status"] == "ok" + + # Nothing here is a function of the feature's name; the old mock was. + for (label, entry) in result.results + @test entry["pvalue"] isa Float64 + @test entry["status"] in ("ok", "boundary") + end + end + else + @info "SKIPPED (by name): live-R estimation tests — R with MASS is not reachable on this machine. CI installs both; a skip there is a provisioning defect, not a pass." + @test_skip Estimation.estimate_models(nb_config(), counts_effect; + sample_metadata = meta_effect, + offset = offsets_effect, taxa_ids = taxa_effect) + end + + # ------------------------------------------------------------------ + # The negative control that matters most: the placeholder statistics cannot come back. + # ------------------------------------------------------------------ + @testset "no placeholder statistics in the execution path" begin + execution_source = read(joinpath(@__DIR__, "..", "..", "src", "analysis", "Execution.jl"), String) + @test !occursin("hash(taxon_id)", execution_source) + @test !occursin("hash(taxa_id)", execution_source) + @test !occursin("mock BH", execution_source) + @test !occursin("Mock results", execution_source) + # And the estimator really is the one being called. + @test occursin("Estimation.estimate_models", execution_source) + + estimation_source = read(joinpath(@__DIR__, "..", "..", "src", "analysis", "estimation.jl"), String) + @test !occursin("hash(", estimation_source) + # A `rand(` here would mean a resampling step without a recorded seed; these fits have + # no random component at all, and the provenance says so. + @test !occursin("rand(", estimation_source) + end +end diff --git a/test/unit/test_execution.jl b/test/unit/test_execution.jl index fd8c48b4..daa328d0 100644 --- a/test/unit/test_execution.jl +++ b/test/unit/test_execution.jl @@ -365,7 +365,7 @@ @test any(prepared_rar .> 0) end - @testset "run_analysis — stub with manifest and diagnostics, hard-stop DANGER banner" begin + @testset "run_analysis — real estimator or an explicit unsuccessful state, hard-stop DANGER banner" begin norm = AnalysisConfig.NormalizationConfig(method="size_factors", epsilon=1e-6) adv = AnalysisConfig.AdvancedConfig(min_samples_per_group=2) config = AnalysisConfig.AnalysisConfig( @@ -388,18 +388,25 @@ drop_policy="drop" ) - # RAdapter stub + # No sample metadata means no design, and no design means no test. What must NOT + # happen is what used to happen here: three rows of statistics whose p-values came + # from `hash(taxon_id)`. An empty result carrying a stated reason is the honest + # answer, and these assertions are the guard that the placeholder stays gone. r_adapter = Execution.RAdapter(method="nb_glm") result_r = Execution.run_analysis(r_adapter, config, prepared; diagnostics=diagnostics, manifest=manifest, sample_ids=sids, taxa_ids=tids) @test result_r isa Execution.ExecutionResult @test result_r.config_id == config.id - @test length(result_r.results) == 3 + @test isempty(result_r.results) + @test result_r.provenance["estimation"]["status"] == "not_run" + @test occursin("descriptive summary", result_r.provenance["estimation"]["reason"]) + @test any(w -> occursin("Estimation not run", w), result_r.diagnostics.warnings) @test result_r.manifest.id == manifest.id - # JuliaAdapter stub + # Same for the Julia adapter: the adapter does not decide whether a design exists. jl_adapter = Execution.JuliaAdapter(method="nb_glm") result_jl = Execution.run_analysis(jl_adapter, config, prepared; diagnostics=diagnostics, manifest=manifest, sample_ids=sids, taxa_ids=tids) @test result_jl isa Execution.ExecutionResult + @test isempty(result_jl.results) # Adapter method mismatch should hard-stop with DANGER banner wrong_adapter = Execution.JuliaAdapter(method="clr_lm") # config is nb_glm