Fix statistical validity and data alignment issues across analyses#85
Merged
Merged
Conversation
Four verified correctness bugs in core statistical code: * linda.R / linda2.R: adaptive zero-handling assigned capitalized "Imputation"/"Pseudo-count" to `zero.handling`, but the dispatch test compares against lowercase "imputation". With the default `adaptive = TRUE`, the imputation branch could therefore never run and the code silently fell back to pseudo-count even when sequencing depth was confounded with the covariates -- exactly the case where imputation is intended. Assign lowercase values to match the test. * linda.R / linda2.R: `outlier.pct = NULL` (the documented "winsorization off" value) caused `1 - NULL` -> numeric(0) and broke winsorization. Guard the call with `!is.null(outlier.pct)` and correct the roxygen docs to state the real default (0.03) and NULL-skips semantics. * generate_taxa_change_test_pair.R: winsorization used apply(MARGIN = 2), capping values within each sample across unrelated taxa rather than per taxon across samples. This distorts each sample's composition and biases every taxon's change score. Switch to per-taxon (MARGIN = 1), matching winsor_feature_table() and the per-taxon pseudocount above. * generate_taxa_change_test_pair.R: subject filtering used colnames(na.omit(matrix)); na.omit drops rows (taxa), so it returned all subjects when any taxon was complete and none when every taxon had an NA -- never the intended per-subject filter. Select complete subject columns via colSums(is.na(...)) == 0. https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
* generate_{alpha,taxa,beta}_volatility_test_long.R: the omnibus group
F-test for >2-level groups used anova() (Type-I sequential) with the
group term entered FIRST, so the test ignored adj.vars even when the
user requested covariate adjustment. Reorder terms so group.var is
entered last; with no interaction in the volatility model this makes
SS(group | adj.vars) the covariate-adjusted test (equivalent to
Type-II/III for a main effect).
* mStat_calculate_adjusted_distance.R: stats::cmdscale() returns
principal coordinates already scaled by sqrt(eig), but the code
multiplied the residualized coordinates by sqrt(eig) again, weighting
axis i by eig_i instead of sqrt(eig_i) and over-emphasizing leading
axes. The buggy version also failed to reduce to the original distance
when covariates explain nothing. Take dist() of the residuals directly.
* mStat_calculate_alpha_diversity.R: Chao1 and ACE are defined only for
integer counts (singleton/doubleton based). Error clearly when given
non-integer/normalized data instead of silently returning meaningless
estimates.
* generate_alpha_volatility_test_long.R: guard the inline covariate
adjustment against NA covariates, which otherwise shorten the residual
vector and misalign/error the assignment back into alpha_df, matching
the existing check in mStat_calculate_adjusted_alpha_diversity().
https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
* mStat_convert_phyloseq_to_data_obj.R: the zero-sum feature filter used feature.tab[rowSums(...) > 0, ] without drop = FALSE, so a single surviving feature collapsed the matrix to a vector, dropping dim and dimnames and breaking every downstream rowname-based alignment (tree, feature.ann, analyses). The DESeq2 and SummarizedExperiment converters already use drop = FALSE; make phyloseq consistent. * mStat_convert_phyloseq_to_data_obj.R: building meta.dat via data.frame() applied check.names = TRUE, silently mangling metadata column names that contain spaces/special characters (e.g. "Body Site" -> "Body.Site") so downstream group.var/condition lookups failed. Add check.names = FALSE. * mStat_normalize_data.R: vegan::rrarefy() draws a random subsample but no seed was set, making "Rarefy" and "Rarefy-TSS" normalization non-reproducible across runs. Add a seed argument (default 123) and wrap the call in mStat_with_local_seed(), which restores the global RNG state so the caller's random stream is unaffected. https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
generate_beta_test_single() computes PERMANOVA via GUniFrac::PermanovaG2, which (through vegan::adonis2) pairs each distance matrix with the `data` rows by POSITION, not by sample name. When the distance object was freshly computed (dist.obj = NULL), the distances were labeled in feature.tab column order while the predictors were taken from data.obj$meta.dat in its own row order, with no realignment step -- only the provided-dist.obj branch reordered the distances to the metadata. A data.obj whose meta.dat rows were not already in feature.tab column order therefore produced PERMANOVA F/R2/p-values computed on shuffled group labels, with no error or warning. NA values in the group/adjustment variables were likewise not dropped consistently from both the distance matrix and the metadata. Align the metadata and every distance matrix to the same complete-case sample set, in a consistent order, for both the freshly-computed and provided-dist.obj paths before calling PermanovaG2. https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
… use * generate_alpha_per_time_test_long.R: the global group reference level was applied as the output label, but the within-timepoint lm fell back to the first factor level present in each subset. When a timepoint was missing the global reference group, the lm baseline differed from the label, mislabelling the contrast. Relevel the subset's group factor to the global reference_level when present (a no-op in the common case) so the fitted baseline matches the label. * generate_beta_test_single.R: warn when time.var has multiple levels but t.level is NULL. In that case PERMANOVA pools all timepoints and treats repeated measurements of the same subject as independent, violating exchangeability and inflating the Type-I error rate. https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
generate_taxa_test_single()'s perform_lm_analysis path reported the residual degrees of freedom as nrow(meta.dat) - length(all.vars(formula)) - 1. all.vars() counts each predictor once, but a k-level factor contributes k - 1 design columns (and interactions add more), so the reported df was overstated for any multi-level factor predictor. The fitted p-values come from lm's summary and were already correct; only this reported df column was wrong. Derive df from the actual model-matrix column count so it accounts for factor expansion. https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa
cafferychen777
had a problem deploying
to
github-pages
GitHub Actions
Failure
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
This PR addresses multiple statistical validity and data handling issues across the microbiome analysis pipeline, including PERMANOVA exchangeability violations, missing data handling, incorrect degrees of freedom calculations, and data alignment problems.
Key Changes
PERMANOVA and Beta Diversity Analysis
time.varhas multiple levels butt.levelis NULL, as pooling timepoints violates PERMANOVA's exchangeability assumption and inflates Type-I errorAlpha Diversity Analysis
model.matrix()to silently drop rows and misalign residualsTaxa Change Analysis
Adjusted Distance Calculation
cmdscale()already returns properly scaled coordinatesLinear Model Analysis
group.varlast, ensuring sequential (Type-I) ANOVA F-tests equal the covariate-adjusted omnibus testData Normalization and Conversion
seedparameter (default 123) to make rarefaction reproducible while preserving caller's RNG statedrop = FALSEto preserve matrix structure when subsetting single rows; addedcheck.names = FALSEto preserve metadata column namesStatistical Testing Functions
outlier.pctis not NULL, preventing unintended behavior whenis.winsor = TRUEbutoutlier.pct = NULL; updated documentation to clarify default behaviorNotable Implementation Details
https://claude.ai/code/session_01CDuBVvoJeszwbQGtXC7uFa