Hello, I think I’ve found a subtle bug in maaslin3 `maaslin_contrast_test()`.
For a bit of context, I had been running maaslin3 multiple times across many contrasts and then removing the duplicate contrasts after the fact. This was somewhat inefficient and I through it would be best to use the maaslin_contrast_test to speed that up and simplify the process.
After switching from the loop-dedup method to maaslin_contrast_test and testing the output, I found that everything is the same except for the pvalues which are much smaller for the same constrasts and same formula.
The root cause is that maaslin_contrast_test uses multcomp::glht z-test in its median-comparison step when everywhere else maaslin3 uses lmerTest::contest t-test (as far as I can tell). The z-test makes p-values much smaller than the t-test.
Also, looking at the documentation for maaslin_contrast_test it states that lmerTest::contest should be used for mixed effects linear.
Look more closely in R/utility_scripts.R, this bug is only happening when median_comparison_abundance = TRUE which is my case.
Here is a reprex using maaslin3’s datasets, I create what I would have thought to be the same pvalues using maaslin3 and maaslin_contrast_test, but they are not the same
library(maaslin3)
taxa <- read.csv(
system.file("extdata", "HMP2_taxonomy.tsv", package = "maaslin3"),
sep = "\\t",
row.names = 1
)
meta <- read.csv(
system.file("extdata", "HMP2_metadata.tsv", package = "maaslin3"),
sep = "\\t",
row.names = 1
)
\# 4 participants per diagnosis x 4 samples: few groups, so small Satterthwaite df
meta <- meta\[
ave(seq_along(meta$participant_id), meta$participant_id, FUN = seq_along) <=
4,
\]
ids <- unlist(lapply(split(meta, meta$diagnosis), function(m) {
head(names(which(table(m$participant_id) == 4)), 4)
}))
meta <- meta\[meta$participant_id %in% ids, \]
meta$diagnosis <- factor(meta$diagnosis, levels = c("nonIBD", "UC", "CD"))
taxa <- taxa\[rownames(meta), colMeans(taxa\[rownames(meta), \] > 0) > 0.5\]
set.seed(1)
fit <- maaslin3(
taxa,
meta,
output = tempfile(),
formula = "\~ diagnosis + (1 | participant_id)",
evaluate_only = "abundance",
warn_prevalence = FALSE,
save_models = TRUE,
plot_summary_plot = FALSE,
plot_associations = FALSE,
verbosity = "ERROR"
)
\# Contrast that is exactly the diagnosisCD coefficient: same fit, same median, same null
L <- matrix(1, 1, 1, dimnames = list("CD_vs_nonIBD", "diagnosisCD"))
ct <- maaslin_contrast_test(fit, L, evaluate_only = "abundance")
main <- subset(
fit$fit_data_abundance$results,
name == "diagnosisCD",
c(feature, coef, pval_individual)
)
both <- merge(
main,
ct$fit_data_abundance$results\[, c("feature", "pval_individual")\],
by = "feature",
suffixes = c("\_maaslin3", "\_contrast_test")
)
both <- both\[order(both$pval_individual_maaslin3), \]
both$is_same <- abs(both$pval_individual_maaslin3 - both$pval_individual_contrast_test) < 0.0001
print(both)
You can also see the issue if you put a browser() at line 979 of R/utility_scripts.R, and then run the reprex. It will show that the variable fit is class lmerModLmerTest and is tested using multcomp::glht.
Proposed fix
I think the fix is pretty straight forward, it would just be replaceing the multcomp::glht with lmerTest::contest in else { ... } at line 978 (link).
Replace this:
summary_out <- summary(multcomp::glht(
fit,
linfct = contrast_vec,
rhs = offsets_to_test\[row_index\]
))$test
with this:
pval <- lmerTest::contest(fit,
matrix(contrast_vec, TRUE),
rhs = offsets_to_test\[row_index\])\[\['Pr(>F)'\]\]
coef <- contrast_vec %\*% lme4::fixef(fit)
sigma <- sqrt((contrast_vec %\*% vcov(fit) %\*%
t(contrast_vec))\[1, 1\])
summary_out <- list(pvalues = pval,
coefficients = coef,
sigma = sigma)
With this change, the reprex above gives identical p-values from maaslin3() and maaslin_contrast_test() for single-coefficient contrasts. I also verified with two other datasets, that I can’t share.
Also, happy to open a PR if that helps things along!
Thanks again for your package. I hope all this makes sense, I’m more of an R developer than a bioinformatics person so please let me kow if anything is not clear.
Cheers,
Eli Daniels