Integrating Cacao Physicochemical-Sensory Profiles via Gaussian Processes Crowd Learning and Localized Annotator Trustworthiness

1. Introduction

Cocoa-based products carry profound cultural significance, offer essential nutritional value, and underpin a multibillion-dollar global industry that supports the livelihoods of millions in rural communities [1]. Given that consumer acceptance is largely driven by the perception of sensory attributes such as aroma, flavor, and texture, comprehensive sensory profiling has become critical, not only for quality assurance and regulatory compliance, but also for the preservation of geographical indication and origin certification schemes [2,3,4]. However, the human evaluations that form the basis of these sensory assessments are intrinsically subjective, with panelist perceptions varying according to factors such as experience, fatigue, and environmental conditions during tasting [5,6].

To address this inherent variability, researchers increasingly integrate diverse instrumental measurements with advanced chemometric and machine learning approaches. Techniques such as ultra-high-performance liquid chromatography and high-resolution mass spectrometry–based sensomics enable the identification of molecular determinants linked to both sensory quality and geographic origin [7]. In parallel, multivariate statistical tools like Orthogonal Partial Least Squares Discriminant Analysis (OPLS-DA) have proven effective in associating physicochemical variables—such as pH, °Brix, and polyphenol content—with specific flavor descriptors [8]. High-throughput metabolomics workflows, including platforms like FlavorMiner, facilitate the extraction of latent flavor signatures from complex, large-scale datasets [9]. Moreover, mixed-effects modeling frameworks have been employed to explicitly account for panelist-level noise and inter-individual variability during critical post-harvest stages such as fermentation and drying [10].

Nevertheless, sensory datasets often remain limited in size and characterized by high levels of noise, primarily due to the substantial time and cost associated with assembling and maintaining trained sensory panels. Moreover, inter-annotator reliability tends to fluctuate depending on the intrinsic chemical complexity of the samples under evaluation [11,12,13,14,15]. Within this evolving methodological landscape, basic physicochemical parameters, such as pH, fat content, polyphenol concentration, and key volatile compounds, have consistently demonstrated strong and reproducible correlations with sensory attributes [6]. As a result, these indicators offer practical and reliable proxies for predicting cocoa quality, particularly in contexts where comprehensive sensory evaluation is unfeasible.

Modern machine learning increasingly addresses the challenges posed by low-quality datasets—characterized by limited size, noise, sparsity, and subjectivity—by employing robust, noise-tolerant models and integrative frameworks. One notable trend is the application of deep learning methods with built-in noise-resilience or regularization mechanisms designed to extract reliable patterns from small and unreliable datasets. These techniques often include ensemble learning, dropout, and variational inference methods. For instance, in the domain of metabolomics-driven prediction in vegetable foods, deep learning models have shown strong performance even when faced with noisy and missing data by effectively integrating heterogeneous sources [16]. In food science applications, the data quality issues are exacerbated by the use of sensor-based and crowd-sourced annotations, often yielding highly subjective and inconsistent datasets. Advanced methods leveraging spectroscopy and hyperspectral imaging integrated with Artificial Intelligence (AI) have been used to address such limitations. These methods include deep architectures that can model high-dimensional, noisy inputs effectively, as shown in grain adulteration detection and cereal quality assessment tasks [17]. Furthermore, food authentication technologies such as electronic nose and tongue systems have increasingly adopted machine learning strategies that not only enhance detection accuracy but also account for sensor noise and low signal-to-noise ratios in complex food matrices [18]. Also, systematic reviews underscore that robustness to overlapping classes, outliers, and subjective assessments is critical in these contexts, often requiring hybrid models or pre-processing stages for noise filtration [19]. These advancements are particularly relevant for sensory profiling applications, where subjective human input is prevalent.

In turn, recent advances in machine learning increasingly rely on crowdsourced or multi-annotator datasets, especially in domains where expert labeling is prohibitively expensive or time-consuming [20,21]. Traditional approaches often adopt label aggregation methods—such as Majority Voting (MV) or Expectation-Maximization (EM)—to create a single consensus label from multiple annotations [22]. Still, these techniques assume homogeneous annotator reliability and independence, which rarely holds in real-world scenarios. To overcome these limitations, modern frameworks now model annotator behavior as a function of the input space, using probabilistic tools like Gaussian Processes (GPs) [23], or deep neural networks that learn both the task and annotator reliability simultaneously [24]. One prominent example is the Correlated Chained Gaussian Processes with Generalized Cross Entropy (CCGP-GCE), which jointly models label noise and inter-annotator dependencies while leveraging a robust loss function to mitigate the impact of outliers and adversarial annotators [20].

Likewise, a growing trend in crowd learning research is the estimation of annotator-specific trustworthiness that varies across the input domain. Unlike earlier models that assigned fixed reliability scores, recent methods estimate annotator accuracy dynamically, using latent variables and input-conditioned reliability functions [25]. For instance, CCGP-GCE incorporates sigmoid-transformed latent functions to model each annotator’s reliability in a Bayesian framework, enabling soft decisions between trustworthy and noisy labels at each instance [20]. This is particularly effective when combined with the Generalized Cross Entropy loss, which blends the robustness of mean absolute error (MAE) with the fast convergence of cross-entropy (CE) [26]. Nevertheless, key challenges persist: scalability is limited by the computational complexity of Gaussian Process inference, and the lack of ground truth hampers the evaluation of trustworthiness estimations. Moreover, deep learning models that disregard inter-annotator correlations often underperform compared to Bayesian alternatives that explicitly encode annotator dependencies. As such, state-of-the-art methods increasingly favor hybrid approaches that integrate probabilistic modeling, structured priors, and noise-robust objectives to learn from complex, inconsistent crowdsourced data.

In particular, regression tasks involving multiple annotators, labeler variability stems from a range of factors including differing levels of expertise, perceptual biases, and context-dependent interpretations. Rather than treating this variability as random noise, recent studies suggest that it often reflects systematic patterns tied to task structure and data representation [27]. For example, traditional agreement metrics fall short in continuous signal domains, where finer-grained discrepancies among annotators are meaningful and can be predictive [27]. Automated methods have also been introduced to inspect and interpret such variability, revealing that annotator behavior itself may encode relevant information [28]. Furthermore, annotator consistency—commonly used as a proxy for reliability—is not constant across the input space. Instead, it varies with the characteristics of the data, underscoring its inherently non-stationary nature [29].

To address the pressing challenge of integrating heterogeneous and noisy sensory data with physicochemical profiles in food science, we propose a novel Multi-Annotator Regression framework based on Correlated Chained Gaussian Processes (MAR-CCGP). Unlike prior approaches, MAR-CCGP explicitly models input-dependent annotator trustworthiness while simultaneously estimating continuous perceptual scores from physicochemical inputs. This dual modeling capability is especially critical in sensory evaluation contexts, where annotations are sparse, subjective, and non-stationary. Our proposal unfolds across the following key stages:

-

Construction of the LUKER-CACAO database, aligning standardized physicochemical measurements with sensory annotations from multiple expert panelists.

-

Application of the MAR-CCGP framework to learn latent ground truth sensory scores and context-dependent annotator reliability through a shared latent factor model from noisy multiple annotator regression data.

-

Development of a localized annotator trust score, leveraging the model’s posterior distribution to assess reliability per annotator and per sample.

We conduct controlled experiments on both real-world and semi-synthetic datasets to benchmark performance and interpretability. Our experiments span the proprietary LUKER-CACAO dataset (five physicochemical inputs, and eight sensory descriptors rated by five annotators) and multiple semi-synthetic benchmarks derived from UCI regression datasets, where structured, region-specific noise profiles were simulated to emulate non-stationary annotator behavior. We benchmark MAR-CCGP against GPR-GT (GP on ground truth) [30], GPR-AVG (GP on consensus averages) [23], and LKAAR (a localized kernel alignment model for annotator relevance) [31]. Results demonstrate that MAR-CCGP consistently achieves superior predictive performance—often approaching the oracle—even in the absence of ground-truth labels. It recovers localized annotator trust profiles and outperforms consensus-only and partially localized models by capturing both inter-annotator correlations and input-conditioned noise structures. Overall, MAR-CCGP offers a suitable and interpretable solution for learning from subjective, sparse, and inconsistent annotations, with significant implications for robust modeling in food quality control.

Section 2 details the materials and methods employed in the study. Section 3 introduces the experimental setup. Section 4 presents and discusses the results. Finally, Section 5 provides the concluding remarks and future perspectives.

2. Materials and Methods

2.1. Casa Luker - Cacao Physicochemical-Sensory Dataset (LUKER-CACAO)

The dataset used in this study was developed in collaboration with Casa Luker (Luker Chocolate https://lukerchocolate.com/en/) and contains a collection of sensory and physicochemical measurements derived from cacao-based product evaluations. These data originate from routine quality control and research procedures applied to a wide variety of product types such as cocoa liquors, dark chocolate, milk chocolate, white chocolate, and cocoa powders.

The original dataset was compiled from two distinct sources: 15 quarterly technical reports containing physicochemical measurements, and 704 sensory evaluation session files. The physicochemical data were obtained through standardized instrumental analyses, while the sensory sessions involved trained panelists who provided evaluations of cacao production samples according to established sensory protocols. To construct a unified and analyzable dataset, all reports were systematically parsed using automated data extraction pipelines. The relevant tabular data were consolidated into two structured matrices: one representing the physicochemical descriptors (input features) and the other capturing the sensory evaluations (multi-annotator outputs). The sensory matrix adopts a structured format in which each row corresponds to a product sample and each column represents a specific sensory descriptor provided by an individual annotator, reflecting the inherent variability and sparsity of crowd-sourced annotations.

A subset of variables from both data modalities was subsequently selected to ensure consistency and relevance. For the physicochemical domain, five variables were retained based on their completeness across samples and their potential to capture key aspects of product structure and composition:

-

Moisture: It refers to the amount of water present in a cocoa-based product, expressed as a percentage of the product’s total weight (%). Moisture content significantly impacts physical attributes such as texture, shelf-life, and microbial stability. More critically, it modulates the release and perception of flavor compounds during consumption. Variations in moisture influence the volatilization of aroma molecules and alter the way flavors are experienced in the mouth, thereby reshaping the sensory profile in terms of intensity, balance, and mouthfeel [32]. In the dataset, moisture values range from 0 to approximately 200%.

-

Fat content: It represents the proportion of lipids present in cocoa-based products, expressed as a percentage of total weight (%). Its modulation influences physical properties like viscosity, structure, and film formation capacity. These physical changes affect how the product interacts with the mouth during consumption, thereby altering the sensory profile by modifying lubrication, mouthfeel, and perception of flavor release [33]. In the dataset, values range from 0 to 60%.

-

Granulometry: It measures the size and distribution of solid particles in cocoa-based products, expressed in micrometers ($\mu$m). Granulometry influences how particles interact and pack together, altering viscosity, flow behavior, and ultimately the perception of mouthfeel during consumption. Changes in granulometry reshape the sensory profile by modifying sensations like smoothness, thickness, and creaminess, which are critical for consumer acceptance [34]. Granulometry values range from 0 to 58 $\mu$m.

-

Plastic viscosity: Measures the resistance of the product to flow after yielding has occurred, expressed in Pascal-seconds (Pa·s). It also reflects how easily the material continues to deform under applied shear during oral processing. Variations in plastic viscosity affect the sensory profile by altering the perceived thickness, smoothness, and creaminess during consumption. These changes shape the overall mouthfeel, influencing whether the product is experienced as rich, velvety, or fluid [34]. Observed values range from 0 to approximately 10.5 Pa·s.

-

Yield stress: Represents the minimum force required to initiate flow in the product, expressed in Pascals (Pa). It is closely related to the structural integrity of the product before deformation starts. Variations in yield stress affect the sensory profile by altering initial mouthfeel sensations such as firmness and body, shaping the consumer’s perception of texture at the start of consumption [34]. Yield stress values range from 0 to approximately 62 Pa.

Likewise, from the sensory domain, eight attributes were selected based on their relevance to flavor perception and the density of available annotations. These attributes represent core dimensions of sensory quality. In trained descriptive analysis of cocoa-based products, sensory attributes are evaluated by expert panels using intensity scales ranging from 1 (absence) to 10 (maximum perceived intensity). These attributes capture key sensory dimensions affecting product quality, processing control, and consumer acceptance [35]:

-

Acidity: It refers to the perception of sourness resulting from organic acids formed during fermentation. When present at appropriate levels, acidity can enhance brightness and complexity; however, excessive acidity is considered a defect, particularly in fine chocolate [36].

-

Bitterness: It reflects the presence of alkaloids (primarily theobromine) and polyphenols, compounds inherent to cocoa. While some degree of bitterness is characteristic and desirable, excessive levels can disrupt sensory balance and negatively impact consumer acceptance [37].

-

Aroma: It encompasses volatile compounds responsible for cocoa’s characteristic smells (fruity, floral, roasted), highly influenced by fermentation and roasting [32].

-

Astringency: It refers to the drying, puckering sensation caused by interactions between polyphenols and salivary proteins. While moderate astringency can contribute positively to mouthfeel and complexity, excessive levels are perceived as unpleasant and may negatively impact sensory acceptance. [37].

-

Sweetness: It reflects sugar content, critical to balancing bitterness and acidity for overall flavor harmony.

-

Hardness: It describes the resistance during biting or deformation, influenced by fat content, tempering, and particle size.

-

Melting speed: It reflects how quickly the product liquefies in the mouth, depending on fat composition and tempering. Faster melting generally enhances flavor release and mouthfeel [35].

-

Global impresssion: It summarizes overall product quality, integrating flavor, aroma, and texture into a single judgment.

Of note, the analytical procedures used to obtain these physicochemical and sensory measurements are standardized and traceable to established methods. Table 1 summarizes the official protocols used for each selected measurement type.

Next, to ensure adequate representation and consistency within the multi-annotator framework, the analysis was further restricted to a subset of five annotators, identified by the codes 135, 154, 155, 160, and 179. These annotators were selected based on their high annotation coverage across the chosen sensory attributes, enabling more reliable modeling of expert-specific patterns and behaviors. Namely, after merging we built an input–output, multi-annotator dataset with 1886 samples, five physicochemical features, and eight sensory outputs from five experts. Figure 1 presents the annotation coverage heatmap, grouping columns by sensory variable and ordering annotators within each block; black squares indicate present labels and white squares missing ones. This visualization highlights that certain variables include complete samples lacking annotations from any expert. Additionally, it reveals systematic patterns in annotator behavior—for instance, the fourth column in each block consistently shows fewer labels provided for attributes such as bitterness and aroma. Figure 3 further quantifies these patterns, providing a detailed overview of inconsistencies in annotator behavior: panel (a) shows the distribution of how many annotators labeled each sample—most receive labels from two to three experts—while panel (b) depicts the distribution of the standard deviation in annotator scores per sample–variable pair, revealing a concentration of low disagreement but also a long tail of cases with high annotator variability, emphasizing the challenge of modeling subjective sensory assessments.

For each selected sensory attribute, we constructed a dedicated dataset by intersecting the cleaned sensory and physicochemical records. The data were restricted to the previously selected subset of annotators and physicochemical features. To ensure that the learning models are trained on complete and reliable input–output pairs, we only retained samples that satisfied two conditions: (i) no more than one missing value among the physicochemical features, and (ii) no more than one missing annotation across the selected annotators for the target attribute.

This filtering strategy strikes a balance between strict completeness and reasonable sample retention, allowing us to preserve useful samples without compromising model integrity. Table 2 summarizes the percentage of available labels for each annotator across all sensory attributes. Table 3 shows the completeness of each physicochemical input feature, along with the final number of retained samples for each sensory regression task.

Figure 2.

Figure 2.

Figure 3. Inconsistencies in annotator behavior. (a) Shows, for each sample, how many annotators provided at least one label across the eight sensory attributes, revealing a typical coverage of two to three annotators per sample.(b) Presents a kernel density estimate of the standard deviation in annotator scores for each sample–variable pair, quantifying how much annotators typically disagree. Most annotations exhibit standard deviations below 2.0, but a notable long tail indicates the presence of substantial inter-annotator variability in a subset of evaluations.

Figure 3. Inconsistencies in annotator behavior. (a) Shows, for each sample, how many annotators provided at least one label across the eight sensory attributes, revealing a typical coverage of two to three annotators per sample.(b) Presents a kernel density estimate of the standard deviation in annotator scores for each sample–variable pair, quantifying how much annotators typically disagree. Most annotations exhibit standard deviations below 2.0, but a notable long tail indicates the presence of substantial inter-annotator variability in a subset of evaluations.

After applying the completeness constraints and constructing the final datasets for each sensory attribute, a final imputation step was performed to handle any remaining missing values in the physicochemical input variables. We employed the Iterative Imputer from scikit-learn, which estimates missing values by sequentially modeling each variable as a function of the others. This approach preserves inter-variable correlations and provides a flexible, data-driven completion method suitable for small datasets [43]. No imputation was applied to the sensory labels. Missing annotations were treated as unobserved and marginalized during training.

The overall data curation workflow is summarized in Figure 3. This diagram outlines the end-to-end process from raw data acquisition through parsing, variable selection, and completeness filtering, culminating in the construction of the LUKER-CACAO dataset. It highlights the progressive reduction and refinement of data needed to obtain high-quality, multi-annotator datasets suitable for learning reliable mappings between physicochemical properties and sensory perceptions.

Figure 4. Pipeline for the construction of the LUKER-CACAO database.

Figure 4. Pipeline for the construction of the LUKER-CACAO database.

2.2. Correlated Chained Gaussian Processes (CCGP)

Let a Gaussian Process (GP) be a collection of random variables $f\left(\mathbf{x}\right)\in \mathbb{R}$ indexed by the input samples $\mathbf{x}\in {\mathbb{R}}^P$ such that any finite number of them follows a joint Gaussian distribution. Regarding this, a GP is completely specified by its mean ${\mu }_x=\mathbb{E}\left\{f\left(\mathbf{x}\right)\right\}$ (${\mu }_x=0$ for simplicity) and covariance $\kappa (\mathbf{x},{\mathbf{x}}^{\prime })=\mathbb{E}\left\{\left(f(\mathbf{x}-{\mu }_x)(f\left({\mathbf{x}}^{\prime }\right)-{\mu }_{x^{\prime }})\right)\right\}$, where ${\kappa }_f:{\mathbb{R}}^P\times {\mathbb{R}}^P\to \mathbb{R}$ is a kernel function, ${\mathbf{x}}^{\prime }\in {\mathbb{R}}^P$, and $\mathbb{E}\left\{·\right\}$ stands for the expectation operator [30]. Then: $f\left(\mathbf{x}\right)\sim \mathcal{GP}(0,\kappa (\mathbf{x},{\mathbf{x}}^{\prime }))$.

Now, let $\mathcal{D}=\{\mathbf{X}\in {\mathbb{R}}^{N\times P},\mathit{y}\in {\mathbb{R}}^N\}$ denote an input–output dataset consisting of N samples and P input features, where $\mathbf{X}$ represents the input matrix and $\mathit{y}$ is the corresponding vector of continuous regression targets. We define the latent function evaluations over the dataset as $\mathbf{f}={[f\left({\mathbf{x}}_1\right),f\left({\mathbf{x}}_2\right),\dots ,f\left({\mathbf{x}}_N\right)]}^{\top }\in {\mathbb{R}}^N$, where $\mathbf{f}\sim \mathcal{GP}(\mathbf{0},{\mathbf{K}}_{ff})$, being ${\mathbf{K}}_{ff}\in {\mathbb{R}}^{N\times N}$ a kernel matrix computed by evaluating the covariance function ${\kappa }_f({\mathbf{x}}_n,{\mathbf{x}}_{n^{\prime }})$ over all pairs of input points ($n,n^{\prime }=\{1,2,\dots ,N\}$).

Remarkably, GPs inherit the properties of multivariate Gaussian distributions, allowing for linear combinations of latent functions and facilitating the construction of additive models. When coupled with an appropriate link function, this framework can be naturally extended to accommodate non-Gaussian observation models, as in generalized linear model settings [44]. Moreover, in the general case of GP-based models for supervised learning, the observed data $\mathcal{D}$ is modeled by constructing a joint distribution that combines a conditional likelihood with one or more latent functions governed by independent GP priors. If each parameter of the likelihood (e.g., noise variance, mean shift) is itself modeled as a latent function over the input space, the resulting structure is known as a Chained Gaussian Process (CGP) [45]. Formally, the joint distribution over $\mathit{y}$ and J latent functions (LFs) $\widehat{\mathbf{f}}={\left[{\mathbf{f}}_1,\dots ,{\mathbf{f}}_J\right]}^{\top }\in {\mathbb{R}}^{NJ}$, conditioned on the input matrix $\mathbf{X}$, yields:

$$p(\mathit{y},\widehat{\mathbf{f}}\mid \mathbf{X})=\prod _{n=1}^{N}p(y_n\mid {\theta }_1\left({\mathbf{x}}_n\right),\dots ,{\theta }_J\left({\mathbf{x}}_n\right))\prod _{j=1}^J\mathcal{N}({\mathbf{f}}_j\mid \mathbf{0},{\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{f}}_j}),$$

(1)

where the likelihood parameters ${\theta }_j\left({\mathbf{x}}_n\right)$, with $j\in \{1,2,\dots ,J\}$, are nonlinear mappings from GP priors, e.g., ${\theta }_j\left(\mathbf{x}\right)=h_j\left(f_j\left(\mathbf{x}\right)\right)$, with $h_j:\mathbb{R}\to {\mathcal{M}}_j$ mapping each $f_j\left(\mathbf{x}\right)$ to the appropriate domain ${\mathcal{M}}_j$. In addition, ${\mathbf{f}}_j={\left[f_j\left({\mathbf{x}}_1\right),\dots ,f_j\left({\mathbf{x}}_N\right)\right]}^{\top }\in {\mathbb{R}}^N$ is an LF vector that follows a GP prior and ${\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{f}}_j}\in {\mathbb{R}}^{N\times N}$ is a kernel-based covariance matrix.

GPs, being nonparametric models, inherently scale poorly with dataset size due to their computational complexity of $\mathcal{O}\left(N^3\right)$ for a regression task defined over $\mathcal{D}$. To address this limitation, we adopt a widely used approximation strategy that introduces a set of $M\ll N$ inducing variables, denoted by ${\mathbf{u}}_j={[f_j\left({\mathbf{z}}_1^j\right),\dots ,f_j\left({\mathbf{z}}_M^j\right)]}^{\top }\in {\mathbb{R}}^M$, where ${\mathbf{Z}}_j=[{\mathbf{z}}_1^j,\dots ,{\mathbf{z}}_M^j]\in {\mathbb{R}}^P$. This sparse approximation reduces the computational cost of inference and learning to $\mathcal{O}\left(N{M}^2\right)$. Under this formulation, the CGP prior yields [46]:

$$p({\mathbf{f}}_j,{\mathbf{u}}_j)\sim \mathcal{N}\left(\left[\begin{array}{c}{\mathbf{f}}_j\\ {\mathbf{u}}_j\end{array}\right]|\mathbf{0},\left[\begin{array}{cc}{\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{f}}_j}& {\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{u}}_j}\\ {\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{f}}_j}& {\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{u}}_j}\end{array}\right]\right).$$

(2)

Here, ${\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{u}}_j}\in {\mathbb{R}}^{N\times M}$ denotes the cross-covariance computed by evaluating $k_j(·,·)$ between input samples and inducing points. Similarly, ${\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{u}}_j}\in {\mathbb{R}}^{M\times M}$ represents the covariance matrix between inducing points. By conditioning on the inducing variables ${\mathbf{u}}_j$, the marginal distribution over the latent function values ${\mathbf{f}}_j$ is given by:

$$p({\mathbf{f}}_j\mid {\mathbf{u}}_j)=\mathcal{N}\left({\mathbf{f}}_j|{\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{u}}_j}{\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{u}}_j}^{-1}{\mathbf{u}}_j,{\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{f}}_j}-{\mathbf{K}}_{{\mathbf{f}}_j{\mathbf{u}}_j}{\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{u}}_j}^{-1}{\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{f}}_j}\right),$$

(3)

and the prior on ${\mathbf{u}}_j$ is:

$$p\left({\mathbf{u}}_j\right)=\mathcal{N}({\mathbf{u}}_j\mid \mathbf{0},{\mathbf{K}}_{{\mathbf{u}}_j{\mathbf{u}}_j}).$$

(4)

In general, the posterior distribution $p(\widehat{\mathbf{f}},\widehat{\mathbf{u}}\mid \mathit{y})$ is intractable due to the non-conjugacy of the likelihood with respect to the priors ($\widehat{\mathbf{u}}={[{\mathbf{u}}_1^{\top },\dots ,{\mathbf{u}}_J^{\top }]}^{\top }\in {\mathbb{R}}^{MJ}$). Then, we approximate the posterior using variational inference with a parameterized distribution $p(\widehat{\mathbf{f}},\widehat{\mathbf{u}}\mid \mathit{y})\approx \tilde{q}(\widehat{\mathbf{f}},\widehat{\mathbf{u}})$ [23]:

$$\tilde{q}(\widehat{\mathbf{f}},\widehat{\mathbf{u}})=p\left(\widehat{\mathbf{f}}\right|\widehat{\mathbf{u}})\tilde{q}\left(\widehat{\mathbf{u}}\right)=\prod _{j=1}^{J}p\left({\mathbf{f}}_j\right|{\mathbf{u}}_j)\tilde{q}\left({\mathbf{u}}_j\right),$$

(5)

where $\tilde{q}\left({\mathbf{u}}_j\right)$ is a variational Gaussian distribution over the inducing variables and:

$$\tilde{q}\left(\widehat{\mathbf{u}}\right)=\prod _{j=1}^J\tilde{q}\left({\mathbf{u}}_j\right)=\prod _{j=1}^J\mathcal{N}\left({\mathbf{u}}_j\right|{\widehat{\mu }}_j,{\widehat{\mathbf{V}}}_j),$$

(6)

with ${\widehat{\mu }}_j\in {\mathbb{R}}^M$ and ${\widehat{\mathbf{V}}}_j\in {\mathbb{R}}^{M\times M}$. These parameters are optimized by maximizing an Evidence Lower Bound (ELBO), which provides a tractable surrogate to the marginal likelihood. Assuming that the data points ${\mathbf{x}}_n$ are independently and identically distributed, such an ELBO can be formulated as [47]:

$$\mathcal{L}=\sum _{n=1}^N{\mathbb{E}}_{\tilde{q}\left({\mathbf{f}}_1\right),\dots ,\tilde{q}\left({\mathbf{f}}_J\right)}\left[logp(y_n\mid {\theta }_1\left({\mathbf{x}}_n\right),\dots ,{\theta }_J\left({\mathbf{x}}_n\right))\right]-\sum _{j=1}^J{\mathbb{D}}_{\mathrm{KL}}\left(\tilde{q}\left({\mathbf{u}}_j\right)\parallel p\left({\mathbf{u}}_j\right)\right),$$

(7)

where ${\mathbb{D}}_{\mathrm{KL}}(·\parallel ·)$ is the Kullback-Leibler divergence and $\tilde{q}\left({\mathbf{f}}_j\right)=\int p({\mathbf{f}}_j\mid {\mathbf{u}}_j)\tilde{q}\left({\mathbf{u}}_j\right)d{\mathbf{u}}_j$.

Still, the standard CGP assumes independent GP priors for each likelihood parameter. The latter assumption is often unrealistic in multi-annotator scenarios, where annotators may share common biases or information sources [31]. To capture these dependencies, we employ a shared latent-factor structure [48]:

$${\breve{f}}_j\left({\mathbf{x}}_n\right)=\sum _{q=1}^Q{w}_{j,q}{\vartheta }_q\left({\mathbf{x}}_n\right),$$

(8)

where ${\breve{f}}_j:{\mathbb{R}}^P\to \mathbb{R}$, ${\vartheta }_q(·)\sim \mathcal{GP}(0,{\kappa }_q(·,·))$, and $w_{j,q}\in \mathbb{R}$. Let $\breve{\mathbf{\theta }}={[{\mathbf{\theta }}_1^{\top },\dots ,{\mathbf{\theta }}_J^{\top }]}^{\top }\in {\mathbb{R}}^{NJ}$ be the sample-dependent model parameter vector, $\breve{\mathbf{f}}={[{\breve{\mathbf{f}}}_1,\dots ,{\breve{\mathbf{f}}}_J]}^{\top }\in {\mathbb{R}}^{NJ}$, ${\breve{\mathbf{f}}}_j={[{\breve{f}}_j\left({\mathbf{x}}_1\right),\dots ,{\breve{f}}_j\left({\mathbf{x}}_N\right)]}^{\top }\in {\mathbb{R}}^N$, and $q\in \{1,2,\dots ,Q\}.$ herefore, the Correlated Chained Gaussian Process (CCGP) framework naturally emerges when the joint distribution in equation 1 is reformulated to explicitly model dependencies across outputs through a chained, conditionally correlated structure:

$$p(\mathit{y},\breve{\mathbf{f}},\breve{\mathbf{u}}\mid \mathbf{X})=p(\mathit{y}\mid \breve{\mathbf{\theta }})\prod _{j=1}^{J}p({\breve{\mathbf{f}}}_j\mid \breve{\mathbf{u}})p\left(\breve{\mathbf{u}}\right).$$

(9)

For each ${\vartheta }_q(·)$, we introduce pseudo-variables ${\breve{\mathbf{u}}}_q={[{\vartheta }_q\left({\breve{\mathbf{z}}}_1^q\right),\dots ,{\vartheta }_q\left({\breve{\mathbf{z}}}_M^q\right)]}^{\top }\in {\mathbb{R}}^M$, by evaluating ${\vartheta }_q(·)$ at ${\breve{\mathbf{Z}}}_q=[{\breve{\mathbf{z}}}_1^q,\dots ,{\breve{\mathbf{z}}}_M^q]\in {\mathbb{R}}^{M\times P}$. Likewise, $\breve{\mathbf{u}}={[{\breve{\mathbf{u}}}_1^{\top },\dots ,{\breve{\mathbf{u}}}_Q^{\top }]}^{\top }\in {\mathbb{R}}^{QM}$, then:

$$p({\breve{\mathbf{f}}}_j\mid \breve{\mathbf{u}})=\mathcal{N}\left({\breve{\mathbf{f}}}_j|{\mathbf{K}}_{{\breve{\mathbf{f}}}_j\breve{\mathbf{u}}}{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}^{-1}\breve{\mathbf{u}},{\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{f}}}_j}-{\mathbf{K}}_{{\breve{\mathbf{f}}}_j\breve{\mathbf{u}}}{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}^{-1}{\mathbf{K}}_{\breve{\mathbf{u}}{\breve{\mathbf{f}}}_j}\right),$$

(10)

$$p\left(\breve{\mathbf{u}}\right)=\mathcal{N}(\breve{\mathbf{u}}\mid \mathbf{0},{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}})=\prod _{q=1}^Q\mathcal{N}({\breve{\mathbf{u}}}_q\mid \mathbf{0},{\mathbf{K}}_{{\breve{\mathbf{u}}}_q{\breve{\mathbf{u}}}_q});$$

(11)

where ${\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}\in {\mathbb{R}}^{QM\times QM}$ is block-diagonal, with ${\mathbf{K}}_{{\breve{\mathbf{u}}}_q{\breve{\mathbf{u}}}_q}\in {\mathbb{R}}^{M\times M}$ computed from ${\kappa }_q(·,·)$. The covariance matrix ${\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{f}}}_j}\in {\mathbb{R}}^{N\times N}$ has elements ${\sum }_{q=1}^Q{w}_{j,q}^2{\kappa }_q({\mathbf{x}}_n,{\mathbf{x}}_{n^{\prime }})$.

Similarly, ${\mathbf{K}}_{{\breve{\mathbf{f}}}_j\breve{\mathbf{u}}}=[{\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{u}}}_1},\dots ,{\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{u}}}_Q}]\in {\mathbb{R}}^{N\times QM}$, where ${\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{u}}}_q}\in {\mathbb{R}}^{N\times M}$ holds elements $w_{j,q}{\kappa }_q({\mathbf{x}}_n,{\breve{\mathbf{z}}}_m^q)$, with $m\in \{1,2,\dots ,M\}$. Similar to the CGP, the posterior distribution of the CCGP, denoted as $p(\breve{\mathbf{f}},\breve{\mathbf{u}}\mid \mathit{y})$, is generally intractable in closed form. Consequently, it is approximated using a parameterized variational distribution, i.e., $p(\breve{\mathbf{f}},\breve{\mathbf{u}}\mid \mathit{y})\approx \tilde{q}(\breve{\mathbf{f}},\breve{\mathbf{u}})$, following the approach in [23]:

$$\tilde{q}(\breve{\mathbf{f}},\breve{\mathbf{u}})=p(\breve{\mathbf{f}}\mid \breve{\mathbf{u}})\tilde{q}\left(\breve{\mathbf{u}}\right)=\prod _{j=1}^{J}p({\breve{\mathbf{f}}}_j\mid \breve{\mathbf{u}})\prod _{q=1}^Q\tilde{q}\left({\breve{\mathbf{u}}}_q\right),$$

(12)

where $\tilde{q}\left({\breve{\mathbf{u}}}_q\right)=\mathcal{N}({\breve{\mathbf{u}}}_q\mid {\breve{\mu }}_q,{\breve{\mathbf{V}}}_q)$ and $\tilde{q}\left(\breve{\mathbf{u}}\right)=\mathcal{N}(\breve{\mathbf{u}}\mid \breve{\mu },\breve{\mathbf{V}})$. Besides, $\breve{\mu }=[{\breve{\mu }}_1^{\top },\dots ,{\breve{\mu }}_Q^{\top }]\in {\mathbb{R}}^{QM}$ (${\breve{\mu }}_q\in {\mathbb{R}}^M$) and $\breve{\mathbf{V}}\in {\mathbb{R}}^{QM\times QM}$ is a block-diagonal matrix holding covariance blocks ${\breve{\mathbf{V}}}_q\in {\mathbb{R}}^{M\times M}.$ We then approximate the joint posterior over all latent functions and their corresponding inducing variables using a factorized variational distribution. This approximation leads to the following ELBO:

$$\mathcal{L}=\sum _{n=1}^N{\mathbb{E}}_{\tilde{q}\left({\breve{\mathbf{f}}}_1\right),\dots ,\tilde{q}\left({\breve{\mathbf{f}}}_J\right)}\left[logp({\mathit{y}}_n\mid {\theta }_1\left({\mathbf{x}}_n\right),\dots ,{\theta }_J\left({\mathbf{x}}_n\right))\right]-\sum _{q=1}^Q{\mathbb{D}}_{\mathrm{KL}}\left(\tilde{q}\left({\breve{\mathbf{u}}}_q\right)\parallel p\left({\breve{\mathbf{u}}}_q\right)\right),$$

(13)

where:

$$\tilde{q}\left({\breve{\mathbf{f}}}_j\right)=\mathcal{N}\left({\breve{\mathbf{f}}}_j\mid {\mathbf{K}}_{{\breve{\mathbf{f}}}_j\breve{\mathbf{u}}}{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}^{-1}\breve{\mu },{\mathbf{K}}_{{\breve{\mathbf{f}}}_j{\breve{\mathbf{f}}}_j}+{\mathbf{K}}_{{\breve{\mathbf{f}}}_j\breve{\mathbf{u}}}{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}^{-1}(\breve{\mathbf{V}}-{\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}){\mathbf{K}}_{\breve{\mathbf{u}}\breve{\mathbf{u}}}^{-1}{\mathbf{K}}_{\breve{\mathbf{u}}{\breve{\mathbf{f}}}_j}\right).$$

(14)

The first term in equation 13 encourages the variational posterior to explain the observed labels, while the second penalizes deviation from the GP prior over latent factors.

2.3. CCGP-based Crowd Learning and Localized Annotator Trustworthiness

Consider a multi-annotator dataset ${\mathcal{D}}_R=\{\mathbf{X}\in {\mathbb{R}}^{N\times P},\mathbf{Y}\in {\mathbb{R}}^{N\times R}\}$, where $\mathbf{X}$ denotes the input features and $\mathbf{Y}$ is the corresponding label matrix provided by R annotators with unknown and potentially heterogeneous levels of expertise. The entry $y_n^r\in \mathbf{Y}$ represents the label assigned by annotator r to the n-th sample. Since annotators may not label all samples, the matrix $\mathbf{Y}$ can contain missing values. Let $\mathcal{R}\subset \{1,\dots ,N\}\times \{1,\dots ,R\}$ denote the set of observed annotation indices. Then, for each $(n,r)\in \mathcal{R}$, a label $y_n^r$ is observed, while for $(n,r)\notin \mathcal{R}$, the annotation is considered missing and excluded.

Here, we propose a CCGP-based crowd learning framework for multi-annotator regression tasks, referred to as MAR-CCGP. This approach is designed with two primary objectives: (i) to model each annotator’s performance as a localized function of the input space, thereby capturing annotator-specific trustworthiness across different regions; and (ii) to accurately infer the true label $y_{\ast }\in \mathbf{Y}$ for a new input ${\mathbf{x}}_{\ast }\in {\mathbb{R}}^P$. Notably, the method operates in a fully unsupervised manner with respect to annotator reliability, no additional supervision regarding annotator expertise, experience, or consistency is assumed.

Conversely, for real-value outputs, we follow the multi-annotator approach in [49], where each $y_n^r$ is considered as a corrupted version of the estimated hidden ground truth ${\widehat{y}}_n\in \mathbb{R}$, yielding:

$$p(\mathbf{Y}\mid \breve{\mathbf{\theta }})=\prod _{(n,r)\in \mathcal{R}}\mathcal{N}(y_n^r\mid {\widehat{y}}_n,{\breve{v}}_n^r)$$

(15)

where ${\breve{v}}_n^r\in {\mathbb{R}}^{+}$ denotes the error variance associated with the r-th annotator for instance n. In our MAR-CCGP framework, each parameter of the likelihood in $\breve{\theta }$ is linked to an LF ${\breve{f}}_j(·)$, as defined in equation 8. Specifically, our model employs $J=R+1$ LFs: one dedicated to capturing the latent ground truth, and the remaining R to characterize the input-dependent error variances of each annotator:

$$\begin{array}{cc}\hfill {\widehat{y}}_n={\breve{f}}_1\left({\mathbf{x}}_n\right)=& \sum _{q=1}^Q{\omega }_{1,q}{\vartheta }_q\left({\mathbf{x}}_n\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\breve{v}}_n^r=exp\left(f_{l_r}\left({\mathbf{x}}_n\right)\right)=& exp\left(\sum _{q=1}^Q{\omega }_{l_r,q}{\vartheta }_q\left({\mathbf{x}}_n\right)\right);\hfill \end{array}$$

(17)

where $\forall l_r\in \{2,\dots ,J\}$. Note that an exponential transformation is applied to the corresponding LF in equation 17 to ensure that the annotator-specific variance remains strictly positive, i.e., ${\breve{v}}_n^r>0$.

Afterward, an ELBO-based optmization from equation 13 is introduced for our MAR-CCGP, as follows:

$$\mathcal{L}=\sum _{(n,r)\in \mathcal{R}}{\mathbb{E}}_{\tilde{q}\left({\breve{\mathbf{f}}}_1\right),\dots ,\tilde{q}\left({\breve{\mathbf{f}}}_{l_R}\right)}\left[log\left(\mathcal{N}(y_n^r\mid {\widehat{y}}_n,{\breve{v}}_n^r)\right)\right]-\sum _{q=1}^Q{\mathbb{D}}_{\mathrm{KL}}\left(\tilde{q}\left({\breve{\mathbf{u}}}_q\right)\parallel p\left({\breve{\mathbf{u}}}_q\right)\right),$$

(18)

where $\tilde{q}\left({\breve{\mathbf{f}}}_j\right)$ for $j\in \{1,l_1,l_2,\dots ,l_R\}$ denotes the variational marginal over the latent function values, as defined in equation 14; $\tilde{q}\left({\breve{\mathbf{u}}}_q\right)$ represents the variational distribution over the inducing variables (see equation 12); and $p\left({\breve{\mathbf{u}}}_q\right)$ corresponds to the GP prior over the inducing variables, given in equation 11.

In turn, given a new input sample ${\mathbf{x}}_{\ast }$, our goal is to compute the predictive mean and variance for both the estimated ground truth ${\widehat{y}}_{\ast }$ and the corresponding annotator-specific error variances ${\breve{v}}_{\ast }^r$. Specifically, as we defined the ground truth prediction as ${\widehat{y}}_{\ast }={\breve{f}}_1\left(\mathbf{x}\right)$, the posterior distribution over ${\breve{f}}_1\left(\mathbf{x}\right)$ is given by:

$$\tilde{q}\left({\breve{f}}_1\left({\mathbf{x}}_{\ast }\right)\right)=\int p({\breve{f}}_1\left({\mathbf{x}}_{\ast }\right)\mid {\breve{\mathbf{u}}}_q)\tilde{q}\left({\breve{\mathbf{u}}}_q\right)d{\breve{\mathbf{u}}}_q=\mathcal{N}\left({\breve{f}}_1\left({\mathbf{x}}_{\ast }\right)\mid {\breve{\mu }}_{1,\ast },{\sigma }_{1,\ast }^2\right);$$

(19)

where:

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\widehat{y}}_{\ast }\right]=& {\breve{\mu }}_{1,\ast }\in \mathbb{R},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \mathrm{var}\left[{\widehat{y}}_{\ast }\right]=& {\sigma }_{\ast }^2\in {\mathbb{R}}^{+}.\hfill \end{array}$$

(21)

Similarly, due to the exponential transformation in equation 17, the posterior distribution of the annotator-specific variance ${\breve{v}}_{\ast }^r$ follows a log-normal distribution. Its parameters are determined by the predictive mean ${\breve{\mu }}_{l_r,\ast }$ and variance ${\sigma }_{l_r,\ast }^2$ in $\tilde{q}\left({\breve{f}}_{l_r}\left({\mathbf{x}}_{\ast }\right)\right)$, yielding:

$$\begin{array}{cc}\hfill \mathbb{E}\left[{\breve{v}}_{\ast }^r\right]=& exp\left({\breve{\mu }}_{l_r,\ast }+{\displaystyle \frac{{\sigma }_{l_r,\ast }^2}{2}}\right)\in {\mathbb{R}}^{+},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill \mathrm{var}\left[{\breve{v}}_{\ast }^r\right]=& exp\left(2{\breve{\mu }}_{l_r,\ast }+{\sigma }_{l_r,\ast }^2\right)\left(exp\left({\sigma }_{l_r,\ast }^2\right)-1\right)\in {\mathbb{R}}^{+}.\hfill \end{array}$$

(23)

Finally, the proposed MAR-CCGP framework enables the assessment of localized annotator trustworthiness through a probabilistic reliability score derived from the model’s posterior predictions. For each annotator r and sample ${\mathbf{x}}_n$, the trustworthiness score ${\mathrm{T}}_r\left({\mathbf{x}}_n\right)$ is defined as:

$${\mathrm{T}}_r\left({\mathbf{x}}_n\right)=exp\left(-{\displaystyle \frac{|y_n^r-{\widehat{y}}_n{|}_2^2}{2{\left({\breve{v}}_n^r\right)}^{\gamma }}}\right),$$

(24)

where ${\widehat{y}}_n$ denotes the predicted ground truth for input sample ${\mathbf{x}}_n$, and ${\breve{v}}_n^r$ is the estimated annotator-specific variance as in equations 20 and 22, respectively. The exponent $\gamma \in (0,1]$ modulates the sensitivity to the model’s predicted uncertainty, acting as a sublinear scaling factor. Setting $\gamma =0.5$ was found to stabilize the trustworthiness score in regions of low predicted variance, preventing numerical explosions and yielding better alignment with empirical annotator behavior. This formulation corresponds to the scaled likelihood of the observed annotation $y_n^r$ under the model’s uncertainty, offering a principled, data-dependent trust metric for each annotator at a given input location.

3. Experimental Set-Up

To evaluate the effectiveness of the proposed MAR-CCGP framework in modeling information from multiple annotators, we conduct a series of experiments using both semi-synthetic and real-world datasets. The semi-synthetic benchmarks provide a controlled environment with access to ground-truth labels, enabling rigorous assessment of the model’s ability to infer annotator trustworthiness. In contrast, the real-world dataset—comprising sensory evaluations of cacao products (see the LUKER-CACAO dataset description in Section 2.1)—lacks ground-truth annotations, making it an ideal scenario to assess the framework’s capabilities in consensus estimation and uncertainty quantification.

3.1. Semi-Synthetic Datasets Annotation Simulation

We evaluate our MAR-CCGP approach under controlled conditions employing regression datasets from the University of California Irvine (UCI) machine learning repository(see https://archive.ics.uci.edu/datasets and Table 4). Each dataset provides continuous targets and input features, enabling the generation of simulated multi-annotator labels with structured noise profiles.

Moreover, to simulate annotator-specific variability, we adopt the following scheme:

-

The input data are standardized and then projected into a two-dimensional space using Uniform Manifold Approximation and Projection (UMAP) [50], to preserve local data structure by minimizing the cross-entropy between high-dimensional and low-dimensional fuzzy simplicial representations.

-

A K-means algorithm with C clusters is then applied to the UMAP projection to derive pseudo-contextual input space partitions [30], as latent indicators of instance-specific difficulty or domain shifts, used to modulate annotator behavior.

-

Let $c_n\in \{1,\dots ,C\}$ denote the cluster assignment associated with the input ${\mathbf{x}}_n$, and let $y_n$ represent the corresponding ground-truth regression value. The simulated label assigned by annotator r to instance n is then generated as follows:

$$y_n^r=y_n+{ϵ}_n^r,{ϵ}_n^r\sim \mathcal{N}({ϵ}_n^r\mid 0,{\tilde{\sigma }}_{r,c_n}^2);$$

(25)

where ${\tilde{\sigma }}_{r,c_n}^2\in {\mathbb{R}}^{+}$ is the annotator-specific error variance for cluster $c_n$. To ensure interpretability and reproducibility, each annotator’s variance profile across clusters is defined in terms of a fixed Signal-to-Noise Ratio (SNR) in dB. Then, the multiple annotator target matrix for regression tasks $\mathbf{Y}\in {\mathbb{R}}^{N\times R}$, is built as in equation 25.

By varying the SNR across annotators and data clusters, we construct distinct variance profiles wherein certain annotators exhibit high reliability in specific regions of the input space and low reliability in others. This setup is intended to emulate real-world, non-stationary labeling behavior. For concreteness in our semi-synthetic data experiments, we fix the number of clusters to $C=4$ and the number of annotators to $R=5$. The SNR values are varied according to the profiles shown in Figure 5. One annotator is designated as an expert, consistently exhibiting an SNR of $10\mathrm{dB}$ across all data clusters.

3.2. Quality Assessment, Method Comparison, and Training Details

To evaluate whether the model-derived reliability scores accurately reflect annotator behavior, we aggregate both the predicted trustworthiness scores and empirical performance metrics across annotators and clusters. Specifically, for each annotator, we compute the average trustworthiness score within each cluster (see equation 24), defined as:

$${\overline{\mathcal{T}}}_r\left({\Omega }_{\tilde{c}}\right)={\displaystyle \frac{1}{|{\Omega }_{\tilde{c}}|}}\sum _{{\mathbf{x}}_n\in {\Omega }_{\tilde{c}}}{\mathcal{T}}_r\left({\mathbf{x}}_n\right),$$

(26)

where ${\Omega }_{\tilde{c}}$ denotes the set of input instances in cluster $\tilde{c}\in C$. For comparison, we compute the empirical coefficient of determination between annotator labels and model-predicted ground truth in the same regions:

$$R_r^2\left({\Omega }_{\tilde{c}}\right)=1-{\displaystyle \frac{{\sum }_{\forall n:{\mathbf{x}}_n\in {\Omega }_{\tilde{c}}}{\parallel y_n^r-y_n\parallel }_2^2}{{\sum }_{\forall n:{\mathbf{x}}_n\in {\Omega }_{\tilde{c}}}{\parallel y_n^r-\overline{y}\parallel }_2^2}},$$

(27)

where $\overline{y}=\frac{1}{|{\Omega }_{\tilde{c}}|}{\sum }_n{y}_n,$∀$n:{\mathbf{x}}_n\in {\Omega }_{\tilde{c}}.$

We benchmarked the proposed MAR-CCGP against three representative baselines, each embodying a distinct strategy for handling multi-annotator learning. As outlined in Table 5, the GPR-GT model serves as an optimistic upper bound, leveraging a GP regressor [30] trained directly on ground truth targets. This approach bypasses the inherent complexities of annotator-induced variability, thereby providing a best-case performance scenario. In contrast, the GPR-AVG represents a naive lower-bound reference, where annotations are averaged per instance to form pseudo-targets for GP training [23]. This method assumes uniform annotator reliability and fails to account for noise or systematic biases. A more nuanced alternative is offered by LKAAR, which simultaneously models annotator-specific noise and bias while learning an input-dependent reliability function via kernel alignment [31]. Notably, both LKAAR and our MAR-CCGP are among the few methods capable of modeling annotator consistency as a function of the input space.

Also, to quantitatively assess the predictive performance of the MAR-CCGP model for the semi-synthetic datasets, we include three standard regression measures: Mean Squared Error (MSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE):

$$\begin{array}{cc}\hfill \mathrm{MSE}(\mathit{y},\widehat{\mathit{y}})& ={\displaystyle \frac{1}{N}}\sum _{n=1}^{N_{\mathrm{test}}}{(y_n-{\widehat{y}}_n)}^2,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \mathrm{MAE}(\mathit{y},\widehat{\mathit{y}})& ={\displaystyle \frac{1}{N}}\sum _{n=1}^{N_{\mathrm{test}}}|y_n-{\widehat{y}}_n|,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill \mathrm{MAPE}(\mathit{y},\widehat{\mathit{y}})& ={\displaystyle \frac{100}{N}}\sum _{n=1}^{N_{\mathrm{test}}}\left|{\displaystyle \frac{y_n-{\widehat{y}}_n}{y_n+\breve{ϵ}}}\right|,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill R^2(\mathit{y},\widehat{\mathit{y}})& =1-{\displaystyle \frac{{\sum }_{n=1}^N{(y_n-{\widehat{y}}_n)}^2}{{\sum }_{n=1}^{N_{\mathrm{test}}}{(y_n-\overline{y})}^2}},\hfill \end{array}$$

(31)

where $y_n\in \mathit{y}$ is the true target (available in semi-synthetic scenarios) and ${\widehat{y}}_n\in \widehat{\mathit{y}}$ is the predicted ground truth. $\breve{ϵ}>0$ is a small constant added to avoid division by zero when $y_n\approx 0$ (set to $\breve{ϵ}={10}^{-6}$ in our experiments), $\overline{y}={\displaystyle \frac{1}{N}}{\sum }_{n=1}^N{y}_n$ denotes the sample mean of the ground truth values.

Table 5. Comparative summary of baseline and advanced models for multi-annotator regression.

Method	Acronym	Description
Gaussian Process on Ground Truth [30]	GPR-GT	Supervised GP regression trained on true outputs. Serves as an oracle upper bound and assumes full access to ground truth. Does not model annotators or uncertainty.
Gaussian Process on Average Annotations [23]	GPR-AVG	GP trained on the per-instance average of annotator targets. Assumes annotators are unbiased and neglects individual reliability. Serves as a baseline that models only the consensus.
Localized Kernel Alignment-based Annotator Relevance [31]	LKAAR	Jointly estimates annotator bias and variance, and embeds annotator consistency as a kernelized function over the input space. Provides localized reliability estimates.
Multi-Annotator Regression based on Correlated Chained Gaussian Process (ours)	MAR-CCGP	Proposed model. Captures latent inter-annotator correlations and input-dependent noise via correlated-chained latent functions and sparse variational GPs. Produces localized consistency-trustworthiness estimates.

Table 5. Comparative summary of baseline and advanced models for multi-annotator regression.

Method	Acronym	Description
Gaussian Process on Ground Truth [30]	GPR-GT	Supervised GP regression trained on true outputs. Serves as an oracle upper bound and assumes full access to ground truth. Does not model annotators or uncertainty.
Gaussian Process on Average Annotations [23]	GPR-AVG	GP trained on the per-instance average of annotator targets. Assumes annotators are unbiased and neglects individual reliability. Serves as a baseline that models only the consensus.
Localized Kernel Alignment-based Annotator Relevance [31]	LKAAR	Jointly estimates annotator bias and variance, and embeds annotator consistency as a kernelized function over the input space. Provides localized reliability estimates.
Multi-Annotator Regression based on Correlated Chained Gaussian Process (ours)	MAR-CCGP	Proposed model. Captures latent inter-annotator correlations and input-dependent noise via correlated-chained latent functions and sparse variational GPs. Produces localized consistency-trustworthiness estimates.

All experiments were conducted on the Kaggle platform using GPU-enabled notebooks. Each session provided access to a single NVIDIA Tesla P100 GPU (16GB VRAM), 30 GB of RAM, and 4 vCPUs from an Intel Xeon CPU @ 2.20GHz. The code-base was developed in Python 3.11.11, and all models were executed in the default Kaggle environment, with package versions listed in the notebook dependencies. Our MAR-CCGP model was implemented using GPflow 2.10.0 atop TensorFlow 2.18.0. Dimensionality reduction and clustering were handled by cuML 25.2.1 (RAPIDS), while preprocessing relied on scikit-learn 1.2.1. Additional tools included NumPy, Pandas, and Matplotlib for analysis and visualization.

All GP-based approaches and the regression components within LKAAR—used sparse variational inference with $M=200$ inducing points and a mini-batch size of 128. Inducing inputs were initialized via k-means clustering. A squared exponential kernel was used in all cases. Training was carried out using the Adam optimizer with a learning rate of ${10}^{-2}$ for up to 1000 steps. If the ELBO did not improve for 20 consecutive iterations, the learning rate was halved to a minimum of ${10}^{-6}$. Early stopping was triggered if no improvement was observed over 500 iterations.

Each experiment was repeated over 15 randomized train/test splits (70/30). For the semi-synthetic datasets, splits were stratified by the UMAP-based cluster labels to preserve input-space diversity across folds. In the cacao-based real-world dataset (LUKER-CACAO), where no cluster labels were available, standard random splits were applied. Random seeds were fixed for all runs to ensure reproducibility. Code and full experiment notebooks for semi-synthetics datasets are publicly available at: https://github.com/UN-GCPDS/python-gcpds.luker_multiple_annotators. Access to the LUKER-CACAO dataset is restricted owing to copyright limitations established by the Casa Luker organization.

Figure 7. LUKER-CACAO dataset MAR-CCGP experimental set-up.

Figure 7. LUKER-CACAO dataset MAR-CCGP experimental set-up.

4. Results and Discussion

4.1. Semi-Synthetic Datasets Results

We first evaluate all models using the semisynthetic benchmark datasets described in Table 4. These datasets allow controlled experimentation with known ground truth, enabling a detailed assessment of both predictive accuracy and the quality of learned annotator reliability.

To simulate non-stationary annotator behavior, we followed the procedure outlined in Section 3.1. Ground truth labels were corrupted using structured noise profiles defined by a SNR matrix, assigning fixed reliability levels per annotator and input cluster. As shown in Figure 5, Annotator 1 is modeled as a consistent expert (SNR = 10 dB across all clusters), while Annotators 2–5 exhibit cluster-dependent behavior. For instance, Annotator 3 is highly accurate in Cluster 2 (6 dB) but unreliable in Cluster 1 (–3 dB), simulating realistic local variability in annotator performance. Figure 8 illustrates the UMAP-based clustering for the Bike Sharing dataset, and Figure 9 shows the resulting noisy annotations alongside the ground truth.

All models were evaluated over 15 randomized 70/30 train/test splits, stratified by the clusters used in the annotator simulation. We report the mean and standard deviation of test-set metrics: MSE, MAE, MAPE, and $R^2$. Results are presented in Table 6.

Across all datasets, MAR-CCGP consistently achieves the best performance in terms of both error metrics (MSE, MAE, MAPE) and explanatory power ($R^2$), highlighting its ability to model input-dependent annotator reliability and inter-annotator correlations. This advantage stems from its capacity to learn both latent ground truth and annotator noise structure jointly, leveraging the semi-parametric latent factor model and shared latent processes. GPR-GT, trained using uncorrupted ground truth, serves as a strong upper bound, yet MAR-CCGP often matches or approaches its performance—particularly in more complex datasets—despite only observing noisy annotations. In contrast, GPR-AVG, which averages noisy labels, suffers from bias introduced by context-dependent errors, yielding significantly higher error rates. LKAAR, while modeling local annotator reliability, fails to capture dependencies between annotators or fully account for the structured noise patterns, leading to lower predictive accuracy and less stable performance. These results support the hypothesis that modeling non-stationary, correlated annotation noise yields more robust predictors, even under limited and noisy supervision.

To evaluate MAR-CCGP’s ability to recover localized annotator trustworthiness, we applied the assessment framework described in Section 2.3. Specifically, for each dataset, we partitioned the input space into clusters and computed two cluster-level matrices: (1) the empirical coefficient of determination $R_r^2\left({\Omega }_{\tilde{c}}\right)$, which quantifies the agreement between each annotator’s labels and the inferred ground truth, and (2) the average localized trustworthiness score ${\overline{\mathcal{T}}}_r(\Omega \tilde{c})$ derived from the model’s posterior predictive distribution. These metrics enable a direct comparison between empirical annotator behavior and the model’s estimated reliability across different regions of the input space.

Figure 10 presents the cluster-wise empirical $R_r^2$ scores (left) alongside the corresponding MAR-CCGP trustworthiness estimates (right) for the Bike Sharing dataset. The left panel clearly shows that Annotator 1 is the most reliable across all clusters, achieving the highest $R_r^2$ in every region — as expected given its uniform SNR profile. MAR-CCGP’s estimated trustworthiness on the right reproduces this behavior closely, assigning consistently high trust scores to Annotator 1 across all clusters. In contrast, Annotators 2–5 exhibit substantial local variations: for instance, Annotator 2 is most accurate in Cluster 1 but becomes much less accurate in Clusters 3 and 4, and this is accurately reflected in its estimated trustworthiness profile. Similarly, Annotator 5 is most reliable in Cluster 4 ($R_r^2=0.360\pm 0.017$), which matches the model’s highest trustworthiness estimate ($0.805\pm 0.007$) in that region. Overall, these heatmaps reveal a strong correspondence between empirical annotator performance and the model’s estimated reliability.

To quantify this alignment, we computed the squared Pearson correlation between the flattened empirical and model-derived trustworthiness matrices, obtaining a strong value of $0.922$. This indicates that the model-internal trust estimates not only reflect qualitative trends but also align quantitatively with observed annotator behavior. Figure 11 further complements this analysis by visualizing per-sample trustworthiness across the input space, illustrating spatially localized patterns consistent with the annotator-specific noise profiles imposed during simulation.

To summarize these findings across all datasets, we averaged the cluster-wise scores per annotator, yielding the aggregated heatmaps shown in Figure 12. Annotator 1 emerges as consistently the most trustworthy across all datasets, exactly as simulated by its uniform SNR profile. This trend is well recovered by the MAR-CCGP model, which also correctly reflects the relative strengths and weaknesses of the remaining annotators. To quantify this overall alignment, we computed the squared Pearson correlation between the flattened empirical and model-estimated trustworthiness scores across all datasets, obtaining a high value of $0.853$. This strong correspondence underlines the proposed method’s effectiveness in capturing true annotator behaviors even under diverse, noisy conditions, and confirms that MAR-CCGP provides robust, input-dependent trustworthiness estimates that closely match empirical annotator reliability.

4.2. LUKER-CACAO Dataset Results

Having demonstrated the predictive performance and trustworthiness inference capabilities of MAR-CCGP on semi-synthetic data, we now turn to a real-world application using the LUKER-CACAO dataset described in Section 2.1. Unlike the controlled benchmarks, this sensory evaluation dataset contains separate, curated input–output tables for each of eight sensory attributes (acidity, bitterness, aroma, astringency, sweetness, hardness, melting speed, and global impression), along with physicochemical features measured per sample (moisture, fat content, granulometry, plastic viscosity, yield stress). Importantly, the dataset contains ratings from five annotators without ground-truth labels, making this an ideal setting for MAR-CCGP to estimate a consensus profile and explore localized annotator trustworthiness.

Following the same experimental setup as the semi-synthetic case (training and testing over 15 random repetitions with 70/30 splits), we applied MAR-CCGP to each sensory variable separately. This setup enables a fair and consistent comparison across attributes and ensures that observed trends reflect annotator behavior rather than sampling variability.

Figure 13 shows the estimated mean sweetness profile learned by the model across all samples, along with its 95% predictive credible intervals. Individual annotator ratings are overlaid for direct comparison. The model’s mean estimate lies centrally among the noisy annotator ratings, especially where agreement is high, and predictive intervals widen in the left-hand region of the plot where annotators strongly disagree. Conversely, intervals narrow as ratings converge at higher sweetness levels, demonstrating the model’s ability to represent both consensus and uncertainty appropriately.

To further explore annotator-specific behavior, we investigated correlations between physicochemical properties and MAR-CCGP–inferred trustworthiness scores. Figure 14 presents a bipartite graph of these relationships. Nodes on the left represent the five physicochemical features, while nodes on the right represent the five annotators. Edge color indicates the Pearson correlation between a feature’s value and an annotator’s estimated trustworthiness across all samples. Warm edges show positive correlations and cool edges show negative correlations. The graph reveals that several annotators—most notably 154 and 160—exhibit negative correlations with fat content and moisture, suggesting that as these features increase, sweetness ratings become less reliable for these annotators. This is consistent with the intuition that fat and moisture may mask sweetness perception. Conversely, granulometry and yield stress show weaker or even positive correlations for some annotators (e.g. Annotator 135), implying that these panelists evaluate sweetness more reliably under these physical conditions.

To visualize trustworthiness across the physicochemical input space, Figure 15 depicts UMAP projections of the cacao samples, color-coded by MAR-CCGP’s estimated trustworthiness per annotator. The top-left panel depicts the estimated mean sweetness profile projected into the UMAP space, serving as a baseline. The remaining panels show each annotator’s trust scores across this space. Annotators 154 and 155 achieve high trustworthiness across most regions (consistently warm yellow hues), indicating stable and reliable performance. Annotators 160 and 179, by contrast, display region-specific drops in trustworthiness, especially for clusters of samples with high moisture or fat content. These localized patterns highlight systematic variations in annotator behavior across product types.

Finally, Table 7 and Figure 16 provide a detailed summary of estimated annotator trustworthiness across all sensory variables and product types. Consistent with the qualitative UMAP patterns, Annotators 154 and 155 stand out as the most consistently trustworthy panelists, yielding mean trust scores above $0.80$ for most attributes with low standard deviations across different product types. This is evident not only in sweetness but also for aroma, melting speed, hardness, and global impression, where both annotators maintain strong and stable reliability. Annotator 160, by contrast, exhibits markedly lower mean trustworthiness across many attributes (typically in the $0.60$–$0.70$ range), with larger standard deviations — especially for perceptually challenging attributes like astringency and bitterness — indicating highly variable performance that depends on the product’s physicochemical profile. Annotator 179 also shows a similar pattern of reduced consistency, attaining good trust scores for a few localized product types but generally underperforming relative to the most consistent annotators. Interestingly, some sensory variables, such as bitterness and astringency, present inherently greater variability in trust scores across all panelists, suggesting that these attributes may be more difficult to evaluate reliably across the diverse product set. In contrast, more mechanically defined properties like hardness and melting speed exhibit higher and more uniform trustworthiness estimates across all annotators, which could reflect the greater perceptual salience of these sensory dimensions. Taken together, these findings demonstrate that MAR-CCGP successfully disentangles systematic annotator-specific behavior from product-specific effects, yielding granular trustworthiness estimates that highlight where particular panelists excel or struggle. This level of insight can inform sensory panel calibration, personalized training, and targeted quality control by drawing attention to specific attributes and product types that require additional guidance or closer monitoring.

Together, these findings demonstrate MAR-CCGP’s capacity to extract nuanced, input-dependent trustworthiness estimates even in real-world sensory settings without ground-truth supervision. The model’s trust scores closely align with expected behavior derived from physicochemical features, and its uncertainty estimates highlight regions where panel agreement is weakest. Importantly, these insights offer actionable guidance for panelist training and calibration, allowing sensory scientists to identify panelists who struggle with specific product types and potentially target these areas for further calibration or retraining.

4.3. Limitations

While the MAR-CCGP framework provides a robust solution for integrating multi-annotator sensory data with physicochemical profiles, several limitations remain. First, the model assumes that annotator reliability can be effectively captured through input-dependent variance, which may oversimplify annotators’ behavior in scenarios involving complex biases or strategic labeling. Second, the proposed trustworthiness score relies on probabilistic estimations derived from model outputs, but its interpretation may be challenging in the absence of external validation mechanisms or behavioral ground truth. Moreover, although the method performs well in semi-synthetic and real-world settings, it is currently limited to low-dimensional input spaces with a manageable number of annotators; scalability to high-dimensional sensory domains or crowd-scale scenarios may require further optimization or approximation strategies. Additionally, the framework depends on Gaussian process inference, which incurs high computational costs as the dataset size increases. Finally, while the LUKER-CACAO dataset enables a valuable real-world demonstration, its proprietary nature restricts reproducibility and broader benchmarking by the community.

5. Conclusions

We introduce a Multi-Annotator Regression framework based on Correlated Chained Gaussian Processes, named MAR-CCGP, a novel multi-annotator approach that jointly models continuous sensory scores and input-dependent annotator reliability through a probabilistic Gaussian Process formulation. The key conceptual innovation lies in disentangling true perceptual signals from annotator-specific noise using a latent consensus function and localized trust estimation. Unlike previous methods that assume uniform or global annotator performance, MAR-CCGP learns region-specific annotator trustworthiness scores, enabling interpretability and robustness in the face of sparse, noisy, and subjective supervision. The model is particularly suited for domains like food science, where expert annotations are limited and often exhibit contextual biases.

Our results demonstrate the effectiveness of MAR-CCGP across both real-world and semi-synthetic settings. On the proprietary LUKER-CACAO dataset, the model achieved strong predictive performance for multiple sensory attributes and provided meaningful trust scores that aligned with empirical annotator behavior, especially in dimensions such as bitterness and aroma. In controlled experiments with structured SNR profiles, MAR-CCGP outperformed consensus-only and local-weighting baselines in both RMSE and reliability estimation. Importantly, the model’s ability to recover annotator-specific patterns in low-trust clusters highlights its utility in curating more reliable datasets and guiding future annotation efforts. These findings suggest that MAR-CCGP not only enhances regression accuracy under label noise but also supports informed decision-making through interpretable reliability scores.

Future research could explore scaling MAR-CCGP to high-dimensional input domains and larger annotator pools by integrating sparse approximations or deep kernel learning methods [51,52]. Another promising direction involves extending the framework to heteroscedastic multi-output settings, allowing the simultaneous modeling of correlations between multiple sensory targets and annotator behaviors. Incorporating behavioral signals or auxiliary metadata from annotators—such as labeling time, confidence, or experience level—could further refine trust estimation [53,54]. Lastly, validating the trust scores through longitudinal studies or expert-in-the-loop experiments may reinforce their adoption in high-stakes domains such as medical diagnostics, environmental monitoring, and consumer preference modeling.