### Install ANCOMBC using BiocManager

Source: https://github.com/frederickhuanglin/ancombc/blob/bugfix/README.md

Installs the latest release version of the ANCOMBC package using BiocManager. Ensure BiocManager is installed first.

```r
if (!requireNamespace("BiocManager", quietly=TRUE))
    install.packages("BiocManager")
BiocManager::install("ANCOMBC")
```

--------------------------------

### Install and Load lmerTest Package

Source: https://github.com/frederickhuanglin/ancombc/blob/bugfix/README.md

Installs and loads the lmerTest package in R, which is useful for debugging mixed-effects models.

```r
install.packages("lmerTest")
library(lmerTest)
```

--------------------------------

### Basic Two-Group Comparison with Sensitivity Analysis using ancombc2

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Performs a differential abundance analysis with a sensitivity analysis to assess robustness to zero imputation. Requires loading ANCOMBC and microbiome packages, and example data. Agglomerates data to the Family level and uses specified fixed and no random effects.

```r
library(ANCOMBC)
library(microbiome)

# ── Load example data ──────────────────────────────────────────────────────────
data(dietswap, package = "microbiome")

# ── Basic two-group comparison with sensitivity analysis ──────────────────────
set.seed(123)
out <- ancombc2(
  data          = dietswap,
  tax_level     = "Family",           # agglomerate to Family level
  fix_formula   = "nationality + timepoint + bmi_group",
  rand_formula  = NULL,               # no random effects
  p_adj_method  = "holm",
  pseudo        = 0,                  # no global pseudo-count
  pseudo_sens   = TRUE,               # enable sensitivity analysis (recommended)
  prv_cut       = 0.10,               # exclude taxa present in < 10 % of samples
  lib_cut       = 1000,               # exclude samples with < 1000 reads
  s0_perc       = 0.05,               # SAM fudge factor at 5th percentile of SEs
  group         = "bmi_group",        # discrete grouping variable
  struc_zero    = TRUE,               # detect structural zeros
  neg_lb        = TRUE,               # use asymptotic lower bound
  alpha         = 0.05,
  n_cl          = 1,
  verbose       = TRUE,
  global        = TRUE,               # omnibus test across bmi_group levels
  pairwise      = TRUE,               # all pairwise comparisons
  dunnet        = TRUE,               # compare each group to reference
  trend         = TRUE,               # monotone trend test
  iter_control  = list(tol = 1e-2, max_iter = 20, verbose = FALSE),
  em_control    = list(tol = 1e-5, max_iter = 100),
  lme_control   = lme4::lmerControl(),
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
  trend_control = list(
    contrast = list(matrix(c(1, 0, -1, 1), nrow = 2, byrow = TRUE)),
    node     = list(2),
    B        = 100
  )
)

# ── Inspect primary results ────────────────────────────────────────────────────
res_prim   <- out$res
# Columns: taxon, lfc_*, se_*, W_*, p_*, q_*, diff_*, passed_ss_*, diff_robust_*
sig_robust <- res_prim[res_prim$diff_robust_bmi_groupoverweight == TRUE, ]
catalog("Robustly DA taxa (overweight vs. lean):")
print(sig_robust[, c("taxon", "lfc_bmi_groupoverweight", "q_bmi_groupoverweight")])

# ── Global test ───────────────────────────────────────────────────────────────
res_global <- out$res_global          # W, p_val, q_val, diff_abn, passed_ss
catalog("\nGlobally DA taxa:")
print(res_global[res_global$diff_robust_abn == TRUE, ])

# ── Pairwise directional test ─────────────────────────────────────────────────
res_pair   <- out$res_pair            # lfc/se/W/p/q/diff/passed_ss per pair

# ── Dunnett's test ────────────────────────────────────────────────────────────
res_dunn   <- out$res_dunn

# ── Trend test ────────────────────────────────────────────────────────────────
res_trend  <- out$res_trend

```

--------------------------------

### ANCOM with Repeated Measures

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Example of using ANCOM for repeated-measures data by specifying a random effects formula. This allows for modeling within-subject correlations.

```r
# ── Repeated-measures example ─────────────────────────────────────────────────
# out_lme <- ancom(
#   data         = pseq_longitudinal,
#   main_var     = "treatment",
#   adj_formula  = "age",
#   rand_formula = "(1|subject_id)",    # random intercept per subject
#   lme_control  = lme4::lmerControl(),
#   struc_zero   = FALSE
# )
```

--------------------------------

### Load ANCOMBC and access help pages

Source: https://github.com/frederickhuanglin/ancombc/blob/bugfix/README.md

Loads the ANCOMBC library and provides quick access to help documentation for key functions like ancombc2, ancombc, ancom, secom_linear, and secom_dist.

```r
library(ANCOMBC)
?ancombc2
?ancombc 
?ancom
?secom_linear
?secom_dist
```

--------------------------------

### Create Matrix in R

Source: https://github.com/frederickhuanglin/ancombc/blob/bugfix/README.md

This R code snippet demonstrates how to create a matrix, which can be used for representing data in statistical analyses.

```R
matrix(c(1, 0, 0, 0,
  -1, 1, 0, 0,
   0, -1, 1, 0,
   0, 0, -1, 1),
   nrow = 4, 
   byrow = TRUE)
```

--------------------------------

### Run ANCOMBC for Bias-Corrected Analysis

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Use this snippet to perform differential abundance analysis with bias correction. Ensure necessary libraries are loaded and data is prepared. The function supports various model formulas and adjustment methods.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)   # baseline samples only

set.seed(123)
out <- ancombc(
  data         = pseq,
  tax_level    = "Family",
  formula      = "age + nationality + bmi_group",  # all covariates including group
  p_adj_method = "holm",
  prv_cut      = 0.10,
  lib_cut      = 1000,
  group        = "bmi_group",     # required for global test and struc_zero
  struc_zero   = TRUE,
  neg_lb       = FALSE,
  tol          = 1e-5,
  max_iter     = 100,
  conserve     = TRUE,            # conservative SE (recommended for small n)
  alpha        = 0.05,
  global       = TRUE,
  n_cl         = 1,
  verbose      = TRUE
)
```

--------------------------------

### Compute Sparse Linear Correlations with secom_linear

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Use `secom_linear` to compute a bias-corrected, sparse taxon-taxon correlation matrix. It supports Pearson or Spearman correlation and allows for multi-ecosystem analyses by accepting a named list of datasets. Ensure necessary libraries are loaded and data is prepared in a phyloseq object.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)

set.seed(123)
res_linear <- secom_linear(
  data        = list(gut = pseq),   # wrap in list; use list(gut=, tongue=) for multi-ecosystem
  taxa_are_rows = TRUE,
  tax_level   = "Phylum",
  pseudo      = 0,
  prv_cut     = 0.50,               # require 50 % prevalence for correlation stability
  lib_cut     = 1000,
  corr_cut    = 0.50,               # flag taxa highly correlated with sampling fraction
  wins_quant  = c(0.05, 0.95),      # Winsorise at 5th / 95th percentile
  method      = "pearson",          # or "spearman"
  soft        = FALSE,              # hard thresholding (TRUE = soft / LASSO-like)
  alpha_grid  = 0,                  # L1 penalty; 0 = no regularisation on top of threshold
  thresh_len  = 20,                 # grid points for CV threshold search
  n_cv        = 10,                 # 10-fold CV
  thresh_hard = 0.3,                # zero out |r| <= 0.3 after CV
  max_p       = 0.005,              # zero out pairs with p > 0.005
  n_cl        = 2,                  # parallel workers
  verbose     = TRUE
)

# ── Outputs ────────────────────────────────────────────────────────────────────
corr      <- res_linear$corr        # full sample Pearson/Spearman matrix
corr_p    <- res_linear$corr_p      # corresponding p-value matrix
corr_th   <- res_linear$corr_th     # sparse matrix after threshold
corr_fl   <- res_linear$corr_fl     # sparse matrix after p-value filter
corr_reg  <- res_linear$corr_reg    # positive-definite regularised matrix
thresh_opt <- res_linear$thresh_opt # CV-selected threshold

# Visualise top correlations
library(igraph)
adj <- abs(corr_fl) > 0
diag(adj) <- FALSE
g <- graph_from_adjacency_matrix(adj, mode = "undirected", weighted = TRUE)
plot(g, main = "Phylum co-occurrence network (secom_linear)")

# ── Multi-ecosystem example (gut + tongue) ────────────────────────────────────
# res_multi <- secom_linear(
#   data = list(gut = pseq_gut, tongue = pseq_tongue),
#   tax_level = "Genus",
#   prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
#   method = "pearson", thresh_hard = 0.3, max_p = 0.005
# )
# cross_corr <- res_multi$corr_th["gut - Bacteroides", "tongue - Veillonella"]

```

--------------------------------

### Extract ANCOM Results

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Retrieve results from the ANCOM analysis, including the W statistic, detection flags for various cutoffs, and raw or adjusted p-values for pairwise regressions. Structural zero indicators are also provided if enabled.

```r
res <- out$res
# Columns: taxon, W, detected_0.9, detected_0.8, detected_0.7, detected_0.6
da_0.7 <- res[res$detected_0.7 == TRUE, ]
cat("DA taxa at W > 0.7:\n")
print(da_0.7)

# Raw pairwise test results
beta_data <- out$beta_data  # n_tax × n_tax coefficients
p_data    <- out$p_data     # n_tax × n_tax raw p-values
q_data    <- out$q_data     # n_tax × n_tax adjusted p-values
zero_ind  <- out$zero_ind   # structural zeros (if struc_zero = TRUE)
```

--------------------------------

### Simulate Microbial Data with sim_plnm

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Generate synthetic microbiome count data using `sim_plnm` based on a Poisson lognormal model estimated from a real dataset. This function is useful for power analysis and benchmarking DA methods. Ensure the ANCOMBC library is loaded and a template dataset is available.

```r
library(ANCOMBC)

# ── Using the built-in QMP dataset as template ─────────────────────────────────
data(QMP)               # 106 samples × 91 OTUs (samples in rows)
dim(QMP)                # [1] 106  91

set.seed(42)
sim_counts <- sim_plnm(
  abn_table    = QMP,
  taxa_are_rows = FALSE,   # QMP has samples in rows
  prv_cut      = 0.05,     # keep taxa present in >= 5 % of samples
  n            = 200,      # generate 200 simulated samples
  lib_mean     = 1e8,      # mean library size
  disp         = 0.5       # dispersion of library size (NegBin size parameter)
)
# sim_counts: taxa × samples matrix
rownames(sim_counts) <- paste0("Taxon",   seq_len(nrow(sim_counts)))
colnames(sim_counts) <- paste0("Sample",  seq_len(ncol(sim_counts)))
dim(sim_counts)          # e.g. [1] 47 200  (after prevalence filtering)
range(sim_counts)        # non-negative integers
```

```r
# ── Downstream: run ANCOM-BC2 on simulated data ────────────────────────────────
meta_sim <- data.frame(
  group = factor(rep(c("ctrl", "trt"), each = 100)),
  row.names = colnames(sim_counts)
)

set.seed(1)
out_sim <- ancombc2(
  data        = sim_counts,
  meta_data   = meta_sim,
  fix_formula = "group",
  p_adj_method = "BH",
  prv_cut     = 0.10,
  pseudo_sens = FALSE,    # skip sensitivity for speed in simulation
  verbose     = FALSE
)
sum(out_sim$res$diff_group_trt, na.rm = TRUE)  # number of DA taxa detected
```

--------------------------------

### Run ANCOM for Non-Parametric Differential Abundance

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Execute the original ANCOM method for differential abundance testing. This function is suitable for both linear models and linear mixed-effects models. Specify the main variable of interest and covariates to adjust for.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)

set.seed(123)
out <- ancom(
  data         = pseq,
  tax_level    = "Family",
  p_adj_method = "holm",
  prv_cut      = 0.10,
  lib_cut      = 1000,
  main_var     = "bmi_group",          # variable of interest (not in adj_formula)
  adj_formula  = "age + nationality",  # covariates to adjust for
  rand_formula = NULL,                 # e.g. "(1|subject_id)" for repeated measures
  lme_control  = NULL,
  struc_zero   = TRUE,
  neg_lb       = TRUE,
  alpha        = 0.05,
  n_cl         = 1,
  verbose      = TRUE
)
```

--------------------------------

### Extract ANCOMBC Results

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Access various result tables from the ANCOMBC output, including log-fold changes, standard errors, p-values, q-values, and differentially abundant taxa indicators. Global test results and structural zero information are also available.

```r
# ── Primary results ────────────────────────────────────────────────────────────
res      <- out$res
lfc      <- res$lfc        # data.frame: taxon + lfc per covariate
se       <- res$se
W        <- res$W
p_val    <- res$p_val
q_val    <- res$q_val
diff_abn <- res$diff_abn   # logical: TRUE if q < alpha

# Differentially abundant taxa for bmi_groupoverweight vs. reference
da_idx <- which(diff_abn$bmi_groupoverweight)
print(lfc[da_idx, c("taxon", "bmi_groupoverweight")])

# ── Global test results ────────────────────────────────────────────────────────
res_global <- out$res_global    # W, p_val, q_val, diff_abn

# ── Structural zeros ──────────────────────────────────────────────────────────
zero_ind <- out$zero_ind        # logical data.frame; TRUE = structural zero in that group

# ── Sampling fractions (log scale) ────────────────────────────────────────────
samp_frac <- out$samp_frac
```

--------------------------------

### secom_linear

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Computes a bias-corrected, sparse taxon–taxon correlation matrix using Pearson or Spearman correlation. It estimates sample-specific biases, subtracts them from log-abundances, Winsorises outliers, and then selects an optimal threshold via cross-validation to induce sparsity. A false-positive guard removes taxon pairs that both correlate strongly with the estimated bias. Supports multi-ecosystem analyses.

```APIDOC
## `secom_linear` — Sparse Estimation of Linear Correlations among Microbiomes

Computes a bias-corrected, sparse taxon–taxon correlation matrix using Pearson or Spearman correlation. It first estimates sample-specific biases (sampling fractions) and subtracts them from log-abundances, Winsorises outliers, then selects an optimal hard/soft threshold via cross-validation to induce sparsity. A false-positive guard removes taxon pairs that both correlate strongly with the estimated bias. Supports multi-ecosystem analyses by accepting a named list of datasets.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)

set.seed(123)
res_linear <- secom_linear(
  data        = list(gut = pseq),   # wrap in list; use list(gut=, tongue=) for multi-ecosystem
  taxa_are_rows = TRUE,
  tax_level   = "Phylum",
  pseudo      = 0,
  prv_cut     = 0.50,               # require 50 % prevalence for correlation stability
  lib_cut     = 1000,
  corr_cut    = 0.50,               # flag taxa highly correlated with sampling fraction
  wins_quant  = c(0.05, 0.95),      # Winsorise at 5th / 95th percentile
  method      = "pearson",          # or "spearman"
  soft        = FALSE,              # hard thresholding (TRUE = soft / LASSO-like)
  alpha_grid  = 0,                  # L1 penalty; 0 = no regularisation on top of threshold
  thresh_len  = 20,                 # grid points for CV threshold search
  n_cv        = 10,                 # 10-fold CV
  thresh_hard = 0.3,                # zero out |r| <= 0.3 after CV
  max_p       = 0.005,              # zero out pairs with p > 0.005
  n_cl        = 2,                  # parallel workers
  verbose     = TRUE
)

# ── Outputs ────────────────────────────────────────────────────────────────────
corr      <- res_linear$corr        # full sample Pearson/Spearman matrix
corr_p    <- res_linear$corr_p      # corresponding p-value matrix
corr_th   <- res_linear$corr_th     # sparse matrix after threshold
corr_fl   <- res_linear$corr_fl     # sparse matrix after p-value filter
corr_reg  <- res_linear$corr_reg    # positive-definite regularised matrix
thresh_opt <- res_linear$thresh_opt # CV-selected threshold

# Visualise top correlations
library(igraph)
adj <- abs(corr_fl) > 0
diag(adj) <- FALSE
g <- graph_from_adjacency_matrix(adj, mode = "undirected", weighted = TRUE)
plot(g, main = "Phylum co-occurrence network (secom_linear)")

# ── Multi-ecosystem example (gut + tongue) ────────────────────────────────────
# res_multi <- secom_linear(
#   data = list(gut = pseq_gut, tongue = pseq_tongue),
#   tax_level = "Genus",
#   prv_cut = 0.5, lib_cut = 1000, corr_cut = 0.5,
#   method = "pearson", thresh_hard = 0.3, max_p = 0.005
# )
# cross_corr <- res_multi$corr_th["gut - Bacteroides", "tongue - Veillonella"]
```
```

--------------------------------

### Validate and Normalize Input Data with data_sanity_check

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Use `data_sanity_check` to validate and normalize various input data types (phyloseq, TreeSummarizedExperiment, matrix, data.frame) before ANCOMBC analysis. It checks formula variables, sample sizes, and configuration constraints. Ensure all required libraries are loaded.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")

check <- data_sanity_check(
  data        = atlas1006,
  tax_level   = "Family",
  fix_formula = "age + sex + bmi_group",
  group       = "bmi_group",
  struc_zero  = TRUE,
  global      = TRUE,
  pairwise    = TRUE,
  dunnet      = TRUE,
  trend       = FALSE,
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
  verbose     = TRUE
)
# Expected output (verbose messages):
# Checking the input data type ... PASS
# The number of groups of interest is: 3
# The sample size per group is: lean = 161, overweight = 190, obese = 178
# Checking other arguments ... PASS

feature_table           <- check$feature_table            # raw OTU/ASV matrix
feature_table_aggregate <- check$feature_table_aggregate  # possibly agglomerated
meta_data               <- check$meta_data                # cleaned metadata
# Booleans may be silently downgraded (e.g. global=FALSE if n_groups < 3)
cat("global enabled:", check$global, "\n")
cat("pairwise enabled:", check$pairwise, "\n")
```

```r
# ── Using matrix/data.frame input ─────────────────────────────────────────────
count_mat <- matrix(rpois(50 * 20, 10), nrow = 50, ncol = 20,
                    dimnames = list(paste0("Taxon", 1:50),
                                    paste0("Sample", 1:20)))
meta      <- data.frame(group = rep(c("A", "B"), each = 10),
                        row.names = paste0("Sample", 1:20))
check2 <- data_sanity_check(
  data       = count_mat,
  meta_data  = meta,
  fix_formula = "group",
  group      = "group",
  verbose    = FALSE
)
```

--------------------------------

### Compute Sparse Distance Correlations with secom_dist

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Use `secom_dist` to compute sparse distance correlations, capturing non-linear dependencies between taxa. This function performs bias correction and uses permutation replicates for p-value estimation. It sparsifies the results by p-value filtering.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)

set.seed(123)
res_dist <- secom_dist(
  data          = list(pseq),
  taxa_are_rows = TRUE,
  tax_level     = "Phylum",
  pseudo        = 0,
  prv_cut       = 0.50,
  lib_cut       = 1000,
  corr_cut      = 0.50,
  wins_quant    = c(0.05, 0.95),
  R             = 1000,            # permutation replicates for p-value estimation
  thresh_hard   = 0.3,             # zero out dcorr <= 0.3
  max_p         = 0.005,
  n_cl          = 2,
  verbose       = TRUE
)

dcorr    <- res_dist$dcorr      # full distance correlation matrix (n_taxa × n_taxa)
dcorr_p  <- res_dist$dcorr_p   # permutation p-value matrix
dcorr_fl <- res_dist$dcorr_fl  # sparse matrix after p-value filtering
y_hat    <- res_dist$y_hat      # bias-corrected log-abundance matrix

# Strongest non-linear associations at Phylum level
upper_tri <- upper.tri(dcorr_fl)
pairs_df <- data.frame(
  taxon1 = rownames(dcorr_fl)[row(dcorr_fl)[upper_tri]],
  taxon2 = colnames(dcorr_fl)[col(dcorr_fl)[upper_tri]],
  dcorr  = dcorr_fl[upper_tri]
)
pairs_df <- pairs_df[pairs_df$dcorr > 0, ]
pairs_df <- pairs_df[order(-pairs_df$dcorr), ]
head(pairs_df)

```

--------------------------------

### secom_dist

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

Extends SECOM to distance correlation, capturing any dependence (linear, monotone, or nonlinear) between two random variables. After bias-correction preprocessing, it computes all pairwise distance correlations and their permutation p-values, then sparsifies by p-value filtering. Distance correlation is always non-negative; values near zero indicate independence.

```APIDOC
## `secom_dist` — Sparse Estimation of Distance Correlations among Microbiomes

Extends SECOM to distance correlation, a measure that captures any dependence (linear, monotone, or nonlinear) between two random variables. After the same bias-correction preprocessing as `secom_linear`, it computes all pairwise distance correlations and their permutation p-values, then sparsifies by p-value filtering. Distance correlation is always non-negative; values near zero indicate independence.

```r
library(ANCOMBC)
library(microbiome)

data(atlas1006, package = "microbiome")
pseq <- phyloseq::subset_samples(atlas1006, time == 0)

set.seed(123)
res_dist <- secom_dist(
  data          = list(pseq),
  taxa_are_rows = TRUE,
  tax_level     = "Phylum",
  pseudo        = 0,
  prv_cut       = 0.50,
  lib_cut       = 1000,
  corr_cut      = 0.50,
  wins_quant    = c(0.05, 0.95),
  R             = 1000,            # permutation replicates for p-value estimation
  thresh_hard   = 0.3,             # zero out dcorr <= 0.3
  max_p         = 0.005,
  n_cl          = 2,
  verbose       = TRUE
)

dcorr    <- res_dist$dcorr      # full distance correlation matrix (n_taxa × n_taxa)
dcorr_p  <- res_dist$dcorr_p   # permutation p-value matrix
dcorr_fl <- res_dist$dcorr_fl  # sparse matrix after p-value filtering
y_hat    <- res_dist$y_hat      # bias-corrected log-abundance matrix

# Strongest non-linear associations at Phylum level
upper_tri <- upper.tri(dcorr_fl)
pairs_df <- data.frame(
  taxon1 = rownames(dcorr_fl)[row(dcorr_fl)[upper_tri]],
  taxon2 = colnames(dcorr_fl)[col(dcorr_fl)[upper_tri]],
  dcorr  = dcorr_fl[upper_tri]
)
pairs_df <- pairs_df[pairs_df$dcorr > 0, ]
pairs_df <- pairs_df[order(-pairs_df$dcorr), ]
head(pairs_df)
```
```

--------------------------------

### ancombc2 — Analysis of Compositions of Microbiomes with Bias Correction 2

Source: https://context7.com/frederickhuanglin/ancombc/llms.txt

The flagship DA function. It estimates and removes sample-specific and taxon-specific biases via an iterative MLE / E-M algorithm, then performs inference using a SAM-inspired stabilised test statistic. Beyond the primary two-group or continuous-covariate test it optionally runs a global omnibus test, pairwise directional test (with mdFDR control), Dunnett's type of comparison, and a user-specified trend test. A built-in sensitivity analysis reruns the full pipeline with three different pseudo-count values (0.1, 0.5, 1) to flag results that are not robust to zero imputation.

```APIDOC
## `ancombc2` — Analysis of Compositions of Microbiomes with Bias Correction 2

### Description
The flagship DA function. It estimates and removes sample-specific and taxon-specific biases via an iterative MLE / E-M algorithm, then performs inference using a SAM-inspired stabilised test statistic. Beyond the primary two-group or continuous-covariate test it optionally runs a global omnibus test, pairwise directional test (with mdFDR control), Dunnett's type of comparison, and a user-specified trend test. A built-in sensitivity analysis reruns the full pipeline with three different pseudo-count values (0.1, 0.5, 1) to flag results that are not robust to zero imputation.

### Method Signature
```r
ancombc2(
  data,
  tax_level = "",
  fix_formula,
  rand_formula = NULL,
  p_adj_method = "holm",
  pseudo = 0,
  pseudo_sens = TRUE,
  prv_cut = 0.10,
  lib_cut = 1000,
  s0_perc = 0.05,
  group,
  struc_zero = TRUE,
  neg_lb = TRUE,
  alpha = 0.05,
  n_cl = 1,
  verbose = TRUE,
  global = TRUE,
  pairwise = TRUE,
  dunnet = TRUE,
  trend = TRUE,
  iter_control = list(tol = 1e-2, max_iter = 20, verbose = FALSE),
  em_control = list(tol = 1e-5, max_iter = 100),
  lme_control = lme4::lmerControl(),
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
  trend_control = list(
    contrast = list(matrix(c(1, 0, -1, 1), nrow = 2, byrow = TRUE)),
    node     = list(2),
    B        = 100
  )
)
```

### Parameters

*   **data** (matrix/data.frame/phyloseq/TreeSummarizedExperiment) - Input microbiome data.
*   **tax_level** (string) - Taxonomic level to agglomerate to (e.g., "Family").
*   **fix_formula** (string) - Fixed effects formula for the model.
*   **rand_formula** (string, optional) - Random effects formula for the model.
*   **p_adj_method** (string) - Method for p-value adjustment (e.g., "holm").
*   **pseudo** (numeric) - Value for global pseudo-count.
*   **pseudo_sens** (boolean) - Whether to enable sensitivity analysis for pseudo-counts.
*   **prv_cut** (numeric) - Prevalence cutoff for taxa.
*   **lib_cut** (numeric) - Library size cutoff for samples.
*   **s0_perc** (numeric) - SAM fudge factor percentile.
*   **group** (string) - Name of the discrete grouping variable.
*   **struc_zero** (boolean) - Whether to detect structural zeros.
*   **neg_lb** (boolean) - Whether to use the asymptotic lower bound.
*   **alpha** (numeric) - Significance level.
*   **n_cl** (numeric) - Number of clusters for parallel processing.
*   **verbose** (boolean) - Whether to print verbose output.
*   **global** (boolean) - Whether to perform a global omnibus test.
*   **pairwise** (boolean) - Whether to perform pairwise directional tests.
*   **dunnet** (boolean) - Whether to perform Dunnett's test.
*   **trend** (boolean) - Whether to perform a trend test.
*   **iter_control** (list) - Control parameters for the iterative algorithm.
*   **em_control** (list) - Control parameters for the EM algorithm.
*   **lme_control** (list) - Control parameters for linear mixed-effects models.
*   **mdfdr_control** (list) - Control parameters for the moderated FDR method.
*   **trend_control** (list) - Control parameters for the trend test.

### Request Example
```r
library(ANCOMBC)
library(microbiome)

data(dietswap, package = "microbiome")

set.seed(123)
out <- ancombc2(
  data          = dietswap,
  tax_level     = "Family",
  fix_formula   = "nationality + timepoint + bmi_group",
  rand_formula  = NULL,
  p_adj_method  = "holm",
  pseudo        = 0,
  pseudo_sens   = TRUE,
  prv_cut       = 0.10,
  lib_cut       = 1000,
  s0_perc       = 0.05,
  group         = "bmi_group",
  struc_zero    = TRUE,
  neg_lb        = TRUE,
  alpha         = 0.05,
  n_cl          = 1,
  verbose       = TRUE,
  global        = TRUE,
  pairwise      = TRUE,
  dunnet        = TRUE,
  trend         = TRUE,
  iter_control  = list(tol = 1e-2, max_iter = 20, verbose = FALSE),
  em_control    = list(tol = 1e-5, max_iter = 100),
  lme_control   = lme4::lmerControl(),
  mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
  trend_control = list(
    contrast = list(matrix(c(1, 0, -1, 1), nrow = 2, byrow = TRUE)),
    node     = list(2),
    B        = 100
  )
)
```

### Response

*   **res** (data.frame) - Primary results including taxon, LFC, SE, W statistic, p-value, adjusted p-value, difference indicator, and robust difference indicator.
*   **res_global** (data.frame) - Results from the global omnibus test.
*   **res_pair** (data.frame) - Results from pairwise directional tests.
*   **res_dunn** (data.frame) - Results from Dunnett's test.
*   **res_trend** (data.frame) - Results from the trend test.

### Response Example (Primary Results)
```r
# Inspect primary results
res_prim <- out$res
sig_robust <- res_prim[res_prim$diff_robust_bmi_groupoverweight == TRUE, ]
cat("Robustly DA taxa (overweight vs. lean):\n")
print(sig_robust[, c("taxon", "lfc_bmi_groupoverweight", "q_bmi_groupoverweight")])
```

### Response Example (Global Test)
```r
# Global test
res_global <- out$res_global
cat("\nGlobally DA taxa:\n")
print(res_global[res_global$diff_robust_abn == TRUE, ])
```
```

--------------------------------

### ANCOM-BC2 Trend Test Contrast Matrix in R

Source: https://github.com/frederickhuanglin/ancombc/blob/bugfix/README.md

This R code defines a contrast matrix for a trend test with 5 ordered categories. It is used to estimate differences between ordered groups and test for a monotonic trend.

```R
contrasts <- matrix(c(0, 1, 1, 1, 1,
        0, 0, 1, 1, 1,
        0, 0, 0, 1, 1,
        0, 0, 0, 0, 1),
        nrow = 5, ncol = 4)
rownames(contrasts) <- c("A", "B", "C", "D", "E")
colnames(contrasts) <- c("B-A", "C-A", "D-A", "E-A")
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.