From 5df9d2b7a859743bef7ab59ea1ce4391816a4dd6 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 01:57:36 +0000 Subject: [PATCH 1/6] Fix statistical bugs in LinDA zero-handling and paired change tests 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 --- R/generate_taxa_change_test_pair.R | 21 +++++++++++++++++---- R/linda.R | 8 ++++---- R/linda2.R | 8 ++++---- 3 files changed, 25 insertions(+), 12 deletions(-) diff --git a/R/generate_taxa_change_test_pair.R b/R/generate_taxa_change_test_pair.R index 9ee79cb..74b6aac 100644 --- a/R/generate_taxa_change_test_pair.R +++ b/R/generate_taxa_change_test_pair.R @@ -198,13 +198,19 @@ generate_taxa_change_test_pair <- message("Zero-handling: Per-taxon half-minimum pseudocount calculated across ALL samples and time points combined. ", "This ensures unbiased change calculations by using the same pseudocount at both time points.") - # Apply winsorization to limit extreme values - otu_tax_agg_filter[, -1] <- apply(otu_tax_agg_filter[, -1], 2, function(x) { + # Apply winsorization to limit extreme values. + # Winsorize each taxon (row) across samples so that outlier samples for a + # given taxon are capped; winsorizing across taxa within a sample (MARGIN = 2) + # would mix unrelated taxa and distort each sample's composition. This matches + # the per-taxon convention used by winsor_feature_table() and the per-taxon + # pseudocount imputation applied above. t() keeps the assignment orientation + # consistent with otu_tax_agg_filter[, -1] (taxa in rows, samples in columns). + otu_tax_agg_filter[, -1] <- t(apply(otu_tax_agg_filter[, -1], 1, function(x) { qt <- quantile(x, probs = c((1 - winsor.qt) / 2, 1 - (1 - winsor.qt) / 2)) x[x < qt[1]] <- qt[1] x[x > qt[2]] <- qt[2] return(x) - }) + })) } # "other": no preprocessing — respect user's pre-processed data @@ -241,8 +247,15 @@ generate_taxa_change_test_pair <- as.matrix() # Prepare metadata for analysis + # Keep only subjects (columns) that have a complete set of change scores. + # value_diff_matrix is taxa (rows) x subjects (columns); selecting complete + # columns retains subjects with no missing change values. (na.omit() on a + # matrix drops rows, so the previous colnames(na.omit(...)) returned every + # subject whenever any taxon row was complete, and none when every taxon had + # an NA -- never the intended per-subject filter.) + complete_subjects <- colSums(is.na(value_diff_matrix)) == 0 aligned_subjects <- mStat_align_subject_metadata_to_matrix( - value_matrix = value_diff_matrix[, colnames(na.omit(value_diff_matrix)), drop = FALSE], + value_matrix = value_diff_matrix[, complete_subjects, drop = FALSE], meta_tab = meta_tab, subject.var = subject.var, keep_vars = c(group.var, adj.vars) diff --git a/R/linda.R b/R/linda.R index 6947527..52bc017 100644 --- a/R/linda.R +++ b/R/linda.R @@ -31,7 +31,7 @@ #' @param mean.abund.filter A real value; taxa with mean abundance less than mean.abund.filter are excluded. Default is 0 (no taxa will be excluded). #' @param max.abund.filter A real value; taxa with max abundance less than max.abund.filter are excluded. Default is 0 (no taxa will be excluded). #' @param is.winsor Boolean. If TRUE (default), the Winsorization process will be conducted for the feature table. -#' @param outlier.pct A real value between 0 and 1; Winsorization cutoff (percentile) for the feature table, e.g., 0.03. Default is NULL. If NULL, Winsorization process will not be conducted. +#' @param outlier.pct A real value between 0 and 1; Winsorization cutoff (percentile) for the feature table, e.g., 0.03. Default is 0.03. If NULL, the Winsorization process will be skipped even when is.winsor is TRUE. #' @param adaptive Boolean. Default is TRUE. If TRUE, the parameter imputation will be treated as FALSE no matter what it is actually set to be. Then the significant correlations between the sequencing depth and explanatory variables will be tested via the linear regression between the log of the sequencing depths and formula. If any p-value is smaller than or equal to corr.cut, the imputation approach will be used; otherwise, the pseudo-count approach will be used. #' @param zero.handling Character. Specifies the method to handle zeros in the feature table. Options are "pseudo-count" or "imputation" (default is "pseudo-count"). If "imputation", zeros in the feature table will be imputed using the formula in the referenced paper. If "pseudo-count", a small constant (pseudo.cnt) will be added to each value in the feature table. #' @param pseudo.cnt A positive real value. Default is 0.5. If zero.handling is set to "pseudo-count", this constant will be added to each value in the feature table. @@ -228,7 +228,7 @@ linda <- function(feature.dat, meta.dat, phyloseq.obj = NULL, formula, feature.d Z[, ind] <- scale(Z[, ind]) # Apply Winsorization to handle outliers if specified - if (is.winsor) { + if (is.winsor && !is.null(outlier.pct)) { Y <- winsor_feature_table(Y, 1 - outlier.pct, feature.dat.type) } @@ -274,12 +274,12 @@ linda <- function(feature.dat, meta.dat, phyloseq.obj = NULL, formula, feature.d if (verbose) { message("Imputation approach is used.") } - zero.handling <- "Imputation" + zero.handling <- "imputation" } else { if (verbose) { message("Pseudo-count approach is used.") } - zero.handling <- "Pseudo-count" + zero.handling <- "pseudo-count" } } if (zero.handling == "imputation") { diff --git a/R/linda2.R b/R/linda2.R index c8aa309..885b520 100644 --- a/R/linda2.R +++ b/R/linda2.R @@ -268,7 +268,7 @@ get_tree_smoothing_info <- function(phy.tree, tax.names, #' @param mean.abund.filter A real value; taxa with mean abundance less than mean.abund.filter are excluded. Default is 0 (no taxa will be excluded). #' @param max.abund.filter A real value; taxa with max abundance less than max.abund.filter are excluded. Default is 0 (no taxa will be excluded). #' @param is.winsor Boolean. If TRUE (default), the Winsorization process will be conducted for the feature table. -#' @param outlier.pct A real value between 0 and 1; Winsorization cutoff (percentile) for the feature table, e.g., 0.03. Default is NULL. If NULL, Winsorization process will not be conducted. +#' @param outlier.pct A real value between 0 and 1; Winsorization cutoff (percentile) for the feature table, e.g., 0.03. Default is 0.03. If NULL, the Winsorization process will be skipped even when is.winsor is TRUE. #' @param adaptive Boolean. Default is TRUE. If TRUE, the parameter imputation will be treated as FALSE no matter what it is actually set to be. Then the significant correlations between the sequencing depth and explanatory variables will be tested via the linear regression between the log of the sequencing depths and formula. If any p-value is smaller than or equal to corr.cut, the imputation approach will be used; otherwise, the pseudo-count approach will be used. #' @param zero.handling Character. Specifies the method to handle zeros in the feature table. Options are "pseudo-count" or "imputation" (default is "pseudo-count"). If "imputation", zeros in the feature table will be imputed using the formula in the referenced paper. If "pseudo-count", a small constant (pseudo.cnt) will be added to each value in the feature table. #' @param pseudo.cnt A positive real value. Default is 0.5. If zero.handling is set to "pseudo-count", this constant will be added to each value in the feature table. @@ -656,7 +656,7 @@ linda2 <- function(feature.dat, meta.dat, phyloseq.obj = NULL, formula, feature. Z[, ind] <- scale(Z[, ind]) # Apply Winsorization to handle outliers if specified - if (is.winsor) { + if (is.winsor && !is.null(outlier.pct)) { Y <- winsor_feature_table(Y, 1 - outlier.pct, feature.dat.type) } @@ -712,12 +712,12 @@ linda2 <- function(feature.dat, meta.dat, phyloseq.obj = NULL, formula, feature. if (verbose) { message("Imputation approach is used.") } - zero.handling <- "Imputation" + zero.handling <- "imputation" } else { if (verbose) { message("Pseudo-count approach is used.") } - zero.handling <- "Pseudo-count" + zero.handling <- "pseudo-count" } } if (zero.handling == "imputation") { From 8aac46173f9e9a61fbefc24c911b4ef039140e7b Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 03:07:38 +0000 Subject: [PATCH 2/6] Fix covariate adjustment and richness-estimator rigor issues * 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 --- R/generate_alpha_volatility_test_long.R | 18 +++++++++++++++++- R/generate_beta_volatility_test_long.R | 4 +++- R/generate_taxa_volatility_test_long.R | 4 +++- R/mStat_calculate_adjusted_distance.R | 11 ++++++++--- R/mStat_calculate_alpha_diversity.R | 16 ++++++++++++++++ 5 files changed, 47 insertions(+), 6 deletions(-) diff --git a/R/generate_alpha_volatility_test_long.R b/R/generate_alpha_volatility_test_long.R index 60993fb..abc3f65 100644 --- a/R/generate_alpha_volatility_test_long.R +++ b/R/generate_alpha_volatility_test_long.R @@ -112,6 +112,19 @@ generate_alpha_volatility_test_long <- function(data.obj, dplyr::select(all_of(adj.vars)) %>% dplyr::mutate(dplyr::across(where(is.character) & !is.factor, factor)) + # Guard against missing covariate values. model.matrix() / lm() drop NA rows + # by default, which would make the residual vector shorter than alpha_df and + # misalign (or error on) the assignment back into alpha_df[[index]], silently + # corrupting the per-subject volatility computed below. This mirrors the + # explicit check in mStat_calculate_adjusted_alpha_diversity(). + if (anyNA(data_subset)) { + stop( + "Adjustment variables contain missing values. Please remove or impute ", + "missing covariates before adjusting alpha diversity.", + call. = FALSE + ) + } + factor_terms <- names(data_subset)[vapply(data_subset, is.factor, logical(1))] contrasts_arg <- if (length(factor_terms) > 0) { lapply(data_subset[factor_terms], stats::contrasts, contrasts = FALSE) @@ -170,7 +183,10 @@ generate_alpha_volatility_test_long <- function(data.obj, valid_terms <- mStat_resolve_variable_terms( data = test_df, - terms = c(group.var, adj.vars) + # group.var is placed last so that the sequential (Type-I) anova() F-test + # for the group term equals SS(group | adj.vars), i.e. the covariate-adjusted + # omnibus test. There is no interaction here, so this matches Type-II/III. + terms = c(adj.vars, group.var) ) model_formula <- mStat_build_formula( diff --git a/R/generate_beta_volatility_test_long.R b/R/generate_beta_volatility_test_long.R index 5c7087d..2970cda 100644 --- a/R/generate_beta_volatility_test_long.R +++ b/R/generate_beta_volatility_test_long.R @@ -170,7 +170,9 @@ generate_beta_volatility_test_long <- # Test the association between volatility and group variable using linear regression valid_terms <- mStat_resolve_variable_terms( data = test_df, - terms = c(group.var, adj.vars) + # group.var last so the sequential (Type-I) anova() F-test for the group + # term equals SS(group | adj.vars), the covariate-adjusted omnibus test. + terms = c(adj.vars, group.var) ) test_result <- lm( mStat_build_formula(response = "volatility", terms = valid_terms), diff --git a/R/generate_taxa_volatility_test_long.R b/R/generate_taxa_volatility_test_long.R index 320663f..442aecd 100644 --- a/R/generate_taxa_volatility_test_long.R +++ b/R/generate_taxa_volatility_test_long.R @@ -193,7 +193,9 @@ generate_taxa_volatility_test_long <- function(data.obj, valid_terms <- mStat_resolve_variable_terms( data = test_df, - terms = c(group.var, adj.vars) + # group.var last so the sequential (Type-I) anova() F-test for the group + # term equals SS(group | adj.vars), the covariate-adjusted omnibus test. + terms = c(adj.vars, group.var) ) model_formula <- mStat_build_formula( diff --git a/R/mStat_calculate_adjusted_distance.R b/R/mStat_calculate_adjusted_distance.R index fea763f..e6eb9c3 100644 --- a/R/mStat_calculate_adjusted_distance.R +++ b/R/mStat_calculate_adjusted_distance.R @@ -108,7 +108,6 @@ mStat_calculate_adjusted_distance <- function (data.obj, } mds_coordinates <- mds_coordinates[, positive_axes, drop = FALSE] - axis_weights <- sqrt(mds_result$eig[seq_len(ncol(mds_result$points))][positive_axes]) res_matrix <- matrix(NA_real_, nrow = nrow(mds_coordinates), ncol = ncol(mds_coordinates)) rownames(res_matrix) <- rownames(mds_coordinates) @@ -124,8 +123,14 @@ mStat_calculate_adjusted_distance <- function (data.obj, res_matrix[, i] <- stats::residuals(model_fit) } - weighted_coordinates <- sweep(res_matrix, 2, axis_weights, `*`) - D.adj <- stats::dist(weighted_coordinates) + # stats::cmdscale() returns principal coordinates already scaled by sqrt(eig) + # (i.e. dist(points) reconstructs the original distance), so the residualized + # coordinates inherit that scaling. Re-multiplying by sqrt(eig) here would + # weight axis i by eig_i instead of sqrt(eig_i), squaring the eigenvalue + # weighting and over-emphasizing the leading axes. Taking dist() of the + # residuals directly correctly reduces to the original distance when the + # adjustment covariates explain nothing. + D.adj <- stats::dist(res_matrix) mStat_attach_dist_metadata(D.adj, data.obj$meta.dat) }) diff --git a/R/mStat_calculate_alpha_diversity.R b/R/mStat_calculate_alpha_diversity.R index 50e06fc..d3e90d2 100644 --- a/R/mStat_calculate_alpha_diversity.R +++ b/R/mStat_calculate_alpha_diversity.R @@ -54,6 +54,22 @@ mStat_calculate_alpha_diversity <- function(x, alpha.name, tree = NULL) { warning("It appears the data may not have been rarefied. Please verify.") } + # Chao1 and ACE are richness estimators defined only for integer abundance + # counts: they are functions of the numbers of singletons/doubletons (and rare + # species). Applying them to normalized/relative-abundance data yields + # statistically meaningless estimates, so require integer counts here. + if (any(c("chao1", "ace") %in% alpha.name)) { + finite_x <- x[is.finite(x)] + if (length(finite_x) > 0 && !isTRUE(all.equal(finite_x, round(finite_x)))) { + stop( + "Chao1 and ACE require integer count data, but non-integer values were ", + "detected. These estimators are based on singleton/doubleton counts and ", + "are not defined for normalized or relative-abundance data. Provide raw ", + "counts, or remove 'chao1'/'ace' from alpha.name." + ) + } + } + x_transpose <- t(x) alpha.obj <- lapply(alpha.name, function(index) { From 8edd4070e0ffa9f12072127294f6d46d3b52abd9 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 04:09:01 +0000 Subject: [PATCH 3/6] Fix phyloseq conversion and make rarefaction reproducible * 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 --- R/mStat_convert_phyloseq_to_data_obj.R | 4 ++-- R/mStat_normalize_data.R | 18 +++++++++++++++--- man/mStat_normalize_data.Rd | 7 ++++++- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/R/mStat_convert_phyloseq_to_data_obj.R b/R/mStat_convert_phyloseq_to_data_obj.R index 96b304d..3f3eb7e 100644 --- a/R/mStat_convert_phyloseq_to_data_obj.R +++ b/R/mStat_convert_phyloseq_to_data_obj.R @@ -47,11 +47,11 @@ mStat_convert_phyloseq_to_data_obj <- function (phylo.obj) { data.obj$feature.tab <- otu_matrix - data.obj$feature.tab <- data.obj$feature.tab[rowSums(data.obj$feature.tab) > 0, ] + data.obj$feature.tab <- data.obj$feature.tab[rowSums(data.obj$feature.tab) > 0, , drop = FALSE] } if (!is.null(phylo.obj@sam_data)) { - data.obj$meta.dat <- data.frame(phylo.obj@sam_data, stringsAsFactors = FALSE) + data.obj$meta.dat <- data.frame(phylo.obj@sam_data, stringsAsFactors = FALSE, check.names = FALSE) if ("sample" %in% colnames(data.obj$meta.dat)) { data.obj$meta.dat <- data.obj$meta.dat[, setdiff(colnames(data.obj$meta.dat), "sample"), drop = FALSE] diff --git a/R/mStat_normalize_data.R b/R/mStat_normalize_data.R index f41337d..6543071 100644 --- a/R/mStat_normalize_data.R +++ b/R/mStat_normalize_data.R @@ -141,6 +141,12 @@ mStat_compute_tmm_scale_factors <- function(otu_tab) { #' \item "DESeq": DESeq normalization #' \item "TMM": Trimmed Mean of M-values (requires edgeR) #' } +#' @param depth Target rarefaction depth for the "Rarefy" and "Rarefy-TSS" +#' methods. If NULL (default), the smallest per-sample total count is used. +#' Ignored by the non-rarefaction methods. +#' @param seed Integer seed used to make rarefaction subsampling reproducible +#' for the "Rarefy" and "Rarefy-TSS" methods. The global RNG state is restored +#' afterwards, so the caller's random stream is not affected. Default is 123. #' #' @return A list with normalized data object and scale factors. #' @@ -184,7 +190,8 @@ mStat_compute_tmm_scale_factors <- function(otu_tab) { mStat_normalize_data <- function(data.obj, method = c("Rarefy-TSS", "Rarefy", "TSS", "GMPR", "CSS", "DESeq", "TMM"), - depth = NULL) { + depth = NULL, + seed = 123) { # Validate input data structure # Ensuring the input is a list is crucial for maintaining the expected data format if (!is.list(data.obj)) { @@ -223,8 +230,13 @@ mStat_normalize_data <- if (all(round(colSums(otu_tab), 5) == 1)) { rarefied_otu_tab <- as.matrix(otu_tab) } else { - # Perform rarefaction using vegan package - rarefied_otu_tab <- t(vegan::rrarefy(t(otu_tab), sample = depth)) + # Perform rarefaction using vegan package. rrarefy() draws a random + # subsample, so wrap it in a local seed to make the result reproducible + # by default while leaving the caller's global RNG stream untouched. + rarefied_otu_tab <- mStat_with_local_seed( + seed, + t(vegan::rrarefy(t(otu_tab), sample = depth)) + ) } # Apply TSS normalization (convert to relative abundance) for Rarefy-TSS only diff --git a/man/mStat_normalize_data.Rd b/man/mStat_normalize_data.Rd index 31209e2..74298de 100644 --- a/man/mStat_normalize_data.Rd +++ b/man/mStat_normalize_data.Rd @@ -7,7 +7,8 @@ mStat_normalize_data( data.obj, method = c("Rarefy-TSS", "Rarefy", "TSS", "GMPR", "CSS", "DESeq", "TMM"), - depth = NULL + depth = NULL, + seed = 123 ) } \arguments{ @@ -42,6 +43,10 @@ Data objects can be created using converters like \item{depth}{Numeric value or NULL. Rarefaction depth for rarefaction workflows. If NULL, uses the minimum sample depth.} + +\item{seed}{Integer seed used to make rarefaction subsampling reproducible +for the "Rarefy" and "Rarefy-TSS" methods. The global RNG state is restored +afterwards, so the caller's random stream is not affected. Default is 123.} } \value{ A list with normalized data object and scale factors. From d39f7532862da3ceb178cc68f06692e05bec955e Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 04:12:20 +0000 Subject: [PATCH 4/6] Fix PERMANOVA distance/metadata misalignment in beta test 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 --- R/generate_beta_test_single.R | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/R/generate_beta_test_single.R b/R/generate_beta_test_single.R index 12b78de..ae55861 100644 --- a/R/generate_beta_test_single.R +++ b/R/generate_beta_test_single.R @@ -121,18 +121,39 @@ generate_beta_test_single <- function(data.obj, } else { # If distance object is provided, use it but ensure it matches the current data message("Using provided dist.obj...") - dist.obj <- mStat_subset_dist(dist.obj, rownames(data.obj$meta.dat)) # Check if all requested distances are available before filtering available_dists <- names(dist.obj) missing_dists <- setdiff(dist.name, available_dists) if (length(missing_dists) > 0) { - stop(paste("The following distances are not available in dist.obj:", + stop(paste("The following distances are not available in dist.obj:", paste(missing_dists, collapse = ", "))) } # Filter dist.obj to only include the distances specified in dist.name dist.obj <- dist.obj[dist.name] } + # GUniFrac::PermanovaG2 / vegan::adonis2 pair each distance matrix with the + # `data` rows by POSITION, not by sample name, and silently mishandle NA + # predictors. Restrict the metadata and every distance matrix to the same + # complete-case samples, in a consistent order, so the F/R2/p-values are + # computed on correctly paired data regardless of the incoming meta.dat row + # order. (Previously only the provided-dist.obj branch realigned, so a freshly + # computed distance with an unsorted meta.dat could be silently misaligned.) + analysis_vars <- c(adj.vars, group.var) + analysis_vars <- analysis_vars[!is.na(analysis_vars) & nzchar(analysis_vars)] + keep_samples <- rownames(data.obj$meta.dat) + if (length(analysis_vars) > 0) { + complete_idx <- stats::complete.cases( + data.obj$meta.dat[, analysis_vars, drop = FALSE] + ) + keep_samples <- keep_samples[complete_idx] + } + if (length(keep_samples) == 0) { + stop("No samples with complete data for the requested group/adjustment variables.", call. = FALSE) + } + data.obj$meta.dat <- data.obj$meta.dat[keep_samples, , drop = FALSE] + dist.obj <- mStat_subset_dist(dist.obj, keep_samples) + # Prepare for PERMANOVA analysis message("Running PermanovaG2 for all distances...") From bbd7f1af0e8f22ea43b223bc23417212f0e16121 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 05:10:24 +0000 Subject: [PATCH 5/6] Harden per-time alpha reference level and PERMANOVA repeated-measures 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 --- R/generate_alpha_per_time_test_long.R | 14 ++++++++++++++ R/generate_beta_test_single.R | 15 +++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/R/generate_alpha_per_time_test_long.R b/R/generate_alpha_per_time_test_long.R index 1aac8e6..7b7a1e4 100644 --- a/R/generate_alpha_per_time_test_long.R +++ b/R/generate_alpha_per_time_test_long.R @@ -144,6 +144,20 @@ generate_alpha_per_time_test_long <- function(data.obj, subset_data.obj <- subset_inputs$data.obj subset_alpha.obj <- subset_inputs$alpha.obj + # Pin the group reference level to the global one so the within-timepoint + # lm uses the same baseline that mStat_restructure_group_term_results + # applies as the label. Without this, lm falls back to the first level + # present in the subset, which can differ from `reference_level` when a + # timepoint is missing the global reference group, mislabelling the + # contrast. Only relevel when the reference level is actually present. + if (!is.null(reference_level) && + reference_level %in% as.character(subset_data.obj$meta.dat[[group.var]])) { + subset_data.obj$meta.dat[[group.var]] <- stats::relevel( + factor(subset_data.obj$meta.dat[[group.var]]), + ref = reference_level + ) + } + subset_meta_tab <- mStat_prepare_alpha_meta_tab( data.obj = subset_data.obj, vars = c(group.var, time.var, adj.vars) diff --git a/R/generate_beta_test_single.R b/R/generate_beta_test_single.R index ae55861..0fea333 100644 --- a/R/generate_beta_test_single.R +++ b/R/generate_beta_test_single.R @@ -113,6 +113,21 @@ generate_beta_test_single <- function(data.obj, # This step ensures we have the necessary distance matrices for the analysis if (!is.null(time.var) && !is.null(t.level)) { data.obj <- mStat_subset_by_meta_values(data.obj, time.var, t.level) + } else if (!is.null(time.var) && is.null(t.level) && + length(unique(data.obj$meta.dat[[time.var]])) > 1) { + # PERMANOVA permutes samples as exchangeable units. Pooling multiple + # timepoints without scoping to one (t.level = NULL) treats repeated + # measurements of the same subject as independent observations, which + # violates exchangeability and inflates the Type-I error rate. Warn rather + # than silently returning an anti-conservative test. + warning( + "time.var '", time.var, "' has multiple levels but t.level is NULL: ", + "PERMANOVA will pool all timepoints and treat repeated measures as ", + "independent, which inflates the Type-I error rate. Specify t.level to ", + "restrict the test to a single timepoint, or supply cross-sectional ", + "(one-row-per-subject) data.", + call. = FALSE + ) } if (is.null(dist.obj)) { From 0cfc3623fb68e1ac00379fdc4f06ec7810b22d56 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 3 Jun 2026 05:19:15 +0000 Subject: [PATCH 6/6] Fix overstated residual df in perform_lm_analysis (other data type) 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 --- R/generate_taxa_test_single.R | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/R/generate_taxa_test_single.R b/R/generate_taxa_test_single.R index 27a75b1..f7e958c 100644 --- a/R/generate_taxa_test_single.R +++ b/R/generate_taxa_test_single.R @@ -107,7 +107,13 @@ perform_lm_analysis <- function(feature.dat, df$reject <- df$padj <= alpha df$baseMean <- NA # Not applicable for linear models df$stat <- df$log2FoldChange / df$lfcSE - df$df <- nrow(meta.dat) - length(all.vars(as.formula(paste("~", formula)))) - 1 + # Residual degrees of freedom must account for the design-matrix rank: a + # k-level factor contributes k - 1 columns, not 1. Counting variables with + # all.vars() overstates df for any multi-level factor predictor. Derive it + # from the actual model matrix instead. (The reported p-values come from + # lm's summary and were already correct; only this df column was wrong.) + df$df <- nrow(meta.dat) - + ncol(stats::model.matrix(stats::as.formula(paste("~", formula)), data = meta.dat)) # Set rownames rownames(df) <- df$feature