Dirichlet-Multinomial Mixture Models for eDNA Metabarcoding Community Structure
eDNAstructure is an R package for fitting Bayesian Dirichlet-Multinomial Mixture (DMM) models to environmental DNA (eDNA) read count data from metabarcoding surveys. Given a sample × taxon count matrix and optional environmental covariates, the model identifies latent ecological communities, estimates their taxonomic compositions, and quantifies how environmental gradients drive community membership — all within a fully Bayesian framework with principled uncertainty quantification.
Stan compiles models to C++ and requires a toolchain on your machine:
- Windows: Install Rtools
- macOS: Run
xcode-select --installin Terminal - Linux: Install
build-essential(Ubuntu/Debian) or equivalent
install.packages("rstan")Verify it works before proceeding:
library(rstan)
example(stan_model, package = "rstan", run.dontrun = TRUE)If you see sampling output without errors, Stan is ready. Full guide: https://mc-stan.org/rstan/articles/rstan.html
install.packages("remotes")
remotes::install_github("pedrobdfp/eDNA_structure", upgrade = "never")During installation, a large amount of black text will appear — this is the Stan model compiling to C++. It only happens once. Every subsequent call to eDNA_dmm() goes straight to sampling with no compilation output.
Installed automatically:
| Package | Purpose |
|---|---|
rstan (≥ 2.21) |
Bayesian inference via Stan |
ggplot2 (≥ 3.4) |
All visualizations |
dplyr, tidyr |
Data manipulation |
vegan (≥ 2.6) |
NMDS ordination |
posterior (≥ 1.4) |
MCMC diagnostics (ESS, Rhat) |
loo (≥ 2.6) |
Leave-one-out cross-validation |
scales |
Axis formatting |
library(eDNAstructure)
library(dplyr) # for pipe and data manipulation
library(ggplot2) # for plot customization
# Option A — use the built-in example dataset
data <- get_example_data()
# Option B — simulate your own dataset with known ground truth
data <- simulate_eDNA_survey(
n_communities = 4,
n_species = 40,
samples_per_community = 5,
seed = 2026
)
# Inspect raw species composition before fitting
plot_true_compositions(
data$counts,
metadata = data$covariates,
facet_var = "TrueCommunity"
)
# Fit the model with a given number of communities (K)
fit <- eDNA_dmm(
counts = data$counts,
covariates = data$covariates[, c("Depth", "Distance_shore")],
K = 4
)
print(fit)
summary(fit)
# Or select K using LOO cross-validation
loo_result <- eDNA_loo(data$counts, data$covariates[, c("Depth", "Distance_shore")],
K_range = 2:5)
loo_result$plot
# The loo_result object stores all fitted models — no need to refit
# Extract the K=4 model directly:
fit <- loo_result$fits[["K4"]]
# Or using the K value as a number:
K_best <- 4
fit <- loo_result$fits[[paste0("K", K_best)]]
# Confirm what you have:
print(fit)
##You can also plot the results!
# Structure plot — one bar per sample, colored by community membership probability
eDNA_dmm_structure(fit, metadata = data$covariates,
facet_var = "TrueCommunity", sort_var = "Depth")
# NMDS ordination colored by community assignment
eDNA_dmm_nmds(fit)$plot
# Prior vs posterior distributions for covariate effects
eDNA_dmm_beta(fit)$plotFor a complete walkthrough — including step-by-step simulation, data formatting, K selection, all visualization options, parameter recovery, and troubleshooting — see the full tutorial vignette. It is designed to be read start to finish and assumes no prior familiarity with Bayesian mixture models.
The primary input to eDNA_dmm() is a sample × taxon matrix of non-negative integer read counts:
- Rows = samples (stations, replicates, individuals, etc.)
- Columns = taxa or ASVs — taxonomic annotation is not required
- Values = raw integer read counts (do not normalize)
data$counts[1:3, 1:5]
# Sp_1 Sp_2 Sp_3 Sp_4 Sp_5
# STN_001 412 310 121 73 0
# STN_002 389 275 98 61 14
# STN_003 52 41 487 312 208If your data are in long format, convert them first:
library(tidyr)
count_matrix <- long_df |>
pivot_wider(names_from = taxon, values_from = reads, values_fill = 0) |>
tibble::column_to_rownames("sample_id") |>
as.matrix()A sample × covariate data frame in the same row order as the count matrix. Covariates are Z-score standardized internally by default.
head(data$covariates)
# sample_id TrueCommunity Depth Distance_shore
# 1 STN_001 1 82 198
# 2 STN_002 1 79 204
# 3 STN_003 2 11 197The core function. Fits a Dirichlet-Multinomial Mixture model via Stan and returns an edna_dmm_fit object.
fit <- eDNA_dmm(
counts = my_counts, # sample × taxon integer count matrix
covariates = my_covs, # sample × covariate data frame, or NULL
K = 4, # number of latent communities to fit
scale_covariates = TRUE, # Z-score standardize covariates (strongly recommended)
chains = 1, # number of MCMC chains (see note on label switching below)
iter = 4000, # total iterations per chain (including warmup)
warmup = 2000, # warmup iterations to discard
adapt_delta = 0.95, # HMC target acceptance rate; increase to 0.99 if divergences
max_treedepth = 12, # increase to 14–15 if "max treedepth exceeded" warnings
seed = 13, # random seed for reproducibility
conc = 0.5, # Dirichlet prior concentration: < 1 = sparse communities
alpha_shape = 5, # Gamma prior shape for overdispersion parameter alpha
alpha_rate = 2, # Gamma prior rate (prior mean = shape/rate = 2.5)
verbose = TRUE # print sampling progress
)The returned edna_dmm_fit object contains:
| Element | Description |
|---|---|
sample_info |
Data frame: posterior membership probabilities and MAP assignment per sample |
pi_mean |
Matrix [K × S]: posterior mean community compositions |
beta_summary |
Data frame: covariate coefficient summaries with ESS and reliability |
alpha_mean |
Scalar: posterior mean overdispersion |
stan_fit |
Raw rstan::stanfit object for advanced diagnostics |
On single chains: Mixture models suffer from label switching across chains — "Community 1" in chain A may map to "Community 2" in chain B, making multi-chain Rhat diagnostics meaningless. A single long chain sidesteps this. Use within-chain ESS (reported by
summary()) as your convergence criterion.
Produces a STRUCTURE-style plot: one vertical bar per sample, divided into colored segments by posterior community membership probability.
p <- eDNA_dmm_structure(
fit,
metadata = my_metadata, # data frame with additional sample variables
sample_id_col = "sample_id", # column in metadata matching sample IDs in fit
facet_var = "TrueCommunity", # facet panels by this variable (e.g. site, year, depth)
sort_var = "Depth", # sort samples within each panel by this variable
community_colors = NULL, # named hex vector (e.g. c("Community 1" = "#E63946"))
# or NULL for automatic HCL palette
bar_width = 0.9, # bar width (0–1); 1 = no gaps between bars
x_text = FALSE, # show sample ID labels on x-axis?
base_size = 11, # base font size in points
title = NULL, # plot title; NULL = auto-generated
subtitle = NULL, # plot subtitle; NULL = auto-generated
legend_position = "bottom" # "bottom", "right", "left", "top", or "none"
)Returns a ggplot2 object — save with ggsave() or extend with additional ggplot2 layers.
Runs NMDS on community dissimilarities and plots samples colored by their MAP community assignment. Point size reflects assignment certainty: larger points are more confidently assigned to a single community.
result <- eDNA_dmm_nmds(
fit,
k = 2, # NMDS dimensions (2 or 3); increase if stress > 0.2
nmds_axes = c(1, 2), # which two axes to display; e.g. c(1,3) for axes 1 and 3
distance = "bray", # dissimilarity metric passed to vegan::vegdist()
use_edna_index = TRUE, # apply eDNA index transform before computing distances
trymax = 100, # maximum random NMDS starts (more = less risk of local optima)
seed = 42, # random seed for NMDS
show_ellipse = TRUE, # draw 95% confidence ellipse per community?
ellipse_type = "t", # ellipse type: "t" (robust) or "norm" (normal-based)
community_colors = NULL, # named hex vector or NULL for automatic palette
size_range = c(1.5, 5), # point size range: c(min, max) mapped to
# 50% certainty (smallest) → 100% certainty (largest)
alpha = 0.85, # point transparency (0 = invisible, 1 = opaque)
base_size = 13,
title = NULL,
subtitle = NULL,
legend_position = "right"
)
result$plot # ggplot2 object
result$nmds # vegan::metaMDS object (access stress value, species scores, etc.)Overlays the prior and posterior distributions for each softmax regression coefficient. A posterior pulled away from the prior is evidence that the covariate genuinely predicts community membership.
result <- eDNA_dmm_beta(
fit,
layout = "joint", # "joint": communities overlaid per covariate panel
# "separate": one row per community, one column per covariate
covariates_to_plot = NULL, # character vector of covariate names to include, or NULL for all
show_intercept = FALSE, # include the intercept term?
beta_prior_sd = 1.0, # prior SD — must match the Stan model (default: Normal(0,1))
n_prior_samples = 4000, # prior draws for the density curve (more = smoother)
community_colors = NULL, # named hex vector or NULL
prior_color = "grey60", # fill color for the prior density
prior_alpha = 0.35, # prior density transparency
posterior_alpha = 0.55, # posterior density transparency
show_annotations = NULL, # NULL = auto (shown for K=2 only); TRUE or FALSE to override
base_size = 13,
title = NULL,
subtitle = NULL
)
result$plot # ggplot2 object
result$table # data frame: mean, 90% CI, P(direction), ESS, reliability per coefficientFits models across a range of K values and compares them using Leave-One-Out cross-validation. Returns an elbow plot and a comparison table to guide K selection.
loo_result <- eDNA_loo(
counts = my_counts,
covariates = my_covs,
K_range = 2:5, # integer vector of K values to evaluate
scale_covariates = TRUE,
chains = 1,
iter = 4000,
warmup = 2000,
adapt_delta = 0.95,
seed = 13,
conc = 0.5,
alpha_shape = 5,
alpha_rate = 2,
verbose = TRUE
)
loo_result$plot # LOO-ELPD elbow plot (higher = better; look for the elbow)
loo_result$loo_table # data frame: K, LOO-ELPD, SE
loo_result$loo_compare # loo::loo_compare() output
loo_result$fits # named list of edna_dmm_fit objects, one per KVisualizes the observed species frequencies per sample as stacked bars — the same layout as eDNA_dmm_structure(), allowing direct before/after comparison. Most useful before fitting to inspect the raw community signal, and with simulated data where true community labels are known.
p <- plot_true_compositions(
counts,
metadata = my_metadata, # data frame for faceting and sorting
sample_id_col = "sample_id",
facet_var = "TrueCommunity", # facet by known or hypothesized grouping
sort_var = "Depth",
top_n = 20, # show top N taxa individually; rest collapsed to "Other"
bar_width = 0.95,
base_size = 11,
title = NULL,
subtitle = NULL,
legend_position = "none" # default none — too many taxa for a useful legend
)Returns a ggplot2 object.
Visualizes the posterior mean taxonomic composition of each latent community as stacked bars — the model's estimate of what each community "looks like" in species space. Colors match those used in eDNA_dmm_structure() and plot_true_compositions() for direct comparison.
rp <- eDNA_dmm_compositions(
fit,
top_n = 20, # show top N taxa; rest collapsed to "Other"
base_size = 13,
title = NULL,
subtitle = NULL,
legend_position = "right", # "right", "bottom", "left", "top", or "none"
bar_width = 0.7 # bar width (0–1)
)Returns a ggplot2 object. The x-axis labels show community numbers (1, 2, 3…). Pair with eDNA_dmm_structure() to connect community identities to sample assignments.
Returns the built-in simulated dataset: 20 samples × 40 taxa across 4 communities separated by depth and distance from shore. Generated by simulate_eDNA_survey() with known ground truth, so fitted parameters can be compared to the true values.
data <- get_example_data()
# data$counts — integer matrix [20 × 40]
# data$covariates — data frame: sample_id, TrueCommunity, Depth, Distance_shore
# data$community_compositions — true composition matrix [4 × 40]
# data$metab_df — raw simulated metabarcoding reads
# data$sample_metadata — full simulation metadataeDNAstructure includes a mechanistic simulation pipeline for generating eDNA datasets with known community structure. Use it for method validation, power analysis, or teaching. The full tutorial (vignettes/tutorial.Rmd) walks through the pipeline in detail.
# Full pipeline in one call
sim <- simulate_eDNA_survey(
n_communities = 4, # number of communities K
n_species = 40, # number of species S
samples_per_community = 5, # sampling stations per community
community_covariate_means = NULL, # K × P matrix of covariate means per community
# NULL = default 2-covariate depth × shore design
covariate_sds = NULL, # K × P SDs (NULL = 5 for all)
mean_read_depth = 10000, # mean reads per sample
bio_reps = 1, # biological replicates per station
seq_reps = 1, # sequencing technical replicates per bio rep
spillover = 0.15, # fraction of species frequency leaking between communities
shedding_error = 0.3, # lognormal SD on per-organism eDNA shedding
decay_rate = 0.1, # exponential distance decay of eDNA signal
seed = 42
)
# sim$counts — ready for eDNA_dmm()
# sim$covariates — ready for eDNA_dmm()Or run each step individually for full control:
community_mat <- generate_community_compositions(
n_communities = 4,
n_species = 40,
n_dominant_range = c(2, 3), # 2–3 high-frequency dominant species per community
n_unique_low_freq = 4, # species present only in one community
community_groups = c(1,2,1,2), # group structure for shared species
n_group_shared = 4, # species shared within each group
spillover = 0.15, # cross-community frequency leakage
seed = 1
)
contrib_obj <- generate_contributors(
community_compositions = community_mat,
samples_per_community = 5,
n_contributors_range = c(20, 100), # organisms per sample
n_distribution = "Negative Binomial", # overdispersed organism counts
size_range = c(1, 10), # body size range (affects eDNA shedding)
distance_range = c(0, 100), # distance from sampler (affects eDNA decay)
seed = 1
)
eDNA_obj <- generate_eDNA(
contributors_list = contrib_obj$contributors_list,
shedding_rate = 1000, # baseline molecules shed per unit size
beta = 0.75, # allometric exponent (metabolic scaling)
decay_rate = 0.1, # exponential distance decay
shedding_error = 0.3, # lognormal noise on shedding (SD on log scale)
bottle_volume = 0.1, # fraction of local eDNA pool captured per bottle
bio_reps = 1
)
metab_df <- simulate_metabarcoding(
eDNA_obj,
mean_read_depth = 10000, # mean total reads per sample
read_depth_sd = 0.2, # lognormal SD for read depth variation
error_sd = 0.05, # Gaussian noise added to species frequencies
rep = 1 # sequencing technical replicates per bottle
)
sample_metadata <- generate_sample_covariates(
contributors_list = contrib_obj$contributors_list,
community_covariates = matrix( # community-specific covariate means
c(80, 200, 10, 200, 80, 20, 10, 20),
nrow = 4, byrow = TRUE,
dimnames = list(NULL, c("Depth", "Distance_shore"))
)
)| Function | Purpose |
|---|---|
simulate_eDNA_survey() |
Full pipeline in one call |
generate_community_compositions() |
K community frequency vectors over S species |
generate_contributors() |
Organisms shedding eDNA per sample |
generate_eDNA() |
Shedding, exponential decay, bottle sub-sampling |
simulate_metabarcoding() |
Amplification bias and multinomial read counts |
generate_sample_covariates() |
Environmental metadata drawn from community-specific distributions |
For sample i, the DMM marginalizes over a latent community assignment z_i:
- Compositions: π_k ~ Dirichlet(conc · 1_S) for k = 1…K
- Membership: P(z_i = k) = softmax(β_0k + β_1k · x_1i + … + β_Pk · x_Pi), community K = reference
- Counts: x_i | z_i = k ~ DirichletMultinomial(N_i, α · π_k)
The global overdispersion α absorbs both technical (PCR, sequencing) and ecological compositional variance. Marginalizing over z_i makes inference exact.
Why only one chain?
Label switching: "Community 1" in chain A may map to "Community 2" in chain B. Multi-chain Rhat values are pathological even when each chain converges perfectly. One long chain avoids this. Check within-chain ESS instead (printed by summary()).
I have divergent transitions. What do I do?
Increase adapt_delta toward 0.99. If they persist, try lower K or verify your count matrix has no all-zero rows or columns.
Can I use raw ASVs instead of taxonomy-collapsed counts? Yes. The model treats each column as a compositional unit and does not use taxonomy. ASVs give finer resolution; taxa collapse dimensionality and often converge faster.
How do I include year as a covariate? Pass it as a numeric column. But if you have only a few discrete years, the linearity assumption may be too strong — consider fitting without year and testing it post-hoc via multinomial regression on the posterior assignments.
The first run takes forever — is something wrong? No. Stan compiles the model to C++ on the first call after installation (1–2 minutes). All subsequent calls skip compilation. This is normal behavior for any rstan-based package.
If you use this package in published research, please cite:
Brandão-Dias et al. (year). Multinomial mixture models from environmental DNA reveal depth stability and dynamic surface turnover of marine vertebrate communities. Under review.
Brandão-Dias et al. (year). eDNAstructure: Dirichlet-Multinomial Mixture Models for eDNA Metabarcoding Community Structure. R package version 0.1.0. https://github.com/pedrobdfp/eDNA_structure
Please also cite Stan:
Carpenter B. et al. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1).
MIT © eDNAstructure authors