| 12 Oct 2026
Technical note: MixSIAPy – Bayesian stable isotope mixing models in Python with multi-backend inference and formula-based compositional regression
Abstract. Bayesian stable isotope mixing models quantify source contributions, but environmental studies increasingly require testing how multiple interacting drivers alter entire source compositions. The widely used R package MixSIAR restricts covariate structures, while post hoc regression of estimated proportions neither propagates posterior uncertainty nor preserves compositional dependence. Its reliance on a separately installed JAGS engine also constrains direct evaluation of gradient based samplers and accelerator-aware computation. We here introduce MixSIAPy, an open-source Python framework addressing these statistical and computational limitations while retaining established MixSIAR model structures. Its principal methodological contribution is formula-based compositional regression embedded directly in the mixing-model likelihood. Multiple categorical and numerical predictors and their interactions act jointly on source compositions in isometric log-ratio coordinates, with regression, grouping, source, and error parameters estimated in one posterior. A single PyMC model can be fitted through four interchangeable No-U-Turn Sampler backends, enabling CPU sampling and optional GPU execution. A shared engine supports scripted, local graphical, and hosted web workflows with reproducible outputs. Across 20 paired benchmarks, 19 converged in both implementations and showed close posterior correspondence; MixSIAPy reached an acceptable posterior faster in 13 cases. Backend performance was model dependent, and GPUs became advantageous only for some models at high chain parallelism. In an additional crop-water-uptake case study, the model identified a biomass associated shift from shallow to middle soil water in linseed, whereas broader effects remained uncertain. MixSIAPy thus enables interacting environmental controls to be tested within source-partitioning models instead of through separate analyses of uncertain estimated proportions.
Competing interests: At least one of the (co-)authors is a member of the editorial board of Hydrology and Earth System Sciences.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.
Summary of the manuscript
The manuscript presents MixSIAPy, a Python package for Bayesian stable isotope mixing models with multiple inference backends and formula-based compositional regression. Comparisons with MixSIAR and a crop water-source application demonstrate its use. The formula interface allows source proportions to be examined in relation to multiple predictors and their interactions. The benchmarks also illustrate how computational performance varies across model configurations and inference backends. The software is a useful contribution. I recommend minor revision to clarify the validation scope, interpretation of results and application details.
Major comments
(1) Please distinguish the features benchmarked against MixSIAR from the additional interaction functionality illustrated in the crop example. Cite any existing checks of the latter and clarify the scope of the validation in Sections 2.2 and 4.1.
(2) The overlap metric in Eq. (4) can equal one even when interval widths differ. Please acknowledge this limitation and avoid equating high overlap with equivalent uncertainty. The standardized mean difference should likewise be described as a descriptive measure, rather than a test against Monte Carlo error.
Minor comments
(1) State in the Figure 3b caption that timing excludes earlier unsuccessful runs, consistent with Section 3.
(2) Reconcile the normal sampling settings: the main text reports 3000 retained draws in total, whereas manual Section 8.2 lists 5000 per chain for three chains. Specify the settings actually used.
(3) Clarify how the 12 datasets reported in Section 3 correspond to the apparent 11 dataset families in Table A1.
(4) Specify that the 13 of 19 timing comparisons refer to model configurations, several of which share the Alligator dataset.
(5) Add cross-references to the manual for contrast coding, intercept handling, the ILR basis and priors, and identify the error model and priors used in the crop example.
(6) Use consistent notation for the baseline term: η_0 in Eq. (2) and β_0 in the accompanying definition.
(7) Explain the zero-discrimination assumption and revise p. 17 so that it does not imply discrimination uncertainty was propagated in this particular fit.
(8) Check the covariance-matrix dimensions in manual Section 9.1: a J-dimensional tracer vector requires a J × J matrix, rather than K × K.
(9) In Figure 4, distinguish predictions at zero random effects from population-averaged predictions. Define the tail-probability estimator and ensure the significance symbols reflect its Monte Carlo resolution; zero sampled tail events do not imply zero probability.
(10) Identify the software release or commit and link directly to the scripts and environment information for reproducing the benchmarks and Figures 2–4.