Multimetab: Bayesian Metabolomics Regression

• Multimetab is an R package that implements a Bayesian variable selection framework for metabolite responses, directly modeling censoring, skewness, and correlated predictors.
• It combines a skew‐normal censored mixture model with spike‐and‐slab regression and a Markov random field prior to leverage nutritional similarities among dietary predictors.
• Posterior inference via MCMC yields inclusion probabilities, coefficient estimates, and model diagnostics that demonstrate improved power and fit in metabolomics cohort studies.

Searching arXiv for the primary paper and closely related work to ground the article in current literature. Multimetab is an R package that implements a Bayesian variable selection framework for metabolite response variables. It was introduced for identifying food items associated with blood metabolites in two cohorts of healthcare professionals, and it is explicitly designed for high-dimensional regression with strongly correlated predictors and non-normal metabolite outcomes measured by mass spectrometry. Its core methodological contribution is the combination of a skew-normal censored mixture model with a Markov random field prior that exploits substantive information on nutritional similarities among dietary factors, with posterior inference delivered through MCMC and summarized through posterior inclusion probabilities, coefficient estimates, and model-fit diagnostics (Clark-Boucher et al., 8 Sep 2025).

1. Problem setting and analytical niche

The motivating setting for multimetab consists of $30$ food intake variables used as predictors and $244$ compounds measured by mass spectrometry used as metabolic responses. The modeling problem is complicated by missingness, left-censoring, and skewness in the metabolite outcomes, together with strong correlation among dietary predictors. Existing methods are described as capable of addressing such factors in low-dimensional settings, but not as designed for high-dimensional regression involving strongly correlated predictors and non-normal outcomes (Clark-Boucher et al., 8 Sep 2025).

Within this formulation, the predictors of interest are dietary variables, while additional covariates enter as confounders. The stated goal is not generic association screening, but variable selection under joint accommodation of censoring, skewness, and prior structural information about the predictors. This places multimetab in a distinct methodological niche within metabolomics: it is not a workflow manager, a database, or a general multi-omics container, but a specialized inferential package for regression with metabolite response variables.

A common simplification is to treat below-detection values as an imputation problem and right-skewness as a transformation problem. Multimetab is structured against that simplification. Its defining premise is that censoring, point mass values, and skew-normal behavior should be modeled directly rather than removed by preprocessing alone (Clark-Boucher et al., 8 Sep 2025).

2. Response model and likelihood construction

For subject $i$, the framework defines $Y_i$ as the observed metabolite value and $Y_i^*$ as the true underlying value, observed only if it exceeds a detection limit $\psi$. The latent continuous metabolite level is represented as

$$V_i = \beta_0 + \sum_{j=1}^p \beta_j X_{ij} + \sum_{t=1}^s \alpha_t C_{it} + \epsilon_i,$$

with

$Y_i^* = U_i V_i,$

and

$$Y_i = \begin{cases} Y_i^*, & Y_i^* \geq \psi \ \text{PMV (missing)}, & Y_i^* < \psi . \end{cases}$$

Here, $U_i=0$ represents biological zeros, while $244$0 with $244$1 represents technical zeros arising from censoring (Clark-Boucher et al., 8 Sep 2025).

The error term is modeled as skew-normal,

$244$2

where $244$3 is a scale parameter and $244$4 is a skewness parameter learned for each metabolite. This construction is intended to account directly for heavy right skewness and outliers prevalent in mass spectrometry metabolite data. The paper contrasts this with ad hoc log transformation or imputation (Clark-Boucher et al., 8 Sep 2025).

For each observation, the likelihood is given as

$244$5

where $244$6 if $244$7 is observed and $244$8 if it is censored, and

$244$9

This likelihood combines explicit censoring with skew-normal response behavior in a single model rather than treating them as separate preprocessing stages.

3. Variable selection and structured prior information

Variable selection in multimetab is based on a spike-and-slab prior. For dietary predictor $i$0,

$i$1

where the latent binary indicator $i$2 determines whether predictor $i$3 is included. Inference is summarized through posterior inclusion probabilities, which function as the primary variable-selection output (Clark-Boucher et al., 8 Sep 2025).

The central extension beyond standard spike-and-slab regression is a Markov random field prior on the selection indicators,

$i$4

Here, $i$5 is a pairwise relationship matrix encoding prior knowledge such as nutritional similarity among food variables, $i$6 controls overall sparsity, and $i$7 controls the influence of the relationship structure. When $i$8, the formulation reduces to an independent Bernoulli prior (Clark-Boucher et al., 8 Sep 2025).

This prior is designed to encourage joint selection of related predictors. In settings where predictors are nutritionally similar or highly correlated in intake, the model attempts to reduce signal dilution by allowing prior dependence across selection indicators. The paper further introduces a prior-based method for hyperparameter specification of $i$9, intended to avoid computationally intensive data-driven grid search and hyperprior specifications that can create identifiability or instability problems (Clark-Boucher et al., 8 Sep 2025).

4. Computation, software implementation, and outputs

Posterior computation in multimetab uses Gibbs sampling with data augmentation for latent $Y_i$0, $Y_i$1, and $Y_i$2. The spike-and-slab structure follows Kuo and Mallick (1998) for Gibbs updates of the coefficient parameters, and the skew-normal errors are handled through latent-variable representations following Sahu et al. (2003). Selection indicators $Y_i$3 are updated from Bernoulli full conditionals determined by $Y_i$4, $Y_i$5, and $Y_i$6 (Clark-Boucher et al., 8 Sep 2025).

The package is available at https://github.com/dclarkboucher/multimetab. Its described features include the full skew-normal censored mixture model, optional MRF or independent priors on selection indicators, automatic construction of the diet relationship matrix $Y_i$7 from user-supplied hierarchical groupings, data preprocessing tools for censored or missing metabolite values, and MCMC posterior sampling of all model parameters. Reported outputs include posterior inclusion probabilities, posterior estimates of regression coefficients, and model-fit diagnostics such as posterior predictive checks and ELPD/WAIC (Clark-Boucher et al., 8 Sep 2025).

Two practical constraints delimit the package’s scope. First, it is designed for the analysis of one metabolite at a time as a response. Second, although guidance is given for weakly informative defaults, the user is advised to set the slab variance $Y_i$8 using empirical variance from marginal regressions, estimate the detection limit $Y_i$9 as the minimum observed nonzero value for each metabolite, and perform thinning and convergence checks for MCMC chains.

Component	Specification	Role
Response model	Skew-normal censored mixture model	Handles point mass values, censoring, and skewness
Variable selection	Spike-and-slab with $Y_i^*$0	Produces posterior inclusion probabilities
Structured prior	MRF prior with $Y_i^*$1	Encourages joint selection of related predictors
Computation	Gibbs sampling with data augmentation	Posterior inference
Output	PIP, $Y_i^*$2, posterior predictive checks, ELPD/WAIC	Interpretation and model assessment

5. Empirical application and reported performance

The empirical application uses two cohorts: the Mind Body Study from NHSII, comprising women, and the Men’s Lifestyle Validation Study from HPFS, comprising men. The reported sample sizes are $Y_i^*$3 for MBS and $Y_i^*$4 for MLVS. The analysis includes $Y_i^*$5 metabolites with less than $Y_i^*$6 point mass values and $Y_i^*$7 foods or food groups, with detailed nutritional groupings used to define the relationship matrix $Y_i^*$8 (Clark-Boucher et al., 8 Sep 2025).

In these data, the framework identified $Y_i^*$9 unique metabolites associated with at least one food variable and $\psi$0 total food–metabolite associations, with $\psi$1 foods having at least one association. The largest number of signals was reported for fish, nuts, coffee, and fermented dairy. Both positive and negative associations were identified, whereas nuts, red meat, and carotenoid vegetables were reported as having only negative associations (Clark-Boucher et al., 8 Sep 2025).

The role of the MRF prior is emphasized in both simulations and cohort analyses. MRF-based analyses detected additional biologically plausible metabolite–diet links not found by models using independent priors, especially involving dairy products and triacylglycerols, and no associations were found by independent priors only. The simulation studies and real-data analysis are described as showing gain in power without increase in FDR. Model fit also favored the skew-normal specification over normal error models, both by visual posterior predictive checks and by higher ELPD/WAIC; the MRF prior provided additional, though smaller, improvements in fit (Clark-Boucher et al., 8 Sep 2025).

The reported biological interpretation is deliberately conservative. Associations involving coffee and caffeine metabolites, fish and phospholipids with fish-derived fatty acids, and dairy intake and triacylglycerols containing odd-chain fatty acids are presented as consistent with previous research, while other findings are described as potentially novel.

6. Position within metabolomics infrastructure and terminological scope

Within the broader metabolomics software and data ecosystem, multimetab occupies a narrower role than longitudinal workflow platforms or literature-derived summary databases. MeTime, for example, is an open-source R package for reproducible longitudinal metabolomics data analysis built around a metime_analyser S4 container, modular add_*, mod_*, calc_*, meta_*, get_*, and write_* functions, and automatic HTML/PDF reporting; it is architected for multiple datasets and extensibility to multivariate or multimodal analyses, including what its description terms “multimetab” analysis (Marella et al., 8 May 2026). By contrast, multimetab is a targeted Bayesian regression package centered on metabolite response variables and structured variable selection (Clark-Boucher et al., 8 Sep 2025).

A second distinction concerns the difference between within-study regression and across-study synthesis. The Aging Metabolite Database compiles aging-related information for almost $\psi$2 non-unique metabolites from $\psi$3 publications and combines subsets of human plasma studies to produce summary statistics and meta-analytic outputs for common metabolites (Milner et al., 2024). Multimetab does not perform this kind of literature database aggregation; its inferential unit is the cohort-level association analysis between dietary predictors and a single metabolite response.

A third distinction concerns data representation. Work on numerical representations of metabolic systems argues that raw instrument output should not be conflated with scientifically interpretable data, and distinguishes multiple levels ranging from raw numbers to concentration data and tailor-made biological activity data (Smilde et al., 2019). Multimetab presupposes that metabolite measurements can be modeled statistically as censored and skewed responses, rather than treating them as uncomplicated ratio-scaled concentrations by default. This suggests that the package is especially aligned with settings in which the measurement process and detection limits remain active determinants of inference.

The term “multimetab” is also not semantically uniform across the literature. In the supplied sources it can refer to the present R package (Clark-Boucher et al., 8 Sep 2025), to multivariate or multimodal longitudinal omics analysis (Marella et al., 8 May 2026), and, in a different modeling tradition, to multi-species ecosystem metabolism in COMETS summaries (Dukovski et al., 2020). In current usage, therefore, “Multimetab” is best understood contextually rather than as a single stabilized designation across metabolomics and systems biology.