Using distributional regression to create interpretable organic carbon models from soil color sensors and X-ray fluorescence

Regressão distributiva para criar modelos interpretáveis de carbono orgânico através de sensores de cor do solo e fluorescência de raios-X

ABSTRACT

Soil carbon content is highly variable both in depth and across the landscape. Thus, modeling it is challenging since carbon data is often heteroskedastic, bimodal, and skewed. The generalized additive model for location, scale, and shape (GAMLSS) is capable of handling complex data by fitting not only the mean (location), but also other parameters related to the scale and shape of distributions. This study used GAMLSS to estimate soil organic carbon (SOC) content in 240 soil samples using as covariates soil color and the chemical composition of samples measured by the Nix color sensor and portable X-ray fluorescence (pXRF), respectively, as well as sample depth. For each unit increase of redness (a*), the SOC’s median decreased by 5.3%. Contents of metallic elements found in mafic minerals (e.g., Fe) had a positive effect on SOC content, whereas SOC decreased with the increase of elements related to felsic minerals (e.g., Si and K). The coefficient of variation (CV) increased in redder samples (6.3% per unit increase of a*) and decreased in yellower samples (3.5% per unit increase of b*). Soil mineralogy was capable of representing the skewness (ν) and kurtosis (τ) of SOC’s distribution. Higher Si and K contents in samples (indicative of felsic minerals) were associated with an increase in extreme observations. GAMLSS could quantify the effects of soil color and chemical composition in different parameters of SOC’s distribution.

Index terms:
GAMLSS; proximal sensing; carbon monitoring; Nix color sensor; pXRF

RESUMO

Teores de carbono do solo são altamente variáveis em profundidade e ao longo da paisagem. Assim, é desafiador modelar carbono no solo, pois, comumente, suas distribuições são heteroscedásticas, bimodais e assimétricas. O modelo aditivo generalizado para posição, escala e forma (GAMLSS) é capaz de lidar com distribuições complexas ao estimar não apenas a média (posição), mas também outros parâmetros relacionados à escala e forma da distribuição. Este estudo usou o GAMLSS para estimar os teores de carbono orgânico do solo (SOC) em 240 amostras de solo utilizando como covariáveis a cor e a composição química medidas pelo colorímetro Nix e por fluorescência de raios-X portátil, respectivamente, além da profundidade de amostragem. Para cada unidade de aumento de vermelhidão (a*), a mediana de SOC diminuiu em 5,3%. Teores de elementos abundantes em minerais máficos (e.g., Fe) apresentaram efeito positivo nos teores de SOC, enquanto os teores de SOC diminuíram com o aumento de elementos comuns em minerais félsicos (e.g., Si e K). O coeficiente de variação foi maior em amostras vermelhas (aumentando em 6,3% por unidade de a*) e menor em amostras mais amarelas (diminuindo 3,5% por unidade de b*). A mineralogia do solo foi relevante para representar a assimetria (ν) e a curtose (τ) da distribuição de SOC. Maiores teores de Si e K nas amostras (indicativos de minerais félsicos) foram associados com maior número de valores extremos. O GAMLSS foi capaz de quantificar os efeitos da cor do solo e de sua composição química em diferentes parâmetros da distribuição de SOC.

Termos para indexação:
GAMLSS; sensoriamento próximo; monitoramento de carbono; sensor de cor Nix; pXRF

Introduction

Great effort is put into monitoring soil organic carbon (SOC) (Pouladi et al., 2023). But creating accurate SOC models is challenging given the complex interactions between the pedological, biotic, and climatic factors that govern SOC dynamics. A common strategy is to gather large datasets of environmental covariates related to SOC and rely on machine-learning algorithms (Pouladi et al., 2023). By doing so, good prediction accuracy can be achieved in many cases, but the behavior of these black-box models is difficult to understand, and their limitations are hard to diagnose. Another downside is that the covariates and algorithms used in machine-learning applications vary considerably (Lamichhane et al., 2019), making models hard to reproduce and compare. Ideally, the desired models need to be accurate, simple, interpretable, and should use covariates that are easy to acquire.

Most statistical models are interpretable, allowing modelers to calculate the effects of each covariate on the response variable. The challenge of using such models with soil data is that they focus on modeling the mean of distributions. Other distributional parameters related to scale and shape (e.g., variance, skewness, kurtosis, etc.) are only implicitly modeled, often resulting in violations of the model’s assumptions due to the complexity of soil data. Machine-learning models are a common alternative, since they make no strong assumptions about the data, but these complex models are difficult to explain (Schmidinger & Heuvelink, 2023). However, although yet underused in many fields, there are distributional statistical models in which each parameter of the distribution is explicitly modeled (Kneib et al., 2023); GAMLSS is one such model (Stasinopoulos et al., 2024).

Proximal sensors, on the other hand, have been extensively used in the modeling of soil properties over the last decade (McBratney et al., 2019; Shi et al., 2023; Ribeiro et al., 2024; Yang & Li, 2026). This includes sensors that measure visible and infrared reflectance, X-ray fluorescence, magnetic susceptibility, among others. Proximal sensors are often portable and easy to use and provide quick analysis of samples. Recent studies have introduced the use of the portable color sensor Nix to estimate SOC (Faria et al., 2022; Singh et al., 2025; Stiglitz et al., 2017). The idea is promising because organic matter absorbs radiation in several regions of the electromagnetic spectrum, but most notably in the visible range (350-700 nm), darkening soil color and making color variables good predictors of SOC (Ben-Dor et al., 2018). In addition to color, certain clay minerals such as Fe oxides and phyllosilicates can bind with organic matter, forming complexes and aggregates that promote its stabilization, playing an important role in the persistence of carbon in soils (Wiseman & Püttman, 2006). Determining and quantifying the mineralogy of soils is not trivial, but the total chemical contents of samples are good representations of soil’s mineralogy and parent material, and can be measured by the portable X-ray fluorescence spectrometer (pXRF) (Mancini et al., 2019).

Accurate and simple models built with data that are easy and quick to obtain can facilitate the monitoring of carbon stocks and provide useful insights about the relationship between SOC and soil color, an important morphological property of soils. The goal of this study is to understand how the color and chemical composition of soil samples can explain the occurrence of SOC in soils by using GAMLSS to model the location, scale, and shape of SOC’s distribution. By modeling beyond the mean, we expect to obtain a model capable of estimating SOC content reliably and providing a broader understanding of the effect of color and chemical variables on the distributional parameters of SOC.

Material and Methods

Study area and soil sampling

The study was conducted on the campus of the Federal University of Lavras (UFLA). A total of 240 soil samples were collected in 124 points (Figure 1). Land uses included: pastures, native forest, planted forests (Eucalyptus and Pinus), and coffee plantations. Soil profiles were first described and classified by pedologists. Samples were then collected in different depths depending on the morphology of soil profiles. Sampling was done along with soil descriptions in profiles and trenches ranging from approximately 0.5 to 2 m in depth. When the desired depth could not be reached, samples were collected using an auger. For modeling purposes, samples were grouped into three depth categories: 0-20, 20-60, and 60-100 cm. Intermediate depths represent B horizons of younger soils or transitional horizons of deep, weathered soils. Sampled profiles were classified according to the Brazilian Soil Classification System (SiBCS) (Santos et al., 2025) and the Soil Taxonomy (Soil Survey Staff, 2014) (Table 1). Studied soils include four soil orders, namely Entisols, Inceptisols, Ultisols, and Oxisols, formed from gabbro and gneiss (Mancini et al., 2019). The distribution of texture classes is presented in Figure 2. Climate in the region is Cwa per Köppen’s classification system (Alvares et al., 2013), which corresponds to humid subtropical climate with hot summers and dry winters. Mean annual precipitation is 1,460 mm and mean annual temperature is 20.4°C.

Figure 1:
Location of sampled points in each depth across the study area located in Lavras, Minas Gerais, Brazil. Maps projected with UTM zone 23K, using the WGS 84 datum. PM: parent material.

Table 1:
Classification of studied soil profiles according to the Brazilian Soil Classification System (SiBCS) (Santos et al., 2025) and the Soil Taxonomy (Soil Survey Staff, 2014), parent material (PM) of each profile (Mancini et al., 2019), and number of samples and profiles within each soil class. Soil orders are presented according to the Soil Taxonomy.

SiBCS (Santos et al., 2025)	Soil Taxonomy (Soil Survey Staff, 2014)	Soil order	PM (Mancini et al., 2019)	Samples	Profiles
CAMBISSOLO HÁPLICO	Typic Dystrustepts	Inceptisols	Gabbro	2	1
CAMBISSOLO HÁPLICO	Typic Dystrustepts	Inceptisols	Gneiss	25	13
GLEISSOLO MELÂNICO	Mollic Endoaquents	Entisols	Gabbro	2	1
GLEISSOLO HÁPLICO	Typic Endoaquents	Entisols	Gabbro	4	2
GLEISSOLO HÁPLICO	Typic Endoaquents	Entisols	Gneiss	5	3
LATOSSOLO AMARELO	Xanthic Hapludoxes	Oxisols	Gneiss	4	3
LATOSSOLO VERMELHO	Anionic Acrudoxes	Oxisols	Gabbro	59	31
LATOSSOLO VERMELHO	Anionic Acrudoxes	Oxisols	Gneiss	16	8
LATOSSOLO VERMELHO-AMARELO	Typic Hapludoxes	Oxisols	Gabbro	9	3
LATOSSOLO VERMELHO-AMARELO	Typic Hapludoxes	Oxisols	Gneiss	34	17
NITOSSOLO VERMELHO	Rhodic Kandiudults	Ultisols	Gabbro	7	4
NITOSSOLO VERMELHO	Rhodic Kandiudults	Ultisols	Gneiss	5	3
ARGISSOLO AMARELO	Inceptic Hapludults	Ultisols	Gneiss	8	3
ARGISSOLO VERMELHO	Typic Rhodudults	Ultisols	Gabbro	11	6
ARGISSOLO VERMELHO	Typic Rhodudults	Ultisols	Gneiss	11	4
ARGISSOLO VERMELHO-AMARELO	Typic Hapludults	Ultisols	Gabbro	8	5
ARGISSOLO VERMELHO-AMARELO	Typic Hapludults	Ultisols	Gneiss	30	17

Figure 2:
Texture distribution and the respective parent material of samples collected across the study area located in Lavras, Minas Gerais, Brazil.

Chemical and colorimetric analyses

Soil samples were air-dried and sieved with a 2-mm mesh. Next, SOC content was determined by oxidizing organic matter with potassium dichromate (K₂Cr₂O₇) in the presence of sulfuric acid (H₂SO₄), followed by titration with ferrous ammonium sulfate (Walkley & Black, 1934). Sample color was measured using the portable color sensor Nix (Stiglitz et al., 2017). Nix is an inexpensive and compact color sensor that uses LED as its light source. The light source is equipped within a small cavity to avoid the interference of external light in measurements. The Nix sensor was configured to use a D50 illuminant and to consider a 2° standard observer. Samples were analyzed in triplicates, being rotated after each repetition. The values of color variables used here are the means of triplicates. Nix provides color measurements in many color spaces, but the color space used in this study was CIELAB, which is defined by variables L*, a*, and b*. Variable L* represents color brightness, with lower values indicating darker colors, and higher values indicating brighter colors. Variable a* represents green (negative values) and red (positive values) colors, but since soils rarely present green hue, in this study a* represents the redness of soil samples. Similarly, variable b* represents blue (negative values) and yellow (positive values) colors, but soil samples usually only have yellow hues, thus here b* is a representation of the yellowness of samples.

The chemical contents of samples were obtained by pXRF. The utilized PXRF was a handheld Bruker S1 Titan LE equipped with a Rh tube (50 keV and 100 µA) capable of emitting X-ray energy directed at samples. This energy displaces electrons from their usual orbit and, when other electrons fill the created gap, energy is emitted in a phenomenon called fluorescence (Silva et al., 2021). The levels of emitted energy are unique to each element and are used to identify which elements are present. The number of energy pulses is used to estimate the content of elements in the sample. Samples were analyzed by pXRF for 60 s each, in triplicates, while rotating the sample every repetition. The equipment was configured to operate in “Soil mode”. The precision of the equipment was checked by scanning the reference samples 2710a and 2711a, certified by the National Institute of Standards and Technology, and the reference samples provided by the manufacturer.

Statistical analyses and modeling

Covariates used in the modeling process

All color variables of the CIELAB color space were used in the modeling process, namely L*, a*, and b*. Depth was added as a categorical variable (0-20, 20-60, and 60-100 cm), since SOC is expected to decrease with depth. The depth group 0-20 cm was used as the reference category. In addition to color and depth, the chemical concentrations of elements obtained via pXRF were also tested. Elements with near-zero variance were removed, leaving 27 detected chemical elements (Ag, Al, As, Bi, Ca, Ce, Cl, Cr, Cu, Fe, K, Mn, Mo, Nb, Ni, P, Pb, Rb, Si, Sr, Ta, Ti, Tl, V, Y, Zn, and Zr). Since there were too many chemical elements, a principal component analysis (PCA) was performed to reduce the dimensionality of pXRF data (Jolliffe & Cadima, 2016). The PCs that explained most of the data’s variance were used in the modeling process. These variables were applied to a stepwise-based algorithm to select the optimal variables for the regression structure of each parameter of the modeled distribution, as explained in the following section.

Statistical modeling with GAMLSS

GAMLSS was used as the distributional regression model of choice. Univariate models such as the generalized linear model (GLM) and the generalized additive model (GAM) assume the response variable follows a distribution belonging to the exponential family, but the variance, skewness and kurtosis of the distribution are not explicitly modeled, but implicitly estimated based on the mean (Stasinopoulos et al., 2024). GAMLSS can explicitly model the parameters for location (e.g., mean, median, etc.), scale (e.g., standard deviation, coefficient of variation, etc.), and shape (e.g., skewness, kurtosis, etc.) of distributions, allowing for the use of a very broad range of distribution families (Equation 1). Mathematically, GAMLSS can be expressed as

${\eta }_k\text{=}{g}_k\left({\theta }_k\right)\text{=}{X}_k{\beta }_k\text{+}{\sum }_{j\text{=}1}^{J_k}p{b}_{jk}\left(x_{jk}\right)$ (1)

where η _k is the predictor, i.e., the regression structure for each parameter; g(‧) are known link functions; k indexes the distribution parameters; β are the coefficients fitted for covariates X; and pb are non-parametric smoothing functions, e.g., a P-spline (Eilers et al., 2015), fitted for J variables.

The chosen distribution was the Box-Cox power exponential distribution (BCPE) (Rigby & Stasinopoulos, 2004), a continuous distribution parameterized by four parameters: median (µ), coefficient of variation (σ), skewness (ν), and kurtosis (τ). Parameters µ, σ and τ use log as the link function, while ν uses the identity. To select the appropriate covariates to fit each parameter, a stepwise-based algorithm was used, called Strategy A. It consists of progressively adding one covariate at a time while modeling µ until the quality of the model increases. After finding the optimal variables for µ, the same is repeated for σ while maintaining the variables chosen for µ. After repeating this process for all parameters, a backward step is conducted, whereby covariates are removed one by one for each parameter to test whether the model would improve. The quality of the model in the stepwise process was measured by the Akaike’s information criterion (AIC), which is a statistical tool used to select competing models built from the same dataset, evaluating the tradeoff between the accuracy and complexity of models. More details about this procedure can be found in Ramires et al. (2021).

The quality of the final model was evaluated by applying diagnostic tools to the normalized quantile residuals (Dunn & Smyth, 1996). Normalized quantile residuals are obtained by calculating the fitted cumulative distribution function (cdf) of residuals and then applying the inverse standard normal cdf to resulting values (Stasinopoulos et al., 2024). Residuals were plotted against fitted values to check for abnormal deviations and were compared against the theoretical quantiles of the normal distribution (QQ-plots) to check for their normality. The worm-plot (Buuren & Fredriks, 2001), which is a detrended QQ-plot (by subtracting the 1:1 line), was used to check for local problems in the adequacy of residuals. The moment bucket plot was used to evaluate the adequacy of the modeled skewness and kurtosis (Bastiani et al., 2022). The bucket plot can help to visually diagnose whether excess kurtosis and skewness are zero per the Jarque-Bera test (Jarque & Bera, 1987), since residuals are expected to be normally distributed.

RESULTS AND DISCUSSION

Principal component analysis (PCA) and variable selection

The two first PCs of the PCA applied to the chemical contents of samples accounted for 42.6% of the variance of the data (Figure 3). The variance explained by other PCs was 10% or less. The two first PCs were therefore selected to be included in the modeling process. PC1 was correlated with metallic elements such as Fe, Cu, Pb, Ti, and others. Most samples that received positive PC1 scores were Oxisols. PC2 was mostly correlated with Sr, K, and Ca (Figure 3). Most samples with positive PC2 scores were Ultisols and Inceptisols, especially those with negative PC1 scores (i.e., with lower contents of Fe, Cu, and other metallic elements). PCs 1 and 2 separated weathered soils formed from gabbro, with higher contents of Fe, from those formed from gneiss, that tend to have higher Si contents in their composition (Figure 3).

Figure 3:
Principal component analysis (PCA) of the elemental contents of 240 soil samples measured by portable X-ray fluorescence (pXRF). Soil orders provided according to the Soil Taxonomy classification. PM: parent material.

The header of Table 2 shows which variables were used in the regression structure of each parameter. The stepwise-based algorithm selected the color variables L*, a*, and b*, and PC1 for the median’s regression structure (µ), with variable b* requiring a non-linear smoothing function. For the coefficient of variation (σ), variables a*, b*, and PC1 were selected. For the regression structures of skewness and kurtosis (ν and τ), only PC2 was selected, and a smoothing function was required for kurtosis.

Table 2:
Descriptive statistics of covariates selected to fit a GAMLSS model to predict soil organic carbon (SOC) using 240 soil samples collected across the study area, located in Lavras, Minas Gerais, Brazil. SD: standard deviation; IQR: interquartile range; L*, a*, b*: color variables of the CIELAB color space; PC: principal component; µ: median; σ: coefficient of variation; ν: skewness; τ: kurtosis.

Statistic	Response variable	------------------------- Covariates -------------------------
	SOC (%)	L* (µ)	a* (µ, σ)	b* (µa, σ)	PC1 (µ, σ)	PC2 (ν, τb)
Mean	1.85	46.86	15.38	23.13	0.00	0.00
SD	1.52	9.17	4.25	4.39	2.70	2.05
Min.	0.06	21.91	4.06	11.70	-4.27	-3.37
Q25	0.70	39.84	12.90	20.35	-2.10	-1.36
Median	1.22	47.05	16.27	23.10	-0.71	-0.37
Q75	2.67	53.72	18.52	26.32	1.42	1.02
Max.	6.84	72.02	23.11	33.41	9.07	9.63
IQR	1.97	13.85	5.59	5.95	3.51	2.37
Skewness	1.18	0.05	-0.57	-0.05	1.08	1.48
Kurtosis	0.70	-0.40	-0.36	-0.31	0.95	3.56

• a A P-spline smoothing function was used with variable b* in the median’s regression structure; b A P-spline smoothing function was used with variable PC2 in the kurtosis’ regression structure.

SOC content presented a positively skewed, heavy-tailed distribution with mean and median values of 1.85 and 1.22%, respectively (Table 2). SOC contents were lower in samples with high brightness and high values of a* and b* (Figure 4), which occurred mostly on samples collected below 20 cm. Samples with dark and pale colors (lower values of L*, a*, and b*) presented higher SOC contents and were found more often at the surface (0-20 cm).

Figure 4:
CIELAB color parameters measured in 240 soil samples collected in three depths across the study area, located in Lavras, Minas Gerais, Brazil. The color of each point corresponds to the real color of the sample as measured by the Nix color sensor.

Model residuals

The fitted distribution was adequate for the modeled data (Figure 5). The normalized quantiles of residuals matched well the theoretical quantiles of the normal distribution (Figure 5a). The worm-plot did not show any large deviation from what is expected for a normally distributed variable (Figure 5b). The moment bucket plot showed that most residuals were found within the 95% confidence interval (Figure 5c), denoting that residuals are expected to have zero excess kurtosis and skewness, and therefore are likely to follow a normal distribution. There was a higher density of fitted values below 2.5% SOC, which was due to the right-tailed nature of the data (Figure 5d), but no abnormally large quantile residuals were observed.

Figure 5:
Residual analyses of the fitted GAMLSS created to predict soil organic carbon (SOC). Analyses include (a) a Q-Q plot, (b) a worm-plot (Buuren & Fredriks, 2001), (c) a moment bucket plot (Bastiani et al., 2022), and (d) a plot with normalized quantile residuals plotted against fitted values (d). The gray area of the worm-plot (b) represents 95% confidence interval for the normal distribution. In the bucket plot (c), it represents the 95% interval of the Jarque-Bera test (Jarque & Bera, 1987), wherein excess kurtosis and skewness are expected to be zero.

Model coefficients

Depth presented the biggest effect on the median of SOC contents (Table 3). At the subsurface, SOC decreases around 70% according to the model, as indicated by the exponential coefficients 0.27 or 0.31 when samples belong to depths 20-60 and 60-100 cm, respectively (Table 3). The color variable L*, which represents color brightness, had a negative effect on SOC content. SOC decreases by 2.8% per unit increase of L* (exp(µ)=0.972). Redness (a*) also had a negative effect, showing that SOC decreases by 5.3% per unit increase of a* (exp(µ)=0.947). The yellowness of samples (b*) had a negative effect on SOC contents for values lower than 18 (Figure 6), that is, samples with lower yellowness are expected to have SOC contents below the median. Conversely, a positive effect was observed for b* larger than 28. For intermediate values of b*, the effect was positive between 18 and 23, and negative between 23 and 28.

Table 3:
Coefficients of each parameter of the Box-Cox power exponential (BCPE) distribution modeled by GAMLSS. The model was fitted to predict soil organic carbon (SOC) using 240 soil samples collected in three depths across the study area, located in Lavras, Minas Gerais, Brazil. The depth group 0-20 cm was used as the reference category. L*, a*, b*: color variables of the CIELAB color space; PC: principal component; µ: median; σ: coefficient of variation; ν: skewness; ν, τ: kurtosis.

Coefficient	Link function	Intercept	Depth 20-60 cm	Depth 60-100 cm	L*	a*	b*	PC1	PC2
µ	Loga	2.818	-1.290	-1.171	-0.029	-0.055	pbb	0.044	-c
exp(µ)		16.737	0.275	0.310	0.972	0.947	pb	1.045	-
σ	Log	-1.113	0.460	0.244	-	0.061	-0.036	-0.061	-
exp(σ)		0.328	1.584	1.277	-	1.063	0.965	0.941	-
ν	Identity	0.537	-	-	-	-	-	-	-0.036
τ	Log	0.790	-	-	-	-	-	-	pb
exp(τ)		2.204	-	-	-	-	-	-	pb

• a Parameters obtained using the log link function were presented with their respective exponential to facilitate interpretations, except for τ, which had no interpretable parameter; b “pb” indicates the parameter used a P-spline smoothing function, and therefore interpretable coefficients are not available (refer to Figure 6 to interpret these covariates); c Blank cells indicate that the variable was not selected by the stepwise-based algorithm to model the respective parameter.

Figure 6:
P-splines smoothing functions fitted by a GAMLSS model created to predict soil organic carbon (SOC). The color variable b* was fitted for the regression structure of parameter µ, and variable PC2, for parameter τ (Table 3). Zero indicates no effect of the variable in the predicted median, while values higher or lower than zero indicate positive and negative effects, respectively. Rugged lines in the x-axis represent available samples.

PC1, which represented parent material (Figure 3), had a positive effect on predictions. SOC content increases by 4.5% per unit increase of the PC1 score. Higher PC1 score values are related to higher contents of the elements Fe, Cu, V, Pb, Cr, and Ti (Figure 3), which are found in higher concentrations in mafic minerals (minerals in which Fe, Mg and other metallic elements are more abundant). Conversely, negative PC1 values are related to Si, K, and Ca (Figure 3), indicating higher influence of felsic minerals (Si-rich minerals). As PC1 becomes more negative, SOC content decreases 4.31% per unit of the PC1 score (the inverse of exp(µ)). In general, coefficients indicated that SOC content is higher at the surface and in soils formed from gabbro. Conversely, as the depth, brightness, and redness of samples increase, lower SOC content is to be expected.

Similar to what was observed for the median, depth also had the biggest effect on the coefficient of variation σ (Table 3), which is a measure of the distribution’s dispersion. Dispersion increased by 58.4% for samples collected at depth 20-60 cm, and by 27.7% for those collected at 60-100 cm. The redness of soils (a*) had a positive effect on dispersion, increasing it by 6.3% for every unit increase of a*. The yellowness of soils (b*) and PC1 had a negative effect on dispersion (a 3.5 and 5.9% decrease per unit increase, respectively). Hence, the dispersion of SOC content is expected to be lower at the surface and in yellower soil samples, and higher in redder samples. Brightness (L*) was not selected by the stepwise-based algorithm for the estimation of dispersion (Table 2).

The value of the intercept obtained in the regression structure of the skewness parameter ν suggests that the distribution has positive skewness for any covariate combination (ν<1) (Rigby & Stasinopoulos, 2004). Although PC2 was selected by the stepwise-based algorithm, it was not significant at the 5% level. Nonetheless, this variable was retained in the final model. This decision is consistent with the statistical literature, including Lee et al. (2016), who suggest that removing a non-significant variable solely based on its p-value, after a variable selection procedure, may not be the most appropriate practice.

In the final fitted regression structure for kurtosis (τ), PC2 was the only significant variable, introduced through a smoothing function that is roughly negatively linear (Figure 6). Nonetheless, as PC2 increases, the number of extreme observations in the tails of the distribution tends to rise (Figure 6), since the distribution is leptokurtic for values of τ<2, as stated by Rigby and Stasinopoulos (2004). The findings for the ν and τ regression structures support the use of the BCPE distribution, which effectively captures asymmetric and leptokurtic behaviors identified in Table 2.

The effect of soil color on organic carbon models

Depth had a notable effect on SOC estimations, both on the median (µ) and on the coefficient of variation (σ). The influence of depth was expected since the main sources of SOC are the deposition of plant litter in the soil surface, the decomposition of roots, and the biological activity of the microfauna in the rhizosphere; all of which are concentrated within 20 cm from the soil’s surface (Weil & Brady, 2014). Such logical importance is also confirmed and explored by other studies using color parameters to model SOC, which found depth to be an important variable, such as the one conducted by (Mikhailova et al., 2017), who reported that adding depth in regression models that used color variables as inputs improved carbon estimation accuracy. Similarly, Faria et al., 2022 trained Random Forest models with color variables to predict SOC and found that the most important variable was depth.

The increase of dispersion (σ) with depth is likely due to the erratic occurrence of SOC in deeper soil layers. As opposed to surficial layers, where SOC is almost always more abundant, the SOC present in subsurface layers is added by the growth of plant roots, by the leaching of soluble organic compounds, or through bioturbation (activity of the microfauna). The decrease of SOC with depth depends on the effect of relief, organisms, climate, parent material, and other soil formation processes, all of which are controlled by chemical, physical, and biological factors, including plant species, root architecture and composition, soil porosity and structure, soil fertility, mineralogy, among others. Thus, SOC contents in deeper layers tend to be more uncertain and erratic.

Soils containing organic compounds are generally darker due to organic polymers and humic acids absorbing light across most of the visible wavelengths (Ben-Dor et al., 2018). It is therefore expected that the increase of color brightness (L*) will have a negative effect on SOC content, which has been confirmed by other works that used color variables to estimate SOC (Stiglitz et al., 2017). Redness (a*) is also reduced by the presence of SOC. The red color of soils stems from hematite, an iron oxide with many characteristic absorption peaks in the visible range (Viscarra Rossel et al., 2010 ). SOC, however, reduces the albedo (overall reflectance) across all the visible range; that is, it makes colors darker and paler, therefore lessening the reddening effect of iron oxides (Mancini et al., 2023). At the same time, Fe-rich soils tend to be more weathered and have more clay content due to the residual accumulation of Fe- and Al-rich minerals, although dependent also on the mineralogy of the parent material. Finer particles have higher surface area, requiring more SOC to effectively darken the soil. This may explain why dispersion increases with the increase of redness (6.3% per unit), as the coloring effect of SOC is more erratic in clay-rich, red soil samples due to their Fe-rich mineralogy and finer particle size. Clay-sized hematite (Fe₂O₃) is a strong pigmenting agent and attenuates the darkening effect of SOC, while clay-rich kaolinitic soils are easily darkened by SOC. In future research, the method utilized in this study could be used to model the interaction between the color of minerals and SOC contents.

The positive effect of PC1 (4.5% SOC per unit increase) is related to the higher contents of SOC in soils developed from gabbro. From the Oxisols studied here, those derived from gabbro are richer in Fe, Pb, Ti, and other metallic elements, and presented a positive PC1 score. There is a vast literature about the processes that favor the protection of organic matter in oxide-rich soils (Li et al., 2023). Two relevant processes known to occur are the bonding between the surface of Fe oxides and hydroxyl and carboxyl groups in organic matter (Wiseman & Püttmann, 2006), forming organo-mineral assemblages less accessible to microorganisms, and the physical isolation caused by the micro-aggregation of these assemblages, occluding organic compounds in compartments physically inaccessible to decomposing organisms (Kleber et al., 2021). The increase of 4.5% SOC content per unit increase of the PC1 score seems to indicate this positive effect between SOC content and the protection conferred by oxidic soils developed from mafic parent materials (gabbro, in this case). Similarly, it seems reasonable that predictions would have lower dispersion in samples from soils where SOC is more stabilized, such as those derived from gabbro, richer in oxidic clays, which may explain why the increase of PC1 (indicative of gabbro-derived soils) is related to a decrease in dispersion (5.9% per unit increase).

Conclusions

SOC was roughly 70% lower in subsurface horizons and exponentially decreased in bright red horizons. Samples with positive PC1 scores were related to gabbro and associated with higher SOC contents. Positive scores of PC2 were correlated with samples derived from felsic minerals. Soil sensing and GAMLSS helped explain the effects of color, depth, and parent material in the variation of SOC. This enables more comprehensive modeling, increasing sample coverage and improving understanding of the physical processes governing SOC distributions.

Acknowledgments

The authors would like to thank the funding agencies CNPq (305215/2024-8; 307532/2022-4; 404983/2023-5; 140697/2024-1), CAPES (PROEX N° 88887.924197/2023-00), and FAPEMIG for financially supporting this research. L.R. Nakamura was partially supported by FAPEMIG through the grant RED-00133-21.

Data Availability Statement

Data available upon request to authors.

References

• Alvares, C. A. et al. (2013). Köppen’s climate classification map for Brazil. Meteorologische Zeitschrift, 22(6):711-728.
• Bastiani, F. D. et al. (2022). Bucket plot: A visual tool for skewness and kurtosis comparisons. Brazilian Journal of Probability and Statistics, 36(3):421-440.
• Ben-Dor, E. et al. (2018). Characterization of soil properties using reflectance spectroscopy. In P. S. Thenkabail., J. G. Lyon., & A. Huete. (eds.). Fundamentals, sensor systems, spectral libraries, and data mining for vegetation (2nd ed). CRC Press, (pp.187-247).
• Buuren, S. v., & Fredriks, M. (2001). Worm plot: A simple diagnostic device for modelling growth reference curves. Statistics in Medicine, 20(8):1259-1277.
• Dunn, P. K., & Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics, 5(3):236-244.
• Eilers, P. H. C. et al. (2015). Twenty years of p-splines. Sort-statistics and Operations Research Transactions, 39(2):149-186.
• Faria, A. J. G. de. et al. (2022). Prediction of soil organic matter content by combining data from Nix Pro^TM color sensor and portable X-ray fluorescence spectrometry in tropical soils. Geoderma Regional, 28:e00461.
• Jarque, C. M., & Bera, A. K. (1987). A test for normality of observations and regression residuals. International Statistical Review, 55(2):163-172.
• Jolliffe, I. T., & Cadima, J. (2016). Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2065):20150202.
• Kleber, M. et al. (2021). Dynamic interactions at the mineral-organic matter interface. Nature Reviews Earth & Environment, 2(6):402-421.
• Kneib, T. et al. (2023). Rage against the mean - A review of distributional regression approaches. Econometrics and Statistics, 26:99-123.
• Lamichhane, S. et al. (2019). Digital soil mapping algorithms and covariates for soil organic carbon mapping and their implications: A review. Geoderma, 352:395-413.
• Lee, J. D. et al. (2016). Exact post-selection inference, with application to the lasso. The Annals of Statistics, 44(3):907-927.
• Li, Q. et al. (2023). Interactions between organic matter and Fe oxides at soil micro-interfaces: Quantification, associations, and influencing factors. Science of the Total Environment, 855:158710.
• Mancini, M. et al. (2023). Digital morphometrics and genesis of soils with buried horizons and lithological discontinuities in southeastern Brazil. Geoderma Regional, 32:e00612.
• Mancini, M. et al. (2019). Parent material distribution mapping from tropical soils data via machine learning and portable X-ray fluorescence (pXRF) spectrometry in Brazil. Geoderma, 354:113885.
• McBratney, A. et al. (2019). Pedometrics timeline. Geoderma, 338:568-575.
• Mikhailova, E. A. et al. (2017). Predicting soil organic carbon and total nitrogen in the Russian Chernozem from depth and wireless color sensor measurements. Eurasian Soil Science, 50:1414-1419.
• Pouladi, N. et al. (2023). Digital mapping of soil organic carbon using remote sensing data: A systematic review. Catena, 232:107409.
• Ramires, T. G. et al. (2021). Validation of stepwise-based procedure in GAMLSS. Journal of Data Science, 19(1):96-110.
• Ribeiro, D. et al. (2024). Effectiveness of proximal and remote sensing to trace soils spatial variability: A study case in a Brazilian Cerrado coffee farm. Ciência e Agrotecnologia, 48:e013824.
• Rigby, R. A., & Stasinopoulos, D. M. (2004). Smooth centile curves for skew and kurtotic data modelled using the Box-Cox power exponential distribution. Statistics in Medicine, 23(19):3053-3076.
• Rossel, R. A. V., & Behrens, T. (2010). Using data mining to model and interpret soil diffuse reflectance spectra. Geoderma, 158(1):46-54.
• Santos, H. G. et al. (2025). Sistema brasileiro de classificação de solos (6th ed). Brasília, DF, Brasil: Embrapa.
• Schmidinger, J., & Heuvelink, G. B. M. (2023). Validation of uncertainty predictions in digital soil mapping. Geoderma, 437:116585.
• Shi, T. et al. (2023). Digital mapping of heavy metals in urban soils: A review and research challenges. Catena, 228:107183.
• Silva, S. H. G. et al. (2021). pXRF in tropical soils: Methodology, applications, achievements and challenges. Advances in Agronomy, 167:1-62.
• Soil Survey Staff. (2014). Keys to soil taxonomy (12th ed). USDA-NRCS. Available in: https://www.nrcs.usda.gov/sites/default/files/2022-09/Keys-to-Soil-Taxonomy.pdf
  » https://www.nrcs.usda.gov/sites/default/files/2022-09/Keys-to-Soil-Taxonomy.pdf
• Singh, R. et al. (2025). Enhancing soil organic carbon estimation with generative AI and Nix color sensor. Scientific Reports, 15:40628.
• Stasinopoulos, M. D. et al. (2024). Generalized additive models for location, scale and shape: a distributional regression approach, with applications Cambridge University Press.
• Stiglitz, R. et al. (2017). Using an inexpensive color sensor for rapid assessment of soil organic carbon. Geoderma, 286:98-103.
• Viscarra Rossel, R. A. et al. (2010). Mapping iron oxides and the color of Australian soil using visible-near-infrared reflectance spectra. Journal of Geophysical Research: Earth Surface, 115:04031.
• Walkley, A., & Black, I. A. (1934). An examination of the Degtjareff method for determining soil organic matter, and a proposed modification of the chromic acid titration method. Soil Science, 37:29-39.
• Weil, R., & Brady, N. (2014). The nature and properties of soils (15th ed). New York: Pearson, 1104 pp.
• Wiseman, C. L. S., & Püttmann, W. (2006). Interactions between mineral phases in the preservation of soil organic matter. Geoderma, 134(1-2):109-118.
• Yang, L., & Li, X. (2026). Digital soil mapping in the era of big data and artificial intelligence. Annals of GIS, 441-465.

Editor de seção:

Renato Paiva http://orcid.org/0000-0001-5107-0285

^*

Corresponding author: [e-mail]

Publication Dates

• Publication in this collection
  19 Oct 2026
• Date of issue
  2026

History

• Received
  16 Mar 2026
• Accepted
  13 Aug 2026