Reliable food safety risk assessment via calibrated Mondrian conformal prediction: a machine learning framework with uncertainty quantification

1 Introduction

Food safety represents a critical public health concern with substantial economic implications across global supply chains. The World Health Organization (WHO) estimates that foodborne diseases affect approximately 600 million people annually, resulting in 420,000 deaths and the loss of 33 million disability-adjusted life years (World Health Organization, 2015). Regulatory frameworks, such as the European Union (EU) Rapid Alert System for Food and Feed (RASFF), have been established to facilitate rapid information exchange regarding food safety hazards among member states (Bouzembrak et al., 2018). The RASFF database contains detailed records of food safety notifications, including hazard types, origin countries, product categories, and risk severity classifications, providing a rich data source for predictive analytics (Marvin et al., 2016).

Formally, this study addresses a binary classification problem: predicting whether an RASFF notification is serious (positive class) or not (negative class). The classification of food safety risk severity presents several methodological challenges that limit the applicability of conventional machine learning approaches. The inherent class imbalance in food safety data, in which serious risk cases typically dominate the notification landscape, leads to biased model predictions and inflated accuracy metrics (He and Garcia, 2009). Furthermore, standard classification models produce point predictions without quantifying prediction uncertainty, which is essential for regulatory decision-making, as false negatives may have severe public health consequences (Angelopoulos and Bates, 2023). Additionally, models trained with class balancing techniques often produce miscalibrated probability estimates, undermining the reliability of risk scores for prioritization purposes (Guo et al., 2017).

The motivation for the present work stems from the observation that existing food safety prediction systems lack principled uncertainty quantification mechanisms. While several studies have applied machine learning to RASFF data for hazard prediction and risk assessment (Bouzembrak and Marvin, 2019; Blauw et al., 2021), these approaches provide deterministic predictions without coverage guarantees. Conformal prediction offers a distribution-free framework for constructing prediction sets with finite-sample validity guarantees (Vovk et al., 2005), making it particularly suitable for high-stakes regulatory applications. The Mondrian variant extends these guarantees to conditional coverage within predefined groups (Vovk, 2003), addressing potential disparities in prediction reliability across different geographic regions or hazard categories.

This study introduces a comprehensive framework that synergistically combines gradient boosting classification, post hoc probability calibration, and Mondrian conformal prediction for food safety risk assessment. The framework employs eXtreme Gradient Boosting (XGBoost) as the base classifier with a synthetic minority oversampling technique (SMOTE)-based resampling and class weighting to address imbalance. Platt scaling is applied post hoc to correct probability miscalibration induced by class balancing. Finally, Mondrian conformal prediction ensures group-conditional coverage across geographic regions. Rather than assuming equitable reliability, this study empirically assesses group-wise coverage and characterizes its finite-sample variability. The framework additionally incorporates counterfactual explanation generation to identify actionable modifications that would alter risk classifications.

1.1 Contributions

The major contributions of this work are summarized as follows:

• Integrated framework: This study introduces a novel framework that combines XGBoost classification with SMOTE resampling, Platt scaling calibration, and Mondrian conformal prediction for food safety risk assessment, achieving 90.3% coverage with an average prediction set size of 1.216.
• Post hoc calibration: A Platt scaling module is proposed that reduces the expected calibration error (ECE) by 84% (from 0.267 to 0.043), addressing the probability miscalibration caused by class-balanced training.
• Group-conditional coverage: The Mondrian conformal prediction component provides region-specific coverage guarantees, with four of the five geographic regions achieving the target coverage threshold of 90%.
• Counterfactual explanations: Multimethod counterfactual generation (prototype interpolation, random perturbation, and growing spheres) provides interpretable insights into risk factor modifications.
• Comprehensive evaluation: Extensive experiments on 22,643 RASFF notifications demonstrate the framework’s effectiveness with a Cohen’s kappa of 0.296, a Matthews correlation coefficient (MCC) of 0.297, and an area under the receiver operating characteristic (ROC) curve (AUC-ROC) of 0.763.

The remainder of this article is organized as follows: Section 2 presents a comprehensive review of related work and identifies research gaps. Section 3 describes the proposed framework in detail. Section 4 outlines the experimental setup, including dataset description, baseline methods, and evaluation metrics. Section 5 presents and analyzes the experimental results. Section 6 discusses the findings and their implications. Finally, Section 7 concludes the article and suggests directions for future research.

2 Related work

2.1 Machine learning for food safety

Machine learning applications in food safety have evolved significantly over the past decade, driven by increasing data availability and computational capabilities. Early work focused on applying traditional classifiers to hazard identification tasks. For example, Marvin et al. (2016) demonstrated the potential of Bayesian networks for predicting food safety hazards using RASFF data, achieving moderate accuracy but limited scalability. Bouzembrak et al. (2018) extended this work by incorporating text mining techniques to extract additional features from notification descriptions, improving the classification performance of microbial hazards.

Ensemble methods have shown particular promise for food safety applications due to their robustness to noise and ability to capture complex feature interactions. Bouzembrak and Marvin (2019) compared random forests, gradient boosting, and support vector machines for predicting border rejections and found that ensemble methods consistently outperformed single classifiers. Blauw et al. (2021) applied gradient boosting to the temporal prediction of food safety alerts, incorporating seasonal patterns and origin country risk profiles. More recently, Nogales et al. (2022) explored deep learning architectures for multihazard classification, although the interpretability limitations of neural networks remain a concern for regulatory applications.

2.2 Class imbalance in food safety data

The challenge of class imbalance in food safety data has received considerable attention. Chawla et al. (2002) introduced SMOTE, which has become the de facto standard for addressing minority-class underrepresentation. He and Garcia (2009) provided a comprehensive survey of imbalanced learning techniques, noting that cost-sensitive learning and resampling methods offer complementary benefits. Recent work by Fernández et al. (2018) evaluated SMOTE variants across multiple domains, recommending combined approaches that integrate resampling with algorithmic modifications, such as class weighting.

2.3 Conformal prediction

Conformal prediction, introduced by Vovk et al. (2005), provides a distribution-free framework for constructing prediction sets with guaranteed coverage probabilities. Unlike Bayesian approaches, which require prior specifications, or frequentist confidence intervals, which rely on asymptotic normality, conformal prediction offers finite-sample validity under the sole assumption of data exchangeability (Shafer and Vovk, 2008). This property makes conformal methods particularly attractive for high-stakes applications where coverage guarantees are essential.

The Mondrian variant of conformal prediction extends marginal coverage guarantees to conditional coverage within predefined groups (Vovk, 2003). This extension targets conditional coverage across predefined subpopulations, a requirement increasingly emphasized in fair machine learning (Romano et al., 2020). Importantly, conditional coverage is a finite-sample property, whose realized value in any single split is subject to sampling variability, particularly for small groups (Barber et al., 2021). Barber et al. (2021) characterized the fundamental limits of conditional coverage, showing that non-trivial conditional guarantees require either structural assumptions or group-wise calibration.

Recent advances have focused on improving the efficiency of conformal prediction sets while maintaining coverage guarantees. Romano et al. (2019) proposed conformalized quantile regression for regression tasks, achieving tighter prediction intervals. Angelopoulos et al. (2021) introduced methods for learning optimal conformal predictors that minimize the expected set size. Stutz et al. (2022) demonstrated that training classifiers with conformal-aware objectives can improve downstream conformal prediction efficiency.

Applications of conformal prediction span diverse domains, including medical diagnosis (Lu et al., 2022), drug discovery (Eklund et al., 2015), and autonomous systems (Lindemann et al., 2023). However, applications to food safety remain limited, with the majority of existing work focusing on point predictions without uncertainty quantification.

2.4 Probability calibration

Well-calibrated probability estimates are essential for decision-making under uncertainty; however, modern neural networks and boosted classifiers often produce miscalibrated predictions (Guo et al., 2017). Calibration refers to the alignment between predicted probabilities and empirical frequencies: a perfectly calibrated model predicting 70% probability should be correct 70% of the time across all such predictions.

Platt scaling, introduced by Platt (1999), fits a sigmoid function to transform classifier outputs into calibrated probabilities. Despite its simplicity, Platt scaling remains effective for binary classification and is particularly well suited for correcting systematic over- or under-confidence (Niculescu-Mizil and Caruana, 2005). Zadrozny and Elkan (2002a) proposed histogram binning as a non-parametric alternative, while isotonic regression offers a monotonicity-preserving middle ground (Zadrozny and Elkan, 2002b).

Temperature scaling, a single-parameter variant of Platt scaling, has gained popularity for multiclass calibration (Guo et al., 2017). More recent work has continued to refine post hoc calibration, analyzing the factors that drive miscalibration and proposing targeted scaling schemes. For example, Wang et al. (2023) showed that separating correctly and incorrectly predicted instances during scaling improves calibration. These findings reinforce the idea that the appropriate calibration method is model- and setting-dependent and should be selected empirically for the base classifier at hand. More sophisticated approaches include focal loss calibration (Mukhoti et al., 2020), which jointly optimizes accuracy and calibration during training. Kumar et al. (2019) introduced differentiable calibration objectives that enable end-to-end calibration-aware training.

However, the interaction between class imbalance handling and calibration has received limited attention. Class weighting and oversampling techniques, while improving balanced accuracy, can induce probability miscalibration by distorting the base rate (Wallace and Dahabreh, 2012). This observation motivates the post hoc calibration approach adopted in the present work.

2.5 Counterfactual explanations

Counterfactual explanations identify the minimal input modifications that would alter a model’s prediction, providing actionable insights for end users (Wachter et al., 2018). Unlike feature attribution methods that quantify input importance, counterfactuals answer “what if” questions that directly support decision-making and recourse (Karimi et al., 2022).

Several algorithms have been proposed for counterfactual generation. Mothilal et al. (2020) introduced the diversified counterfactual explanation (DiCE) method, which generates diverse counterfactuals through optimization with diversity constraints. Poyiadzi et al. (2020) proposed a feasible and actionable counterfactual explanation (FACE) method that ensures counterfactual feasibility by restricting paths to observed data manifolds. Growing sphere algorithms (Laugel et al., 2018) progressively expand search radii until valid counterfactuals are found, offering simplicity at the cost of potential suboptimality.

Prototype-based methods leverage class representatives to guide the counterfactual search. Van Looveren and Klaise (2021) combined autoencoders with prototype loss to generate plausible counterfactuals. The present work adapts prototype interpolation for food safety, using class means as target prototypes and binary search for minimal perturbation identification.

2.6 Summary of literature review

Table 1 presents a comprehensive summary of the existing literature, highlighting the methods employed, their advantages, and their limitations.

TABLE 1. Comprehensive literature review summary (2018–2025).

2.7 Identified research gaps

Based on the comprehensive literature review presented in Table 1, several critical research gaps were identified. Table 2 summarizes these gaps and articulates how this study addresses each limitation.

TABLE 2. Summary of research gaps and how the present work addresses them.

2.8 Position of the present work

The preceding literature review reveals several recurring themes and limitations across existing approaches. While significant progress has been made in applying machine learning to food safety prediction, the majority of existing methods suffer from one or more of the gaps identified in Table 2. Specifically, prior work predominantly focuses on improving classification accuracy without addressing the critical need for uncertainty quantification in regulatory decision-making, and techniques for handling class imbalance often induce probability miscalibration, which undermines the reliability of risk scores.

The present study distinguishes itself from prior research in several fundamental ways. Unlike existing food safety prediction systems that provide deterministic point predictions (Bouzembrak and Marvin, 2019; Blauw et al., 2021), the proposed framework delivers prediction sets with statistically guaranteed coverage. Furthermore, while previous applications of conformal prediction have primarily targeted marginal coverage (Shafer and Vovk, 2008), this article employs Mondrian conformal prediction to ensure equitable prediction reliability across geographic regions. The integration of post hoc Platt scaling specifically addresses the calibration distortion induced by SMOTE and class weighting, a challenge that has received limited attention in the conformal prediction literature.

The unique positioning of this work lies in its novel integration of three complementary components: class imbalance handling, probability calibration, and group-conditional conformal prediction. To the best of the authors’ knowledge, Mondrian conformal prediction has not previously been applied to RASFF-based food safety classification by severity with the explicit handling of calibration–imbalance interactions; the contribution is framed as a novel integration rather than an absolute first, and Section 2 positions it against the closest prior work in food safety artificial intelligence, RASFF prediction, and uncertainty-aware learning. This approach not only addresses gaps G1 through G5 (Table 2) but also provides a template for applying uncertainty-aware machine learning in other regulatory domains.

3 Methodology

3.1 Problem formulation

The food safety risk classification task is formalized as follows: Let D = { ( x i , y i ) } i = 1 N denote a dataset of N food safety notifications, where x i ∈ R d represents the d-dimensional feature vector for notification i and y i ∈ { 0,1 } denotes the binary risk label (0 = not serious, 1 = serious). Each notification is additionally associated with a group membership g i ∈ G, where G represents the set of geographic regions.

The objective is to construct a prediction system that produces, for a new notification x, a prediction set C ( x ) ⊆ { 0,1 } satisfying the coverage guarantee (Equation 1):where α ∈ ( 0,1 ) is a user-specified miscoverage rate. For Mondrian conformal prediction, this guarantee is strengthened to hold conditionally within each group (Equation 2):

P y ∈ C x ≥ 1 − α (1)

P y ∈ C x ∣ g = g ′ ≥ 1 − α , ∀ g ′ ∈ G (2)

Additionally, the system should produce calibrated probability estimates p ̂ ( x ) = P ( y = 1 ∣ x ) such that the empirical frequency of positive outcomes matches the predicted probability across all probability levels.

3.2 Proposed framework overview

The proposed framework consists of five sequential stages: data preprocessing and feature engineering, class imbalance handling, base model training, post hoc probability calibration, and Mondrian conformal prediction. Figure 1 illustrates the overall architecture.

3.3 Feature engineering

The feature engineering stage transforms raw RASFF notification records into a structured feature representation suitable for machine learning. The 21 features extracted from each notification are organized into three categories.

Categorical encodings are obtained by applying label encoding to the origin country, origin region, product category, hazard category, and notifying country. Temporal features include the notification year, month, and quarter, along with cyclical transformations using sine and cosine functions to capture periodic patterns Equations 3, 4): Binary seasonal indicators identify the summer (June–August) and winter (December–February) periods. Historical risk rates are computed for each categorical attribute using expanding window calculations to prevent data leakage:where c denotes a categorical value, t denotes the current time point, and n c , < t is the count of prior observations with category c.

month_sin = sin 2 π ⋅ month 12 (3)

month_cos = cos 2 π ⋅ month 12 (4)

rate c , t = 1 n c , < t ∑ i : c i = c , t i < t y i (5)

To prevent target leakage, each historical rate is computed as an expanding cumulative mean with self-exclusion over the chronologically sorted data: the rate assigned to observation i uses only strictly prior observations ( t j < t i ) that share the same category and never observation i’s own label. Because all rates are computed before splitting but are strictly backward-looking, a test observation’s rate depends only on events preceding it in time, regardless of split boundaries. This property was empirically verified: flipping the label of the latest observation leaves every earlier observation’s rate features numerically unchanged, confirming that no future information propagates backward. A sensitivity analysis quantifying the influence of these features is reported in Section 5.12.

3.4 Class imbalance handling

The RASFF dataset exhibits substantial class imbalance, with 79.1% of notifications classified as serious risk. To address this imbalance, the framework employs a dual strategy combining SMOTE with class weighting.

SMOTE generates synthetic minority-class samples by interpolating between existing minority instances and their k-nearest neighbors, as defined in Equation 6 (Chawla et al., 2002):where x i is a minority sample, x j is one of its k-nearest neighbors, and λ ∼ Uniform ( 0,1 ). A sampling strategy of 0.5 is applied, generating synthetic samples until the minority class reaches 50% of the majority class size.

x synthetic = x i + λ ⋅ x j − x i (6)

Class weights are computed inversely proportional to class frequencies (Equation 7):where n c is the number of samples in class c. These weights are incorporated into the XGBoost training objective.

w c = N 2 ⋅ n c (7)

3.5 XGBoost classification

XGBoost (eXtreme Gradient Boosting) is employed as the base classifier due to its strong empirical performance and ability to handle heterogeneous features (Chen and Guestrin, 2016). The model minimizes the regularized objective function given in Equation 8:where ℓ is the logistic loss, w y i is the class weight for sample i, f k represents the k-th tree, and Ω ( f k ) = γ T k + 1 2 λ ‖ w k ‖ 2 is the regularization term penalizing tree complexity.

L θ = ∑ i = 1 N w y i ⋅ ℓ y i , y ̂ i + ∑ k = 1 K Ω f k (8)

Threshold optimization is performed on the calibration set by maximizing the Cohen’s kappa coefficient (Equation 9):where p o is the observed agreement and p e is the expected agreement by chance.

κ = p o − p e 1 − p e (9)

3.6 Platt scaling calibration

Class balancing techniques distort the training distribution, leading to miscalibrated probability estimates. Platt scaling corrects this miscalibration by fitting a logistic regression model to the calibration set (Platt, 1999).

Let s ( x ) = log ( p ̂ ( x ) / ( 1 − p ̂ ( x ) ) ) denote the logit of the uncalibrated probability. The Platt calibrator learns parameters a, b to produce calibrated probabilities, as shown in Equation 10:where σ is the sigmoid function. The parameters are estimated by minimizing log loss on the calibration set.

p ̂ cal x = σ a ⋅ s x + b = 1 1 + exp − a ⋅ s x + b (10)

The expected calibration error (ECE) quantifies the quality of calibration (Equation 11):where B m is the set of samples in probability bin m, acc ( B m ) is the empirical accuracy, and conf ( B m ) is the average predicted probability.

ECE = ∑ m = 1 M | B m | N acc B m − conf B m (11)

3.7 Mondrian conformal prediction

Mondrian conformal prediction extends standard conformal prediction to provide group-conditional coverage guarantees (Vovk, 2003). The calibration procedure computes non-conformity scores and thresholds separately for each group.

For a calibration sample ( x i , y i ) with predicted probabilities p ̂ 0 , p ̂ 1 for classes 0 and 1, respectively, the non-conformity score is defined as Equation 12:representing the probability mass assigned to incorrect classes.

α i = 1 − p ̂ y i (12)

For each group g ∈ G, let I g = { i : g i = g } denote the indices of calibration samples in group g, and let n g = | I g |. The group-specific threshold is computed as shown in Equation 13:where the quantile level incorporates a finite-sample correction.

q ̂ g = Quantile α i : i ∈ I g , ⌈ n g + 1 1 − α ⌉ n g (13)

At prediction time, for a test sample x with group membership g, the prediction set is constructed as in Equation 14:including all classes whose non-conformity scores fall below the group-specific threshold.

C x = c ∈ 0,1 : 1 − p ̂ c x ≤ q ̂ g (14)

Algorithm 1 presents the complete Mondrian conformal prediction procedure.

Algorithm 1 Mondrian conformal prediction for food safety risk assessment.

3.8 Counterfactual explanation generation

To provide actionable insights for regulatory decision-making, the framework incorporates counterfactual explanation generation using three complementary methods.

Prototype interpolation identifies the minimum interpolation coefficient λ* (Equation 15) such that interpolating between the original instance and the target class prototype changes the prediction:where μ c is the mean of the training samples in the target class c, and λ * is found via binary search.

x cf = x + λ * ⋅ μ c − x (15)

Random perturbation generates counterfactuals by importance-weighted random modifications (Equation 16):where S is a randomly selected subset of features with selection probability proportional to feature importance, δ j ∼ N ( 0,0.3 σ j ), and e j is the j-th standard basis vector.

x cf = x + ∑ j ∈ S δ j ⋅ e j (16)

Growing spheres progressively expand the search radius until valid counterfactuals are found (Equation 17):where d ∼ N ( 0 , I d ) is a random direction and r is incrementally increased.

x cf = x + r ⋅ d ‖ d ‖ (17)

Risk factor analysis quantifies the effect size of each feature using Cohen’s d (Equation 18):where x ̄ j ( c ) and s j ( c ) are the mean and standard deviation of feature j for class c, respectively.

d j = x ̄ j 1 − x ̄ j 0 s j 1 2 + s j 0 2 / 2 (18)

4 Experimental setup

4.1 Dataset description

The experimental evaluation utilized the RASFF database, which contains food and feed safety notifications from European Union member states and associated countries. The dataset was filtered to include notifications from 2019 to 2025, resulting in 28,708 records. After removing entries with missing or ambiguous risk decisions, 22,643 notifications remain for machine learning analysis. To ensure transparency in the data preparation procedure, the following rules were applied: The binary label maps “serious” and “potentially serious” to the positive class and “not serious” and “no risk” to the negative class; the ambiguous categories “undecided” and “potential risk” were excluded. Records with unparseable notification dates were dropped. No exact duplicate rows were present; a small number of records sharing a reference number corresponded to follow-up notifications with distinct dates and were retained as separate events, consistent with RASFF reporting practice. Missing categorical fields were imputed to an explicit “unknown”/“other” level, historical rate features with no prior history defaulted to the neutral prior 0.5, and count features defaulted to zero. Temporal ordering was preserved so that all subsequent feature construction and splitting were strictly chronological.

Table 3 summarizes the dataset characteristics. The binary risk classification distinguishes between serious/potentially serious risks (the positive class) and not serious/no risk cases (the negative class). The substantial class imbalance, with 79.1% serious cases, motivates the dual SMOTE and class weighting approach.

TABLE 3. Dataset statistics for the Rapid Alert System for Food and Feed food safety notifications.

The data were temporally split to prevent data leakage and simulate realistic deployment conditions. The first 60% of the notifications form the training set; the next 20% constitute the calibration set for threshold optimization, Platt scaling, and conformal prediction; and the final 20% serve as the held-out test set. Because the calibration partition serves several purposes, a robustness check was performed using a fully disjoint design, in which threshold selection, probability calibration, and conformal quantile estimation each use non-overlapping subsets (60% training/10% validation/10% calibration/20% test). The disjoint design yields essentially identical performance and coverage to the shared-calibration design (with differences below 0.01 in Cohen’s kappa and below one percentage point in coverage), indicating that the multipurpose use of the calibration partition does not materially affect the reported conclusions.

The geographic distribution across origin regions includes the EU (7,862 notifications, 34.7%), Asia (6,241, 27.6%), other (5,631, 24.9%), the Americas (1,790, 7.9%), and Africa (1,119, 4.9%). Because the split is temporal rather than random, the regional composition is not identical across the training, calibration, and test partitions; a chi-square test of independence between region and split confirms a statistically significant association (χ 2 = 77.3, p < 0.001). This distributional shift is an expected consequence of temporal splitting and of evolving notification patterns over time, and it is reported explicitly rather than assumed away; the Mondrian construction is precisely intended to maintain group-wise coverage under such regional variation. Hazard categories span microbiological contamination, mycotoxins, pesticide residues, heavy metals, and allergens, with microbiological hazards representing the largest category. Product categories include meat, seafood, produce, nuts and seeds, spices, grains, and other products.

4.2 Baseline methods

The proposed XGBoost-based framework was compared against several baselines. The random classifier randomly assigns predictions according to class proportions, establishing a lower bound. The majority classifier always predicts the majority class (serious risk), providing an accuracy reference. A random forest with 200 trees and balanced class weights serves as an ensemble baseline. Gradient boosting with 150 estimators provides a comparison with traditional boosting. LightGBM offers an alternative gradient boosting implementation.

Table 4 summarizes the baseline configurations and their calibration set performance.

TABLE 4. Baseline method configurations and calibration set performance.

4.3 Evaluation metrics

Classification performance was evaluated using multiple metrics to account for class imbalance. Accuracy measures overall correctness, but it may be misleading under imbalance. Balanced accuracy averages recall across classes. Precision quantifies the proportion of predicted positives that are true positives. Recall (sensitivity) measures the proportion of actual positives correctly identified. The F1 score provides the harmonic mean of precision and recall.

Cohen’s kappa coefficient measures agreement beyond chance, with values above 0.2 indicating fair agreement. The Matthews Correlation Coefficient (MCC) provides a balanced measure suitable for imbalanced data, ranging from − 1 to + 1. The area under the ROC curve (AUC-ROC) quantifies discriminative ability across all thresholds.

Calibration quality is assessed via the expected calibration error (ECE) and the Brier score. The ECE measures the average gap between the predicted probabilities and the empirical frequencies. The Brier score is the mean squared error of probability predictions (Equation 19):

Brier = 1 N ∑ i = 1 N p ̂ i − y i 2 (19)

Conformal prediction performance was evaluated through coverage (the proportion of prediction sets containing the true label), average set size (an efficiency measure), and the singleton rate (the proportion of prediction sets containing exactly one class).

4.4 Implementation details

All experiments were implemented in Python using the scikit-learn, XGBoost, LightGBM, and imbalanced-learn libraries. Table 5 reports the hyperparameter settings for the proposed framework.

TABLE 5. Hyperparameter settings for the proposed framework.

Bootstrap confidence intervals with 1,000 resamples were computed for key metrics. The classification threshold was optimized on the calibration set by maximizing Cohen’s kappa over the range [ 0.1 , 0.9 ] with a step size of 0.01. All experiments use a fixed random seed of 42 for reproducibility.

5 Results and analysis

5.1 Data characteristics

Figure 2 presents the exploratory data analysis of RASFF notifications. The class distribution (panel A) confirms the substantial imbalance with 79.1% serious risk cases. The regional distribution (panel B) shows the EU and Asia as the dominant origin regions, collectively accounting for 62.3% of notifications. The hazard-type distribution (panel C) reveals microbiological contamination as the most frequent hazard category.

Risk rates vary substantially across regions (panel D), with the Americas exhibiting the highest serious risk rate and Africa showing lower rates. Hazard-specific risk rates (panel E) indicate that mycotoxins and allergens are associated with elevated serious risk probabilities compared with microbiological hazards. Product category analysis (panel F) shows diverse notification volumes across food categories.

Figure 3 illustrates temporal patterns in the data. Monthly notification volumes (panel A) show considerable variability with no obvious trend during the study period. Annual trends (panel B) reveal stable notification volumes from 2021 onward following an initial ramp-up period. Seasonal patterns (panel C) suggest slightly elevated risk rates during the winter months. The temporal data split (panel D) visualizes the training, calibration, and test set boundaries.

5.2 Geographic analysis

Figure 4 presents the geographic distribution of notifications and risk rates. Panel A shows the top 15 origin countries, with Türkiye, Poland, and India leading in notification volume. The color coding indicates regional membership, with EU countries (blue), Asian countries (yellow), and other regions represented. Panel B displays risk rates by country for those with at least 100 notifications, revealing that Iran and Türkiye have the highest serious risk rates.

5.3 Feature analysis

Figure 5 presents feature importance and correlation analysis. The feature importance plot (panel A) reveals that the hazard category rate dominates with an importance value exceeding 0.20, followed by the hazard category and the notifying country encoding. Historical rate features (highlighted in orange) collectively contribute substantial predictive power, validating the temporal feature engineering approach. The correlation matrix (panel B) shows moderate positive correlations among historical rate features (0.26–0.48) and a notable negative correlation between seasonal indicators (−0.32).

5.4 Model comparison

Figure 6 compares model performance across classifiers. Panel A shows the kappa and AUC scores for each model on the calibration set. XGBoost achieves the highest kappa (0.389) and shares the highest AUC (0.808) with LightGBM. The performance gap between gradient boosting methods and random forest is notable, with kappa improvements of 0.047–0.075. Panel B contrasts the proposed XGBoost framework with baseline classifiers, demonstrating substantial improvement over the random (kappa = 0) and majority (accuracy = 0.823) baselines.

5.5 Threshold optimization

Figure 7 illustrates the threshold optimization process. Panel A plots kappa and F1-scores against the classification threshold, revealing an optimal threshold of 0.28 that maximizes kappa. The optimal threshold is substantially below the default of 0.5, reflecting the class imbalance where a lower threshold is needed to achieve adequate recall for the minority class. Panel B shows the precision–recall trade-off, with precision increasing and recall decreasing as the threshold increases.

5.6 Classification performance

Figure 8 presents comprehensive classification results on the test set. The confusion matrix (panel A) shows 3,183 true positives, 367 true negatives, 436 false positives, and 543 false negatives at an optimal threshold of 0.28. The ROC curve (panel B) demonstrates discriminative ability with an AUC of 0.763. The precision–recall curve (panel C) shows an average precision of 0.936, substantially above the baseline (0.82).

Panel D summarizes the performance metrics, with an accuracy of 0.784, a balanced accuracy of 0.656, a precision of 0.880, a recall of 0.854, an F1-score of 0.867, a kappa of 0.296, and an MCC of 0.297. The score distribution (panel E) shows separation between classes but with substantial overlap, explaining the moderate kappa values. Bootstrap confidence intervals (panel F) indicate an accuracy of 0.784 [0.772, 0.796] and a kappa of 0.296 [0.265, 0.331].

Table 6 summarizes the final test set performance metrics.

TABLE 6. Final test set performance metrics with 95% confidence intervals.

The additional metrics make the difficulty of the minority (not serious) class explicit: specificity is 0.457, and the negative predictive value is 0.403, both substantially below the headline accuracy of 0.784, which is inflated by the 79.1% majority class. Because a false negative (a serious notification classified as not serious) carries greater public health consequences than a false alarm, the false-negative rate of 0.146 is reported explicitly. A cost-sensitivity scan over false-negative-to-false-positive cost ratios of 1:1 to 10:1 is provided in (Supplementary Material Table S1). These figures motivate the more cautious framing of the practical value adopted in Section 6.4: the framework is a decision-support prototype rather than a regulatory-ready system.

5.7 Probability calibration

Figure 9 demonstrates the effect of Platt scaling on probability calibration. The reliability diagram (panel A) compares calibration before and after Platt scaling. Before calibration, the model exhibits severe overconfidence, with the calibration curve lying substantially above the diagonal for low predicted probabilities. After Platt scaling, the calibration curve closely tracks the diagonal, indicating well-calibrated probabilities.

Panel B quantifies the calibration improvement, with the ECE decreasing from 0.267 to 0.043 (an 84% reduction) and the Brier score decreasing from 0.222 to 0.127 (a 43% reduction). Panel C shows the probability distribution shift, with precalibration probabilities heavily concentrated near 1.0 and postcalibration probabilities exhibiting better spread.

Because temperature scaling is a widely used single-parameter alternative, both methods were compared on the same calibration data. Temperature scaling (fitted temperature T = 0.855) did not improve upon Platt scaling, yielding a higher expected calibration error (0.056 vs. 0.043). This is consistent with the observation that Platt scaling has two parameters that better accommodate the sigmoidal miscalibration typical of gradient-boosted trees, whereas single-parameter temperature scaling was designed primarily for deep networks. Non-parametric alternatives (isotonic regression and histogram binning) achieved comparable calibration but carry a higher risk of overfitting on a moderately sized calibration set. Platt scaling was therefore retained as the calibration method.

Table 7 summarizes calibration metrics before and after Platt scaling.

TABLE 7. Calibration metrics before and after Platt scaling.

5.8 Conformal prediction results

Figure 10 presents the Mondrian conformal prediction results. The prediction set size distribution (panel A) shows that 78.4% of predictions are singletons (set size = 1), indicating high efficiency, while 21.6% include both classes (set size = 2). The average set size of 1.216 represents a favorable balance between coverage and efficiency.

Coverage calibration (panel B) demonstrates that empirical coverage closely tracks target coverage across α levels from 0.05 to 0.20. The confidence distribution (panel C) shows a mean confidence of 0.811. The coverage–efficiency trade-off (panel D) illustrates that smaller α (higher coverage targets) requires larger prediction sets.

Per-region coverage (panel E) shows Asia (90.9%), the EU (91.1%), and other (92.1%) meeting the 90% target, while the Americas (87.4%) and Africa (77.6%) fall below it in this single split. To distinguish genuine shortfalls from sampling noise, 95% Wilson confidence intervals were computed for each group (Table 8). The 90% target lies within the interval for the Americas ([0.836, 0.904]); therefore, its deviation is not statistically significant. In contrast, for Africa, the target lies outside the interval ([0.718, 0.825]), indicating a genuine shortfall in this particular split that is attributable to the smallest calibration and test sizes. As shown in Section 5.13, the mean coverage for every region across many random resplits converges to the nominal level, indicating that these single-split deviations reflect finite-sample variance rather than a systematic failure of the guarantee. Set sizes by region (panel F) show variation, ranging from 1.10 (the Americas) to 1.33 (other).

TABLE 8. Per-group empirical coverage with 95% Wilson confidence intervals (single-test split). The target ( 1 − α = 0.90 ) lies inside the interval for all groups except Africa.

Figure 11 provides a detailed Mondrian coverage analysis. Panel A compares the achieved versus the target coverage for each region, highlighting the valid (green) and invalid (red) coverage. Panel B shows the coverage gap, with Africa exhibiting a −12.4% gap below the target.

Table 9 summarizes Mondrian conformal prediction results by geographic region.

TABLE 9. Mondrian conformal prediction results by geographic region.

5.9 Error analysis

Figure 12 presents an error analysis across different dimensions. Error rates by region (panel A) show Africa with the highest error rate (approximately 27%), consistent with the coverage shortfall observed in the conformal prediction results. The Americas and Asia exhibit moderate error rates, while the EU and other regions show lower error rates. Error rates by hazard type (panel B) reveal that allergens and heavy metals have lower error rates than mycotoxins and pesticides.

The prediction breakdown (panel C) shows that true positives dominate (3,183), reflecting both the class imbalance and high recall. The error rate by confidence (panel D) demonstrates the expected inverse relationship: low-confidence predictions (0.0–0.2 bins) exhibit error rates exceeding 50%, while high-confidence predictions (0.8–0.9 bins) have error rates below 20%.

5.10 Risk factor analysis

Figure 13 presents a risk factor analysis based on Cohen’s d effect sizes. Panel A shows effect sizes of the top 10 features that distinguish high-risk from low-risk predictions. The hazard category rate exhibits the largest positive effect size (d = 1.055, large), indicating that notifications involving hazard types with historically high-risk rates are substantially more likely to be classified as serious. The hazard category count shows a large negative effect (d = −1.013), suggesting that hazard categories with more historical notifications tend toward lower current risk classifications.

Moderate effects are observed for the origin rate (d = 0.667), hazard category encoding (d = −0.504), and the product category rate (d = 0.502). Panel B shows the directional differences in the standardized feature means between high- and low-risk predictions.

Table 10 summarizes the top risk factors with their effect sizes and significance.

TABLE 10. Top risk factors with effect sizes and statistical significance.

5.11 Counterfactual explanations

Figure 14 presents the results of counterfactual explanation. Panel A shows the distribution of counterfactual generation methods, with prototype interpolation accounting for 64% of successful generations and random perturbation accounting for 36%. Panel B displays the distance distribution, with a mean distance of 1.45 indicating that relatively small perturbations are sufficient to flip predictions.

Panel C identifies the most frequently changed features in counterfactuals, with the hazard category count, the origin rate, and the notifying country encoding requiring the largest modifications. Panel D provides an example counterfactual, illustrating the specific feature value changes needed to convert a high-risk prediction to a low-risk prediction.

From the 25 generated counterfactuals, the analysis reveals that reducing the hazard category rate and increasing the hazard category count are the most actionable modifications for changing risk classifications. This aligns with the risk factor analysis showing these features as the strongest predictors.

5.12 Feature-ablation sensitivity

To address the possibility that the historical rate features dominate the model or introduce leakage-driven optimism, a sensitivity analysis retrained the full pipeline with these features removed. Excluding all four historical rate features changes Cohen’s kappa by less than 0.01 and does not reduce the AUC-ROC. The rates are individually correlated with the label (Pearson r between 0.12 and 0.35); however, they are largely redundant with the label-encoded categorical features, which the gradient-boosted trees can exploit directly. The practical implication of the leakage concern is important: because performance is essentially unchanged when the historical rates are removed, any residual leakage they might carry cannot be inflating the reported results.

5.13 Multitrial coverage stability

Because empirical coverage in a single calibration/test split is subject to sampling variability, the conformal evaluation was repeated over 5,000 random resplits of the pooled calibration and test data, recomputing the group-wise non-conformity thresholds in each trial. The mean coverage for every region converges to the nominal 90 % level (Africa 0.902, the Americas 0.901, Asia 0.900, the EU 0.900, and others 0.901; overall 0.901), with the smallest regions exhibiting the widest trial-to-trial dispersion (standard deviation up to 0.022 for Africa versus 0.010 for the EU). This confirms that the below-target coverage observed for Africa in the single reported split reflects finite-sample variance rather than a systematic breakdown of the group-conditional guarantee and that the coverage is centered on the target rather than skewed.

5.14 Summary dashboard

Figure 15 presents an integrated summary dashboard consolidating key results. The dashboard displays a confusion matrix, a ROC curve (AUC = 0.763), a regional coverage comparison, performance metrics (accuracy, kappa, MCC, and F1-score), a calibration curve (ECE = 0.0427), and a prediction set size distribution.

6 Discussion

6.1 Key findings

The experimental results demonstrate that the proposed framework successfully addresses the identified research gaps while achieving competitive classification performance. The XGBoost-based classifier with SMOTE and class weighting achieved a kappa of 0.296 and an MCC of 0.297 on the held-out test set, indicating fair agreement beyond chance. While these values are moderate, they reflect the inherent difficulty of food safety risk prediction where many factors influencing risk severity may not be captured in structured notification data.

The Platt scaling calibration component reduced the ECE by 84%, addressing the probability miscalibration induced by class balancing techniques. This improvement is critical for regulatory applications where calibrated risk scores enable the principled prioritization of inspection resources. A precalibration ECE of 0.267 indicated that the uncalibrated model would substantially overestimate risk probabilities, potentially leading to inefficient resource allocation.

The Mondrian conformal prediction component achieved 90.3% overall coverage with an average set size of 1.216, demonstrating that uncertainty quantification can be provided without sacrificing efficiency. The 78.4% singleton rate indicated that the framework provides definitive predictions in the majority of cases while appropriately expressing uncertainty through larger prediction sets when warranted. Four of the five geographic regions achieve valid coverage, although the African region falls short due to limited calibration data (n = 232).

6.2 Mapping results to research gaps

The experimental results provide evidence that the proposed framework addresses each identified research gap. Regarding gap G1 (lack of uncertainty quantification), the Mondrian conformal prediction component provides prediction sets with 90.3% empirical coverage, closely matching the 90% target. This is, to our knowledge, among the first applications of conformal prediction to RASFF-based food safety risk assessment, enabling decision-makers to reason with quantified confidence.

Gap G2 (probability miscalibration) is addressed through Platt scaling, which reduces the ECE from 0.267 to 0.043. The reliability diagram (Figure 9) demonstrates that the postcalibration probabilities closely track the diagonal, indicating well-calibrated predictions. The Brier score reduction from 0.222 to 0.127 confirms improved probabilistic accuracy.

Gap G3 (absence of group-conditional coverage) is addressed through the Mondrian extension, which maintains per-region thresholds. Coverage meets the target for the EU (91.1%), Asia (90.9%), and others (92.1%). For the Americas (87.4%), the target lies within the 95% confidence interval; therefore, the deviation is not statistically significant. Africa (77.6%) exhibits a genuine shortfall in the single reported split, driven by the smallest calibration and test sizes. The multitrial analysis (Section 5.13) shows that the mean coverage for all regions, including Africa, converges to the nominal level; thus, the Mondrian construction provides the intended guarantee on average, with single-split variability that is largest for small groups.

Gap G4 (limited interpretability) is addressed through counterfactual explanation generation and risk factor analysis. The counterfactual module identifies the minimal modifications needed to change predictions, providing actionable insights for regulatory intervention. Risk factor analysis reveals that hazard category (d = 1.055) and origin (d = 0.667) rates are the strongest predictors, guiding targeted monitoring efforts.

Gap G5 (no integrated framework) is addressed through the synergistic combination of SMOTE, Platt scaling, and Mondrian conformal prediction. The temporal data split and calibration set design ensure that each component receives the appropriate data for training, calibration, and evaluation.

6.3 Limitations

Several limitations should be acknowledged. Moderate kappa and MCC values (approximately 0.30) indicate that substantial classification uncertainty remains, reflecting factors not captured in structured notification data such as outbreak context, contamination levels, and product volumes. The African region does not achieve valid conformal coverage due to limited calibration data, suggesting that minimum group size requirements should be carefully considered in deployment.

The framework assumes temporal stationarity in risk patterns, which may not hold during food safety crises or regulatory changes. The RASFF database reflects EU reporting practices and may not generalize to other regulatory systems. Additionally, the counterfactual explanations assume feature independence and may suggest modifications that are not practically feasible.

6.4 Practical implications

For regulatory authorities, the framework offers several practical benefits. Calibrated probability estimates enable risk-based prioritization of inspection resources, directing attention to notifications with the highest predicted risk. The conformal prediction sets communicate uncertainty to decision-makers, supporting more nuanced risk communication than binary classifications alone.

The regional coverage analysis identifies geographic areas where prediction reliability may be compromised, informing data collection priorities. The risk factor analysis guides targeted monitoring by identifying hazard categories, origin countries, and product categories associated with elevated risks.

7 Conclusion

This study presented a comprehensive machine learning framework for food safety risk assessment that integrates XGBoost classification, SMOTE-based resampling, Platt scaling calibration, and Mondrian conformal prediction. An experimental evaluation of 22,643 RASFF notifications demonstrated that the framework achieves competitive classification performance (accuracy = 0.784, kappa = 0.296, AUC = 0.763) while providing calibrated probability estimates (ECE = 0.043) and prediction sets with 90.3% coverage and an average set size of 1.216.

The framework addresses critical gaps in existing food safety prediction systems by providing distribution-free uncertainty quantification, correcting calibration distortion from class balancing, and empirically characterizing group-wise coverage across geographic regions. Counterfactual explanation generation and risk factor analysis provide interpretable insights. Given the modest agreement beyond chance (Cohen’s kappa near 0.30) and the reduced minority-class performance, the framework is positioned as a decision-support prototype that requires further validation before regulatory deployment rather than a regulatory-ready system.

Future work may explore several directions. Adaptive conformal prediction methods could adjust to distribution shifts during food safety crises. Multitask learning could jointly predict risk severity and hazard types. Natural language processing could extract additional features from notification text descriptions. Finally, integration with real-time surveillance systems could enable the prospective validation of prediction performance.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frfst.2026.1890831/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Fernando Perez Rodriguez, University of Cordoba, Spain

Reviewed by: Mohamed Farag, University of Bonn, Germany

Jia-Qiang Lv, Tsinghua University, China