-
Notifications
You must be signed in to change notification settings - Fork 31
Description
Hello ANCOM-BC experts - I’d appreciate advice on how to parameterize ANCOM-BC2 so pairwise contrasts for all my requested comparisons show up reproducibly (I’m seeing single-index columns referencing one baseline and missing the two-index pair columns I expect).
Short experimental design
Subset: ps_Chap3_DA_ITS_AT = samples where (Arrival_time == "CA" & Month == "12M") | (Arrival_time == "LA" & Month == "24M") and Treatment %in% c("K","M","KM").
N: 30 samples (6 Treat_AT groups × 5 each).
I am trying to study within-treatment arrival-time comparisons (eg. K treatment CA concurrent-arrival vs K treatment late-arrival). Intially I tried to run Treatment * Arrival_time + Block but model failed. So I combined Treatment & Arrival into a variable and ran Treat_AT + Block instead:
Treat_AT = paste(Treatment, Arrival_time, sep = "_") with enforced levels: K_CA, K_LA, KM_CA, KM_LA, M_CA, M_LA.
Block is Block 1 to 5 (was supposed to be covariate as Block were found to be significant in beta diversity analysis)
Exact ANCOM-BC2 call / parameters (what I used)
res <- ancombc2(
data = ps_Chap3_DA_ITS_AT,
tax_level = <NULL or "Phylum"/"Family"/"Genus">,
fix_formula = "Treat_AT + Block",
rand_formula = NULL,
group = "Treat_AT",
p_adj_method = "BH",
prv_cut = 0.10,
lib_cut = 1000,
s0_perc = 0.05,
pseudo_sens = TRUE,
struc_zero = TRUE,
neg_lb = TRUE,
dunnet = FALSE,
alpha = 0.05,
n_cl = 1,
iter_control = list(tol = 1e-2, max_iter = 20, verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
global = TRUE,
pairwise = TRUE
)
Contrasts I specifically want (within-treatment arrival-time comparisons)
K_CA vs K_LA
M_CA vs M_LA
KM_CA vs KM_LA
(Under my enforced ordering these map to Treat_AT1 vs Treat_AT2, Treat_AT5 vs Treat_AT6, Treat_AT3 vs Treat_AT4.)
Problem / question (brief)
res$res_pair shows lfc_Treat_AT1..lfc_Treat_AT5 and pairwise columns like lfc_Treat_AT2_Treat_AT1, but no Treat_AT6 token (so the M_CA vs M_LA pairwise column such as q_Treat_AT6_Treat_AT5 is missing). I did not set dunnet = TRUE or an explicit reference manually; I forced the factor levels in phyloseq before running.
Questions
Is it expected ANCOM-BC2 parameterizes with a single-reference index even when pairwise = TRUE?
Would releveling Treat_AT (so a different reference) force explicit two-index pairwise columns for all contrasts? Example I might try:
sample_data(ps_Chap3_DA_ITS_AT)$Treat_AT <- relevel(sample_data(ps_Chap3_DA_ITS_AT)$Treat_AT, ref="K_CA")
re-run
Alternatives: should I set group = NULL and compute contrasts externally, or change group usage? Which guarantees explicit two-index pairwise columns?
Any quick snippet to detect which reference/parameterization ANCOM used so I can re-run only affected models?