diff --git a/EXPLAINME.adoc b/EXPLAINME.adoc new file mode 100644 index 00000000..47e051f1 --- /dev/null +++ b/EXPLAINME.adoc @@ -0,0 +1,52 @@ += Occupancy and two-part abundance models: claim receipts +:toc: + +This page records what is implemented in the first slice of issue #18 and what +is intentionally *not* claimed. The remaining end-user occupancy workflow is +tracked in follow-up issue #205; code present today is not to be confused with +that outstanding work. + +== Implemented in the current PR + +[cols="1,2,2",options="header"] +|=== +|Claim |Implementation receipt |Boundary / evidence + +|Selectable ZINB differential model +|`src/analysis/differential.jl` (`fit_two_part_nb`, `differential_abundance`), config enum in `src/server/routes/config.jl` and `frontend/src/components/PipelineStages.tsx` +|Base R likelihood with zero and count components; the group effect is jointly tested by a 2-df likelihood-ratio test. Count offset is TSS/RLE. BH is applied by the shared result path. Backend Julia/R CI must pass before merge. + +|Selectable hurdle negative-binomial differential model +|Same fitter and config/UI path; the positive count likelihood is zero-truncated NB. +|All observed zeros are assigned to the occurrence process. It does not offer a lognormal or Poisson count component yet. + +|Model cautions and results +|`frontend/src/components/DifferentialPanel.tsx`, CSV result types/export, diagnostics in `differential_abundance` +|The UI surfaces a caution and exposes the zero-process coefficient. The joint LRT is not a model-selection test and does not establish that the chosen zero process is biologically real. + +|Repeated-visit occupancy estimator +|`MetaManifold.Differential.fit_occupancy`; tests in `test/unit/test_differential.jl` +|Accepts sites × visits × taxa Boolean detection histories; estimates group occupancy and common detection probability; refuses fewer than two visits. This is a lower-level API, not the standard analysis method or endpoint. + +|Scientific limitations +|`docs/occupancy-models.md`, `README.md` +|No single-visit occupancy claim. Closure and destructive microbiome sampling limitations are stated. Detection covariates are not implemented. +|=== + +== Not yet implemented + +The follow-up issue linked from #18 owns first-class occupancy configuration and +UI, occupancy/detection formula inputs and covariates, identifiability checks, +zero-formula/count-distribution options, model diagnostics/comparisons, reference +parity, scale benchmarks, schema/provenance completion, and a synchronized wiki +roadmap. Until that work lands, do not describe the occupancy estimator as a +ready-to-use workflow. + +== Validation record + +* Frontend typecheck, production build and frontend unit tests are run locally. +* Julia/R backend tests are defined in `test/unit/test_differential.jl` and run + by the repository CI workflow. Their successful CI result is required before + merge; this sandbox does not have Julia or R installed. +* CI also contains other workflows; each current failure must be inspected and + either fixed or shown to be an unrelated, pre-existing non-required lane. diff --git a/README.md b/README.md index 6e67d8c7..429d3778 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,7 @@ MetaManifold wraps standard amplicon sequencing workflows into a single configur ## Overview -MetaManifold consists of a Julia backend (pipeline engine + REST API) and a TypeScript/React frontend. The pipeline runs FastQC, MultiQC, cutadapt, DADA2, swarm, vsearch, and cd-hit-est under the hood, and places ASVs on reference trees with MAFFT, trimAl, IQ-TREE, RAxML and gappa; results are stored in per-run DuckDB databases and served to the frontend as interactive Plotly charts and filterable tables. Pipeline configuration is editable in the web UI at every cascade level (see [Configuration](#configuration)). +MetaManifold consists of a Julia backend (pipeline engine + REST API) and a TypeScript/React frontend. The pipeline runs FastQC, MultiQC, cutadapt, DADA2, swarm, vsearch, and cd-hit-est under the hood, and places ASVs on reference trees with MAFFT, trimAl, IQ-TREE, RAxML and gappa; results are stored in per-run DuckDB databases and served to the frontend as interactive Plotly charts and filterable tables. Pipeline configuration is editable in the web UI at every cascade level (see [Configuration](#configuration)). For analysis-method status and planned work, see [ROADMAP.md](ROADMAP.md); the full occupancy follow-up is tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205). **Pipeline stages** @@ -289,48 +289,38 @@ not exclude samples from this analysis. `analysis.differential.method` chooses the method: -- `nb_glm` (default): a negative-binomial GLM per taxon with a library-size - offset. The effect is log2 fold change; positive values indicate greater - abundance in the second condition. `offset` selects `tss` (default) or `rle` - scaling. `dispersion` selects per-taxon maximum likelihood with - `MASS::glm.nb` (default) or a cross-taxon shrinkage trend from R's - `glmGamPoi::glm_gp`, followed by fixed-theta NB Wald fits. +- `nb_glm` (default): a negative-binomial GLM per taxon using `MASS::glm.nb`, + with a library-size offset. The effect is log2 fold change; positive values + indicate greater abundance in the second condition. `offset` selects `tss` + (default) or `rle`; `dispersion` selects maximum-likelihood theta or a + cross-taxon `glmGamPoi` trend. +- `zinb`: a zero-inflated negative-binomial model with separate logistic + zero and count components, a count-component library-size offset, a joint + 2-df likelihood-ratio test, and BH correction. +- `hurdle_nb`: a logistic occurrence model plus zero-truncated NB counts for + positive observations. It uses the same joint test and offset convention. - `clr_welch`: Welch's t-test on centred log-ratios after zero replacement with `zCompositions::cmultRepl`. The effect is a CLR difference, not a fold change. - `replacement_method` selects `bayesian_multiplicative` (GBM, the default) or - `multiplicative_replacement` (CZM). `replacement_delta` controls the - replacement fraction (default 0.65), and optional positive `bayesian_alpha` - overrides GBM's estimated Dirichlet concentration; `offset` is not used. - [Zero handling, posterior semantics, warnings and references](docs/statistics/zero-handling.md) - explain why every replacement is biased and when a count, occupancy or - zero-inflated model is preferable. -- `ancom_bc`: ANCOM-BC sampling-fraction bias-corrected log-linear model - (Lin & Peddada 2020, *Nat. Commun.* 11:3514). Estimates within-group sample - offsets and between-group sampling-fraction bias via an E-M Gaussian mixture - on `log(counts + 1)`, propagating bias-estimation variance into two-sided Wald - tests. The effect is a bias-corrected log2 fold change; recommended when - asymmetric shifts or microbial load changes would bias standard CLR centering - (requires at least 3 tested taxa; assumes fewer than half are differential). -- `aldex2`: Dirichlet(`count + 0.5`) Monte-Carlo CLR model (Fernandes et al. - 2013, 2014). Draws `mc_samples` (default 128, seeded by the study's `seed`) - posterior instances per sample and reports standardized CLR effect sizes with - consistent-sign posterior predictive Welch's t-test p-values; recommended for - sparse or variable-depth cohorts where deterministic zero imputation is - undesirable. -- `songbird`: compositional baseline-category multinomial regression with - reference-frame invariant CLR-centred differentials summing to zero across - fitted taxa (Morton et al. 2019, *Nat. Commun.* 10:2719; Zhang et al. 2017). - Conditions on read totals without extra-multinomial biological overdispersion; - completely separated single-group taxa are marked `failed` rather than - diverging. - -For every method, `min_prevalence` (default 0) filters taxa by the fraction of -samples with nonzero counts, and Benjamini–Hochberg adjustment is mandatory -across the fitted taxa that produced a p-value. Failed or filtered taxa have no -substitute p-value; inspect their status and note in the results table. Invalid -input or an unavailable required R package produces an explicit error. See -`bench/compositional_differential/README.md` for method-selection guidance, -benchmarks, and literature citations. + `replacement_method` selects Bayesian-multiplicative (GBM, default) or + multiplicative (CZM) replacement; `replacement_delta` sets the replacement + fraction and optional `bayesian_alpha` overrides GBM's concentration. The + offset is not used. See [zero-handling notes](docs/statistics/zero-handling.md). + +ZINB and hurdle models can overfit sparse data; do not select ZINB solely +because counts are sparse. Compare against NB diagnostics. The methods report +a count-process effect and the zero-process coefficient; this is not proof of +biological structural absence. The public +`MetaManifold.Differential.fit_occupancy` lower-level estimator accepts +repeated-visit detection histories (sites × visits × taxa), refuses single-visit +histories and is not selectable in the standard analysis configuration or UI. +It currently has no detection covariates. Do not use it as a substitute for +the full workflow tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205). +See [occupancy model notes](docs/occupancy-models.md). +Across the differential methods, `min_prevalence` (default 0) filters taxa by the fraction of +samples with nonzero counts. Benjamini–Hochberg adjustment includes only taxa +that produced a p-value. Failed or filtered taxa have no substitute p-value; +inspect their status and note in the results table. Invalid input or an +unavailable required R package produces an explicit error. Implementation and reference tests: `src/analysis/differential.jl`, `src/server/routes/analysis.jl`, `bench/compositional_differential/benchmark.jl`, diff --git a/ROADMAP.md b/ROADMAP.md new file mode 100644 index 00000000..bc98c01c --- /dev/null +++ b/ROADMAP.md @@ -0,0 +1,37 @@ +# Analysis roadmap: occupancy and zero processes + +**Updated: 2026-10-22.** This is the repository-side status summary for the +occupancy work. The GitHub Wiki roadmap still needs a synchronized update; see +[follow-up issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205). + +## Available now + +- `nb_glm` and `clr_welch` remain the existing differential methods. +- `zinb` and `hurdle_nb` are available in the differential-abundance config and + UI. They report the count effect and zero/occurrence effect and test group + effects jointly. The UI surfaces model cautions; BH adjustment remains on. +- `MetaManifold.Differential.fit_occupancy` fits repeated-visit Boolean + detection histories (sites × visits × taxa), with group occupancy ψ and a + common detection probability p. This is a low-level estimator only. + +## Not available as a workflow yet + +Occupancy cannot yet be selected in AnalysisConfig or run through the regular +analysis endpoint. The estimator does not accept detection covariates/formulas, +primer or batch effects, or random effects; it refuses single-visit histories. +The ZINB/hurdle models do not yet have user-defined component formulas, model +selection diagnostics, or full scale benchmarks. These constraints and the +closure/destructive-sampling caveat are detailed in +[Occupancy model notes](docs/occupancy-models.md). + +## Remaining work + +Issue [#205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205) +tracks first-class repeated-detection input and configuration, detection and +zero formulas, identifiability validation, diagnostics/reference comparisons, +simulation recovery, benchmarks, provenance/schema completion, and updates to +the Wiki and user/developer roadmaps. Do not describe the lower-level fitter as +a complete occupancy workflow until those acceptance criteria pass. + +Issue [#18](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/18) +records the originating feature request and the first implementation slice. diff --git a/config/defaults/pipeline.yml b/config/defaults/pipeline.yml index 58a559bc..4d235d2e 100644 --- a/config/defaults/pipeline.yml +++ b/config/defaults/pipeline.yml @@ -346,26 +346,15 @@ analysis: # Differential abundance: compares two conditions on the raw read counts, with # mandatory Benjamini-Hochberg adjusted p-values across the fitted family. differential: - # The model fitted across the two conditions. - # - nb_glm: negative-binomial GLM on read counts with a library-size offset; - # the effect is a log2 fold change. - # - clr_welch: Welch's t-test on centred log-ratios over the taxa passing - # min_prevalence after zCompositions::cmultRepl zero replacement; the - # effect is a CLR difference, not a fold change. - # - ancom_bc: bias-corrected log-linear model (Lin & Peddada 2020) that - # estimates sample-specific and between-group sampling-fraction shifts via - # an E-M Gaussian mixture; the effect is a bias-corrected log2 fold change. - # - aldex2: Dirichlet(count + 0.5) Monte-Carlo CLR instances with consistent-sign - # Welch's t-test and Fernandes et al. (2014) standardized CLR effect size. - # - songbird: compositional baseline-category multinomial regression with - # reference-frame invariant CLR-centred differentials (Morton et al. 2019; - # Zhang et al. 2017); the effect is a CLR log-ratio differential, not an - # absolute fold change. + # Per-taxon model: nb_glm (negative-binomial GLM), clr_welch (Welch test on + # zero-replaced CLR values), zinb (zero-inflated NB), or hurdle_nb (logistic + # occurrence plus zero-truncated NB). Two-part models test both group effects + # jointly; their count-process fold change is reported separately. method: nb_glm # How library size enters each model, as a log offset: tss (total reads per # sample) or rle (median-of-ratios over the taxa present in every sample). # rle is refused when no taxon is present in every sample. - # Used by nb_glm only; compositional methods ignore it. + # Used by count models (nb_glm, zinb, hurdle_nb); clr_welch ignores it. offset: tss # nb_glm only: where each taxon's negative-binomial theta comes from. ml # (the default) takes the maximum-likelihood estimate MASS::glm.nb fits per diff --git a/docs/occupancy-models.md b/docs/occupancy-models.md new file mode 100644 index 00000000..46a28661 --- /dev/null +++ b/docs/occupancy-models.md @@ -0,0 +1,87 @@ +# Occupancy and two-part abundance models + +**Status: partial implementation.** `zinb` and `hurdle_nb` are selectable in +differential abundance. `fit_occupancy` is a repeated-visit lower-level API, +not an AnalysisConfig method or the standard endpoint. The full workflow is +tracked in [issue #205](https://github.com/hyperpolymath/MetaManifold-WebUI/issues/205), +a follow-up to #18. + +The differential-abundance configuration supports `zinb` and `hurdle_nb` in +addition to the existing `nb_glm` and `clr_welch` models. The two new methods +are fixed-effect, per-taxon models with a two-part negative-binomial likelihood. +They use the configured `tss` or `rle` library-size offset in the count +component and apply Benjamini–Hochberg correction to the joint tests. + +## Zero-inflated negative binomial (ZINB) + +For taxon count `Y`, the model has a logistic structural-zero probability +`π`, and an NB count distribution with mean `μ` and dispersion `θ`: + +- `P(Y = 0) = π + (1 - π) NB(0; μ, θ)` +- `P(Y = y > 0) = (1 - π) NB(y; μ, θ)` + +Both the zero component and log-mean count component include a group effect; +the latter includes the log library-size offset. The reported count effect is +the log fold change (and its log2 conversion). The result also reports the +zero-component group log-odds coefficient. The p-value is a joint 2-df +likelihood-ratio test against a model without either group effect. The fits use +base R's `stats::optim` and `stats::dnbinom`, so no extra R modeling package is +required. + +## Hurdle negative binomial + +A logistic component models whether a count is zero. Conditional on a positive +count, the abundance follows an NB distribution truncated at zero. All zeros +are assigned to the occurrence component; unlike ZINB, the count component +cannot produce zero. The same joint test, count-process effect, offset, and +multiple-testing policy apply. + +These are not interchangeable explanations of sparse data. NB overdispersion +can itself produce many zeros; selecting ZINB just because the matrix is sparse +can overfit and make the structural-zero component weakly identified. Compare +models using diagnostics and report the chosen model. Results do not currently +include a formal Vuong comparison or DHARMa residual diagnostic. + +## Repeated-visit occupancy estimator + +The lower-level Julia entry point +`MetaManifold.Differential.fit_occupancy(history, groups, taxa; levels)` accepts +a Boolean array with dimensions **sites × repeat visits × taxa**. It fits an +occupancy probability `ψ` with an intercept and binary group effect and a +detection probability `p` shared across visits. `psi_reference` is the fitted +occupancy probability for `levels[1]` (the reference); `psi_contrast` is for +`levels[2]` (the contrast). A per-taxon likelihood-ratio test assesses the +group effect in occupancy; p-values are BH-adjusted. The returned table also +includes the estimated `p`. + +This is intentionally not selected by `analysis.differential.method` and is +not called by the single-observation differential endpoint: that endpoint has +no repeat-visit design or detection-history input. The estimator refuses fewer +than two visits rather than treating technical or biological samples as repeat +visits. The initial implementation does not yet accept detection covariates +such as library size, batch, or primer, and does not implement arbitrary +occupancy/detection formulas or random effects. A deployment that needs these +should collect a defensible repeated-detection design and supply/validate +covariates before interpreting `ψ` and `p`. + +Occupancy assumes closure over repeat visits: a site's true occupancy does not +change during the detection period. Microbiome sampling is often destructive, +and the meaning of a repeated visit may therefore be questionable. Single-visit +occupancy is not empirically separated into ψ and p without additional strong +assumptions. Interpret results cautiously and document the sampling design. + +## Provenance and limitations + +Pipeline configuration retains the selected two-part method, offset, +prevalence threshold, and dispersion option. The current analysis assumes a +single binary group contrast and has no user-supplied zero/detection formula. +The two-part tests use a 2-df LRT, not an individual-component significance +test; their reported count fold change is only one part of the joint hypothesis. + +References: + +- MacKenzie, D. I. et al. (2002). Estimating site occupancy rates when detection + probabilities are less than one. *Ecology* 83(8), 2248–2255. +- Martin, T. G. et al. (2005). Zero tolerance ecology: improving ecological + inference by modelling the source of zero observations. *Ecology Letters* 8, + 1235–1246. diff --git a/docs/statistics/zero-handling.md b/docs/statistics/zero-handling.md index 66095019..117d7105 100644 --- a/docs/statistics/zero-handling.md +++ b/docs/statistics/zero-handling.md @@ -1,134 +1,60 @@ -# Zero handling and dispersion for differential abundance - -## Why zeros need an explicit decision - -A centred log-ratio (CLR) is `log(x_ij) - mean_j(log(x_ij))`, so it is undefined -when any part of the composition is zero. A sequencing zero can mean biological -absence or a taxon that was present below the run's detection limit; the count -table alone does not distinguish those explanations. - -**Every zero replacement is biased. A replaced zero is not a measurement.** -Multiplicative methods avoid some of the distortion caused by adding a -pseudocount, but they cannot recover an unobserved abundance. Do not try several -parameters and report only the run with the desired significance. The selected -method, delta and (when Bayesian) alpha are stored in the resolved run -configuration, the analysis response and the CSV settings line. - -A negative-binomial model does not have this problem: its likelihood is defined -at zero, so `nb_glm` analyses the observed integer counts without a pseudocount -or replacement. - -## Replacement policies for `clr_welch` - -The Analysis configuration expander exposes -`analysis.differential.replacement_method`. Both policies call the pinned R -reference implementation, `zCompositions::cmultRepl`, with `z.warning = 1` and -`z.delete = FALSE`. Prevalence filtering is deliberately done by MetaManifold -before replacement; cmultRepl is never allowed to drop a sparse row or taxon and -silently change the tested family. - -| Setting | Reference call | Behaviour | -| --- | --- | --- | -| `bayesian_multiplicative` (default) | `cmultRepl(..., method = "GBM", output = "prop")` | Geometric Bayesian-multiplicative posterior replacement with a leave-one-out prior mean and an estimated Dirichlet concentration. | -| `bayesian_multiplicative` with `bayesian_alpha > 0` | `cmultRepl(..., method = "user", t = ..., s = alpha, output = "prop")` | The same leave-one-out prior mean with an explicitly selected positive Dirichlet concentration. | -| `multiplicative_replacement` | `cmultRepl(..., method = "CZM", output = "prop")` | Count-zero multiplicative replacement. Before adjustment, a zero receives `delta * 0.5 / sample_total`. | - -For both policies, cmultRepl closes each sample to a composition and scales its -observed parts by one common factor. Consequently, each output row sums to one -and ratios among observed parts are preserved. This is the important difference -from adding the same pseudocount to every cell, which changes observed-part -ratios. It does **not** make the imputed values unbiased. - -A taxon seen in fewer than two samples is refused for Bayesian replacement -because its leave-one-out prior is undefined. CZM needs at least one observation -of every included taxon. An all-zero sample, a never-observed CZM taxon, fewer -than two taxa, or an invalid resulting composition is refused explicitly rather -than repaired with a guessed value. - -### Delta - -`analysis.differential.replacement_delta` is cmultRepl's `frac` and must be -strictly between 0 and 1; the default is 0.65. - -- For CZM it multiplies the `0.5 / sample_total` threshold. -- With `adjust = TRUE` (the reference default used here), either policy caps an - imputed proportion that exceeds its taxon's smallest observed proportion at - `delta` times that proportion. -- Values below 0.01 produce a numerical caution because very small replacements - can underflow. Values at least 0.9 produce a caution because imputed values - can receive substantial mass relative to observed parts. These are numerical - heuristics, not biologically optimal cutoffs. - -Delta is recorded even when a particular table has no zeros, so the configuration -that produced a result remains reproducible. - -### Bayesian posterior, without pretending to sample it - -For a Bayesian replacement with at least one zero, `replace_zeros` also returns -`posterior.prior_mean` and `posterior.concentration`. For sample `i`, these fully -parameterise the pre-adjustment posterior - -```text -Dirichlet(observed_counts[i, :] + concentration[i] * prior_mean[i, :]) -``` - -The analysis diagnostics record the concentration for every sample and this -parameterisation. cmultRepl's subsequent `frac` cap and multiplicative closure -are deterministic adjustments, not draws from that Dirichlet distribution. The -Web UI therefore describes the result as an **adjusted posterior estimate**, not -as a posterior sample or an uncertainty interval. - -Alpha is optional (`null` by default), must be finite and strictly positive when -set, and is ignored by CZM. A larger alpha gives the prior more weight relative -to the observed counts. - -## glmGamPoi dispersion for `nb_glm` - -`analysis.differential.dispersion: glmGamPoi` calls the pinned R -`glmGamPoi::glm_gp` once across all taxa passing the prevalence filter and takes -`overdispersion_shrinkage_list.dispersion_trend`. Each taxon is then fitted with -`MASS::negative.binomial(theta = 1 / dispersion)` and a Wald test with fixed -scale 1. The default `ml` path remains unchanged and uses `MASS::glm.nb`. - -This is plug-in negative-binomial Wald inference using glmGamPoi's shrunken -dispersion trend. It is **not** glmGamPoi's quasi-likelihood test, and it does -not propagate uncertainty in the estimated trend. Selecting it for `clr_welch` -is refused. Missing R packages produce an explicit service-unavailable error; -there is no silent fallback to another estimator. - -## When not to replace zeros - -- Use `nb_glm` when a count model answers the scientific question; raw zeros are - valid observations under that model. -- Raise `min_prevalence` when rare taxa would make imputed values dominate the - composition. The result reports median and maximum imputed mass per sample. -- Consider an occupancy/observation model when detection itself is the object of - inference, or a zero-inflated or hurdle model when structural and sampling - zeros must be separated. Those models are not implemented by `clr_welch`, and - multiplicative replacement is not a substitute for them. +# Zero handling in differential abundance + +**Status:** implemented for `clr_welch`; updated 2026-10-22. + +Zeros are observed values, not missing data. The `clr_welch` method requires a +strictly positive composition before taking logarithms, so it replaces zeros +with `zCompositions::cmultRepl` after filtering taxa by +`analysis.differential.min_prevalence`. Replaced values are imputed inputs to +the log-ratio calculation; they are not measurements and do not establish +structural absence. + +## Replacement options + +- `bayesian_multiplicative` (default) uses GBM, the geometric + Bayesian-multiplicative replacement. It uses a leave-one-out taxon profile as +a prior mean. `bayesian_alpha: ~` uses the estimated Dirichlet concentration; +a finite positive value uses the same prior mean with the requested +concentration. Posterior parameters are recorded in result diagnostics. The +implementation uses posterior estimates, not posterior draws. +- `multiplicative_replacement` uses CZM, the count zero-multiplicative + replacement. Its starting threshold depends on `replacement_delta` and the +sample total. Both methods preserve ratios among observed parts before +multiplicative closure, and both cap high imputations relative to the taxon's +smallest observed proportion. + +`replacement_delta` must be strictly between 0 and 1. Values below 0.01 or at +least 0.9 trigger a caution. It is not a tuning parameter for significance. +Bayesian replacement requires each retained taxon to occur in at least two +samples; CZM requires at least one. A composition needs at least two retained +taxa. Input that cannot produce a valid replacement is refused rather than +silently converted to a pseudocount result. + +## Interpretation and alternatives + +The `clr_welch` estimate is a contrast-minus-reference difference in centred +log-ratio units, not a fold change. Replacement choice and parameters are +returned in `config`; zero counts, imputed-mass summaries, warnings, and (for +GBM) concentration information are returned in `diagnostics.zero_replacement`. + +Every zero-replacement method introduces assumptions and can bias results, +especially when zeros are common or sample sizes are small. Do not choose a +replacement or `delta` to obtain a preferred p-value. Consider the negative +binomial methods on raw counts (`nb_glm`, `zinb`, `hurdle_nb`) when appropriate. +Repeated-detection occupancy is a separate, incomplete workflow; see +[occupancy model notes](../occupancy-models.md) and +[roadmap](../../ROADMAP.md). ## References - Martín-Fernández, J. A., Barceló-Vidal, C. & Pawlowsky-Glahn, V. (2003). - *Dealing with zeros and missing values in compositional data sets using - nonparametric imputation*. Mathematical Geology 35, 253–278. - -- Martín-Fernández, J. A., Hron, K., Templ, M., Filzmoser, P. & - Palarea-Albaladejo, J. (2015). *Bayesian-multiplicative treatment of count - zeros in compositional data sets*. Statistical Modelling 15(2), 134–158. - -- Palarea-Albaladejo, J. & Martín-Fernández, J. A. (2015). *zCompositions — R - package for multivariate imputation of left-censored data under a - compositional approach*. Chemometrics and Intelligent Laboratory Systems 143, - 85–96. -- Ahlmann-Eltze, C. & Huber, W. (2020). *glmGamPoi: fitting Gamma-Poisson - generalized linear models on single cell count data*. Bioinformatics 36(24), - 5701–5702. - -Implementation and direct-reference tests live in -`src/analysis/differential.jl` and `test/unit/test_differential.jl`. The -opt-in 100/1,000/10,000-taxon runtime and allocation benchmark is documented in -`bench/advanced_zero_handling/README.md`. + Dealing with zeros and missing values in compositional data sets using + nonparametric imputation. *Mathematical Geology*, 35, 253–278. +- Martín-Fernández, J. A., Hron, K., Templ, M., Filzmoser, P. & Palarea-Albaladejo, + J. (2015). Bayesian-multiplicative treatment of count zeros in compositional + data sets. *Statistical Modelling*, 15(2), 134–158. +- The implementation uses the documented `zCompositions::cmultRepl` methods; + consult that package's version-matched documentation for exact details. diff --git a/frontend/src/api/types.ts b/frontend/src/api/types.ts index ce8f9951..e668807d 100644 --- a/frontend/src/api/types.ts +++ b/frontend/src/api/types.ts @@ -418,10 +418,12 @@ export interface DifferentialRow { padj: number | null /** nb_glm only. */ dispersion_theta: number | null + /** Group coefficient in the logistic zero/occurrence component (two-part NB only). */ + zero_process_log_odds?: number | null prevalence: number } -export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'ancom_bc' | 'aldex2' | 'songbird' +export type DifferentialMethod = 'nb_glm' | 'clr_welch' | 'zinb' | 'hurdle_nb' /** Which row field holds the effect size, and what to call it. */ export interface DifferentialEffect { @@ -437,12 +439,12 @@ export interface DifferentialResult { groups: { reference: string; contrast: string } n_samples: Record effect: DifferentialEffect - /** Empty for compositional methods, which take no library-size offset. */ + /** Empty for clr_welch; count models report their size factors. */ size_factors: { sample: string; group: string; factor: number }[] config: { method: DifferentialMethod min_prevalence: number - /** nb_glm only. */ + /** Count-process offset for nb_glm, zinb and hurdle_nb. */ offset?: 'tss' | 'rle' /** nb_glm only: where each taxon's theta comes from. */ dispersion?: 'ml' | 'glmGamPoi' @@ -458,6 +460,8 @@ export interface DifferentialResult { } diagnostics: { n_taxa: number; n_tested: number; n_failed: number; n_boundary: number; n_filtered: number + /** Safety note for two-part models. */ + model_warning?: string /** clr_welch only. */ zero_replacement?: { method: string diff --git a/frontend/src/components/DifferentialPanel.tsx b/frontend/src/components/DifferentialPanel.tsx index 6648eae1..f7dca8cf 100644 --- a/frontend/src/components/DifferentialPanel.tsx +++ b/frontend/src/components/DifferentialPanel.tsx @@ -292,6 +292,13 @@ export function DifferentialPanel({ study, runs, option, aggregate }: { + {result.diagnostics.model_warning && ( +
+ Model caution: {result.diagnostics.model_warning} +
+ )} +