Geoadditive Regression Model for CoDa

• The paper introduces a geoadditive regression model for CoDa that respects simplex constraints while integrating nonlinear and spatial effects.
• It employs an isometric log-ratio transformation with a multivariate Gaussian additive framework for coherent compositional predictions and uncertainty quantification.
• Robust simulation studies and soil texture applications validate the model using model comparison metrics like WAIC and CoDa‑R².

Searching arXiv for the cited CoDa geoadditive and related Bayesian compositional regression papers. Use the arXiv search tool to look up:

1. (Martínez-Minaya et al., 11 Aug 2025)
2. (Martínez-Minaya et al., 2023)
3. (Maier et al., 16 Oct 2025)
4. geoadditive compositional data logistic-normal Bayesian spatial regression
5. Bayes Hilbert space additive density regression arXiv A geoadditive regression model for compositional data (CoDa) is a regression model for a compositional response that respects simplex constraints, handles nonlinear covariate effects, captures residual spatial structure, yields coherent predictions summing to one, and supports full Bayesian uncertainty quantification. In the formulation developed for soil texture prediction, a composition on the simplex is represented with orthonormal log-ratio coordinates, a flexible multivariate Gaussian additive model is fitted to those coordinates with nonlinear and spatial terms, and the fitted values are mapped back to the simplex for interpretation (Martínez-Minaya et al., 11 Aug 2025).

1. Definition and motivation

Compositional data are vectors of positive parts carrying only relative information and constrained to sum to a constant. Soil texture is the motivating example: sand, silt, and clay percentages sum to $100\%$. Because of the unit-sum constraint, the parts are not free to vary independently; if one part increases, at least one other must decrease. Standard regression on raw parts is therefore inappropriate because it ignores positivity, constant-sum dependence, and the fact that ratios, not absolute values, carry the meaningful information. Analyzing components separately can even produce impossible predictions, such as predicted percentages summing to more than $100\%$ (Martínez-Minaya et al., 11 Aug 2025).

The geoadditive extension is motivated by settings in which compositional responses vary in space and respond to environmental drivers nonlinearly. Soil texture is a canonical case: compositions depend on terrain, lithology, slope, elevation, and latent spatial structure. A purely linear CoDa model is too rigid. A geoadditive model allows nonlinear effects of covariates like elevation and slope, a smooth spatial field over longitude–latitude, and simultaneous modeling of all composition parts coherently. This is particularly useful for digital soil mapping, where residual spatial pattern remains even after including covariates (Martínez-Minaya et al., 11 Aug 2025).

The resulting model is both methodological and practical. It shows how to carry out Bayesian geoadditive CoDa regression with brms in R and how to evaluate fit using CoDa-specific Bayesian $R^2$ measures, namely BR-CoDa-$R^2$ and BM-CoDa-$R^2$, together with WAIC for model selection and evaluation (Martínez-Minaya et al., 11 Aug 2025).

2. Simplex geometry and log-ratio coordinates

The sample space for a $D$-part composition is not $\mathbb{R}^D$ but the simplex,

$$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$

with $\kappa>0$, and the model sets $\kappa=1$ without loss of generality. The paper adopts Aitchison geometry, which gives the simplex Euclidean structure through appropriate operations and inner product. In this geometry, perturbation and powering are

$100\%$0

and the Aitchison inner product is

$100\%$1

where $100\%$2 is the geometric mean. This geometry justifies working with orthonormal coordinates and makes distance, variance, and regression compatible with relative information (Martínez-Minaya et al., 11 Aug 2025).

The key coordinate map is the isometric log-ratio transformation. Given an orthonormal basis $100\%$3 of $100\%$4,

$100\%$5

The paper uses the Gram–Schmidt ilr basis. Because ilr is an isometric isomorphism from the simplex to $100\%$6, ordinary multivariate Gaussian regression becomes valid after transformation. For the soil-texture case with $100\%$7, the two ilr coordinates are interpreted as the balance of sand versus silt and the balance of the geometric mean of sand and silt versus clay (Martínez-Minaya et al., 11 Aug 2025).

This coordinate representation is central to the model’s coherence. The transformed response retains the compositional geometry while moving estimation into a Euclidean space where multivariate additive Gaussian models are available. The inverse map $100\%$8 then returns fitted values and posterior predictions to the simplex, ensuring valid compositions (Martínez-Minaya et al., 11 Aug 2025).

3. Statistical formulation of the geoadditive Bayesian model

The model adopts the logistic-normal approach: $100\%$9 is logistic-normal if and only if $R^2$0 is multivariate normal. For observation $R^2$1,

$R^2$2

The covariance matrix is shared across observations, with diagonal terms $R^2$3 and off-diagonal terms $R^2$4. The generic Bayesian additive CoDa model is

$R^2$5

Thus each ilr coordinate receives its own additive predictor, and cross-coordinate dependence is captured by the residual covariance matrix $R^2$6 (Martínez-Minaya et al., 11 Aug 2025).

In the Basque Country soil application, the fitted specification is

$R^2$7

Here $R^2$8 is a fixed effect for lithology, $R^2$9 is a random effect, $R^2$0 and $R^2$1 are univariate penalized P-splines, and $R^2$2 is a bivariate tensor-product P-spline over spatial coordinates. The spatial component is therefore not a GP or CAR in this paper; it is a 2D tensor-product smooth in longitude and latitude, implemented separately for each ilr coordinate (Martínez-Minaya et al., 11 Aug 2025).

Implementation is done in R with brms, using its multivariate-response Bayesian GAM framework backed by Stan. The workflow is conceptually simple: transform the composition to ilr coordinates, fit a multivariate Gaussian additive model to those coordinates, and posterior-sample fitted values and back-transform with $R^2$3. Because the response in transformed space is Gaussian, the likelihood is multivariate Gaussian and the link is effectively identity, since $R^2$4. Posterior inference is via Hamiltonian Monte Carlo as implemented in Stan. The priors are a zero-centered Gaussian prior with standard deviation $R^2$5 for fixed effects, Gaussian priors with mean $R^2$6 and standard deviation $R^2$7 for smooth coefficients, a half Student-$R^2$8 prior with $R^2$9 degrees of freedom for standard deviations, and an LKJ prior with regularization parameter set to $R^2$0 for the correlation matrix (Martínez-Minaya et al., 11 Aug 2025).

4. Goodness-of-fit and model comparison in the compositional setting

A principal methodological contribution is the extension of Bayesian $R^2$1 to compositional outcomes through ilr variance decomposition. The paper defines total compositional variance as

$R^2$2

This decomposition is geometrically coherent because it uses the sum of ilr-coordinate variances, which is justified by the isometry of ilr (Martínez-Minaya et al., 11 Aug 2025).

Explained variance in CoDa is defined as the variance of predicted ilr coordinates,

$R^2$3

where $R^2$4. For posterior draw $R^2$5,

$R^2$6

Two residual definitions are proposed. In BR-CoDa-$R^2$7, the residual variance is

$R^2$8

whereas in BM-CoDa-$R^2$9, it is

$D$0

Values near $D$1 indicate that most ilr variance is explained; values near $D$2 indicate little explained variance. The paper notes that BR-CoDa-$D$3 tends to have a narrower posterior than BM-CoDa-$D$4, paralleling remarks in univariate Bayesian $D$5 (Martínez-Minaya et al., 11 Aug 2025).

The paper also uses WAIC,

$D$6

with lower WAIC indicating better expected predictive fit penalized for complexity. WAIC and CoDa-$D$7 are explicitly treated as complementary: WAIC for out-of-sample predictive adequacy, and BR/BM-CoDa-$D$8 for interpretable explained compositional variance. The authors also propose probabilistic model comparison using posterior probabilities such as

$D$9

with $\mathbb{R}^D$0, so that probabilities $\mathbb{R}^D$1 indicate substantial superiority of $\mathbb{R}^D$2, probabilities $\mathbb{R}^D$3 favor $\mathbb{R}^D$4, and intermediate values indicate similar explanatory power (Martínez-Minaya et al., 11 Aug 2025).

5. Simulation evidence and soil texture prediction

The methodology is validated through three simulation studies. In a linear CoDa regression on compositions in $\mathbb{R}^D$5 with $\mathbb{R}^D$6 and covariates $\mathbb{R}^D$7, posterior means and $\mathbb{R}^D$8 credible intervals recovered the true coefficients, standard deviations, and correlations well. In a CoDa GAM on compositions in $\mathbb{R}^D$9 with $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$0, the model recovered both the univariate and bivariate nonlinear effects and residual covariance well. In a model-comparison simulation, the true linear model $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$1 achieved WAIC $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$2, BR-CoDa-$$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$3, and BM-CoDa-$$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$4, whereas a model with only irrelevant covariate $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$5 had WAIC $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$6 and BR/BM-CoDa-$$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$7. In the GAM case, the best model with both smooths and richer basis achieved WAIC $$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$8, BR-CoDa-$$\mathcal{S}^D = \left\{ \boldsymbol{y} = \left[ y_1, \ldots, y_D \right]; \ y_d > 0,\ d = 1,\ldots,D;\ \sum_{d=1}^D y_d = \kappa \right\},$$9, and BM-CoDa-$\kappa>0$0. These results show that the geoadditive CoDa model can recover smooth effects and distinguish correct from misspecified models (Martínez-Minaya et al., 11 Aug 2025).

The empirical application uses $\kappa>0$1 soil samples at $\kappa>0$2 cm depth in the Basque Country from 2010 to 2018, with a three-part response composition of sand, silt, and clay and covariates at spatial resolution $\kappa>0$3. Candidate models use Lithology, Elevation, Slope, a Longitude–Latitude spatial smooth, and a Year random effect. The best model, $\kappa>0$4, includes all five components and attains WAIC $\kappa>0$5, BR-CoDa-$\kappa>0$6, and BM-CoDa-$\kappa>0$7. The second-best model, $\kappa>0$8, omits Year and has WAIC $\kappa>0$9 and BR/BM-CoDa-$\kappa=1$0. The paper states that Year contributed about a $\kappa=1$1 increase in explained compositional variance. Models omitting spatial and environmental terms performed clearly worse; for example, the lithology-only model $\kappa=1$2 had BR/BM-CoDa-$\kappa=1$3 (Martínez-Minaya et al., 11 Aug 2025).

The estimated effects are interpreted on the log-ratio scale and then translated back to compositions. With Surface deposits as the lithology reference level, category $\kappa=1$4 (Decarbonate marls) had the largest positive effect for ilr1, corresponding to sand versus silt, and category $\kappa=1$5 (Limestones and dolomites) had the strongest positive effect for ilr2, corresponding to the $\kappa=1$6-to-clay balance. The paper interprets these as about $\kappa=1$7 higher sand-to-silt ratio and about $\kappa=1$8 higher $\kappa=1$9-to-clay balance relative to the reference. Increasing slope is associated with higher clay and lower sand and silt. Higher elevations correspond to a slight increase in sand and clay and a decrease in silt. The 2D longitude–latitude smooth revealed clear spatial patterns: sand higher in specific areas, silt higher where sand decreases, and clay more homogeneous with weaker spatial gradient. Uncertainty was lower in densely sampled areas and higher in peripheral or sparse regions. The fitted model was then used to predict posterior mean maps for sand, silt, and clay, posterior standard deviation maps, and USDA textural class maps (Martínez-Minaya et al., 11 Aug 2025).

6. Relation to other Bayesian CoDa frameworks and methodological scope

Within Bayesian CoDa regression, the geoadditive ilr-Gaussian model sits alongside Dirichlet and logistic-normal formulations. The paper explicitly states that, compared with Dirichlet regression, this approach is preferable when richer covariance structure is needed, when one wants additive smooths and spatial terms in a standard multivariate Gaussian framework, and when Bayesian multivariate GAM software is desired. Compared with generic logistic-normal models, the specific contribution is the geoadditive extension plus CoDa-specific $100\%$00 measures (Martínez-Minaya et al., 11 Aug 2025).

A closely related development is the logistic-normal Dirichlet model (LNDM), which proposes a general Bayesian latent Gaussian framework for compositional responses in R-INLA. In that framework, the response is transformed with the additive log-ratio transform, the transformed coordinates are modeled through additive predictors of the form

$100\%$01

and spatial, temporal, and other structured additive effects can be inserted naturally within latent Gaussian models. Under the usual interpretation of geoadditive models as additive regressions with structured spatial and possibly other smooth or random terms, this means the method can be used as a geoadditive or spatially structured additive regression model for CoDa, although it is not presented as a single named geoadditive model class in the classical semiparametric sense (Martínez-Minaya et al., 2023).

A further extension arises in additive density regression in Bayes Hilbert space. That framework treats densities as positive objects carrying only relative information under a normalization constraint and formulates a structured additive density-on-scalar regression model

$100\%$02

It does not fit an explicit spatial effect in the reported application, but it is formulated as a general structured additive model with tensor-product bases and penalties. This suggests a functional or infinite-dimensional analogue of geoadditive CoDa regression in which the response is a density-valued composition rather than a finite-dimensional simplex vector (Maier et al., 16 Oct 2025).

The scope of the geoadditive CoDa model is accompanied by clear assumptions and limitations. The formulation requires positive compositions and does not develop a zero-handling method. Regression is performed in ilr coordinates, whose numerical form depends on the chosen orthonormal basis, although the methods are valid under any ilr and interpretability varies by basis. The model assumes that the ilr-transformed response follows a multivariate Gaussian distribution and that a common residual covariance matrix $100\%$03 is adequate across observations. The paper also does not provide a special identifiability derivation beyond standard GAM regularization and centering of smooths, and computational details such as MCMC tuning are not deeply discussed in the supplied text (Martínez-Minaya et al., 11 Aug 2025).