Quality Evaluation of Honey Using Spectral Analysis and Chemometric Regression Methods

1. Introduction

The quality evaluation of food and agricultural products is a main topic with significant meaning and growing interest in the scientific field. The increasing efforts to ensure authenticity in the food industry have led to collective actions from producers and suppliers, together with regulators and researchers, to address the entire production and distribution cycle.

Honey and bee products are a true natural gift to humans. Honey is rightly called a superfood with a wide range of uses – not only delicious, high-nutritional, and healthy for everyday consumption, but also healing for internal and external use. The honey primarily consists of various carbohydrates, and its composition can vary depending on factors such as environmental conditions, climate, plant species, and also beekeeping practices.

Bulgaria has a large beekeeping production and a sustainable part in the honey distribution chain in Europe. Controlling honey quality indicators from the raw material to the product on the store shelf is a key factor in maintaining a good reputation for Bulgarian honey. The task is certainly difficult because of the diverse nature of that biological product. Honey purity, bio origin, and absence of harmful chemical ingredients present a growing challenge to food safety, food fraud detection, and consumer confidence.

The purpose of the qualification is to determine the quality and plant origin, ecological purity, and possible adulterations, as well as to establish whether it was obtained, processed, and stored under optimal conditions. This is achieved through a comprehensive study of the organoleptic and physicochemical properties, the amount of its main components, and various contaminants, combined with microscopic examination of the sediment after centrifugation of the aqueous solutions of honey /pollen analysis/. Pollen analysis is very often used to determine the plant origin of honey. It is based on the fact that honey always contains pollen from the plants from which the nectar was collected, on the one hand, and the differences in the size, shape, and structure of the shells of the pollen grains of individual plants.

Regarding the methods for analyzing honey, there is a trend towards the unification of analytical methods worldwide and the introduction of modern, accurate, and rapid methods with high reproducibility. The International Honey Commission has proposed such methods, which are recognized as official and are used in quality control. These methods are also included in the Bulgarian Regulation on Sampling and Methods of Analysis in the Control of Honey (Bogdanov et al., 1997; Ivanov, 2006).

Traditional methods for determining the physicochemical properties of honey are quite laborious, require special skills, consumables, and good analytical skills. In addition, they take a lot of time. This can be saved with the help of rapid analytical techniques, which can quickly and easily detect non-compliance with regulatory requirements. In recent years, many studies have been conducted to find fast and inexpensive methods for the qualification of food products.

As an alternative to traditional methods, spectroscopic analysis in the near-infrared region (near–infrared spectroscopy – NIRS) can be employed. A significant difference of this method compared to the others is that it analyzes the samples without preliminary chemical treatment, dissolution, or separation of components. The advantage of the method is that it is fast and does not require anything other than a spectrophotometer, software for processing spectral data, and a sample. In addition, the integrity of the sample is preserved, and it can be used for other analyses. Usually, in quantitative analysis, some kind of linear model is sought for the relationship between the NIRS information and the chemical information. The obtained results show good accuracy in determining glucose, fructose, di- and trisaccharides, moisture, electrical conductivity, optical activity, etc., which are important criteria for evaluating honey.

The main areas in which authors work on the topic of honey analysis are as follows: floral and geographic origin detection [Bogdanov, Maionea], authenticity [Zhang, 2023], adulteration recognition [Furong Huang, 2020; Md Abdullah Al Noman,2025; Hategan], quality assessment [Truong,2022; Zhang, 2021], and determination of chemical content [Silva; Koraqi, Nunes]. Information concerning the chemical content and presence of heavy metals in honey is crucial for producing high-quality honey for the food chain market.

[Bogdanov] reviewed multitude of known physico-chemical methods used for the determination of the botanical origin of honey. Using statistical methods such as principal component analysis (PCA); linear discriminant analysis (LDA) and cluster analysis (CA) the authors concluded, that a good classification of unifloral honeys in regard to different honey quality parameters - sugars; electrical conductivity; optical rotation; acidity, proline; and nitrogen content could be achived.

Routine physicochemical methods are used in the classical approach for the authentication of botanical origin, together with pollen and melissopalynological analysis. These methods are mostly validated and widely used in the routine control of honey, but this approach is laborious and depends on highly specialized personnel.

In the available scientific literature on the topic for assessing the origin of honey employment of several different analytical techniques such as gas chromatography coupled to mass spectrometry, vibrational spectroscopy, NMR, and physico-chemical properties has been discussed.

[Silva] This study used the inductively coupled plasma optical emission spectrometry technique (ICP OES) to determine the concentration of metals (Al; Ca; Cu; Fe; Li; K; Mg; Mn; Na and Zn) in honey samples to determine geographic origin and authenticity by use of statistical techniques and Machine Learning. This approach proved very efficient in discriminating the origin of the honey samples.

The study of [Kassa] evaluated heavy metal concentrations in honey samples collected from 5 regions in Ethiopia to assess the safety of the honey for consumption. Heavy metal content was quantified using an Atomic Absorption Spectrophotometer after sample digestion with an acid solution. Physicochemical parameters including moisture content; electrical conductivity; pH; ash content; and temperature were analyzed. The results showed that all honey samples met international quality standards.

[Koraqi H. ] This study examined the standard physicochemical properties; bioactive components and antioxidant activity of Kosovan honeys according to their floral source. The Kohonen Neural Network (KNN); and Principal Component Analysis (PCA) were used to identify and group botanical origin and the results showed the usefulness of two methods for dimensionality reduction and detecting the structure and general regularities.

[Maionea] in their review the application of multivariate data analysis and machine learning techniques for exploratory and prediction analysis of honey origin are reported. Most of the studies reviewed base their methods on atomic spectra and physicochemical properties as descriptive variables for honey. Sensorial data obtained from electronic tongue and nose; and color histograms of honey images also show high discriminative power for ascertaining honey origin. Principal component analysis (PCA); discriminant analysis (DA), and cluster analysis are the preferred techniques for performing exploratory and predictive analyses for the purpose of identifying origin.

Spectroscopic sensing techniques combined with chemometric approaches have been widely used for the quality evaluation of food and agro-products. The deep learning approach combined with NIR spectroscopic sensing techniques has shown great potential for quality evaluation of food and agro-products. Current advances in deep learning-based qualitative analysis include variety identification, geographical origin detection, and bruise detection, whereas quantitative analysis includes multiple-component content prediction for fruits, grains, and crops.

To obtain the information from spectral data, many linear and nonlinear chemometric methods, including partial least squares (PLS), artificial neural network (ANN), and support vector machine (SVM), have been developed to quantify the properties of food and agro-products. However, these nonlinear models usually produce many parameters and are confronted with the risk of overfitting.

Artificial neural networks (Schuhfried et al., 2016), decision and regression trees (Popek et al., 2017), regression models (Wei and Wang, 2014), and support vector machines (Gan et al., 2016) are among the popular machine learning algorithms for data classification, which are used in honey analyses. The CNNs are the first developed and widely used deep learning approaches for spectral analysis. In recent years, researchers have focused on deep learning methods. The critical characteristic of deep learning is the data-driven strategy for extracting hidden features from raw data without human engineering and prior knowledge.

[Nunes] aimed to determine the chemical profile of floral honey by ultraviolet and visible spectrophotometry (UV–vis) and near-infrared (NIR) vibrational spectroscopy. The analysis includes also the total phenolic content (TPC), considering two harvest years, geographic and botanical origins, and seasonality of samples collected in southern Brazil. The authors concluded that the total phenolic content (TPC) of honey showed great variation by geographic origin. Their work showed that the UV–vis range of 260–360 nm coupled to PCA can distinguish honey with higher and lower TPC. The NIR spectral data set revealed discrepancies in honey according to the harvest seasons. Nunes et. Al. proved UV–vis and NIR spectroscopies coupled to PCA to be a useful approach to distinguish honey samples as to possible adulterations and geographic origin associated with their higher and lower phenolic contents.

The study of Truong et al. [Truong] used analysis of spectral information (based on important spectral bands) and showed the good ability of Vis-NIRS to capture unique chemical information related to potency and purity. In the same study, Methylglyoxal (MGO) is set as the primary antibacterial in manuka honeys, but there are other chemical markers contributing to its value. Chemometric models using Partial Least Squares (PLS), Support Vector Machine (SVM), and discriminant analysis, PLS-DA, and SVM-DA (RBF) were developed to predict manuka honey potency and purity. The authors achieved a prediction of potency with 74 % accuracy by a PLS regression model. For prediction of botanical origin, PLS-DA and SVM-DA classifiers - 89 % overall accuracy, but for mono-floral manuka honeys, especially 92–97 % accuracy.

Furong Huang used Near-infrared spectroscopy (NIR) and mid-infrared (MIR) spectroscopy data for the detection of adulteration in honey together with an SVM classifier [Furong Huang; 2020]. The pre-processing step of raw data includes First Derivative (FD) transform, Second Derivative (SD) transform, Multiple Scattering Correction (MSC), and Standard Normal Variate Transformation (SNVT), as well as low-level and intermediate-level data fusion. Then, Grid Search (GS); Genetic Algorithm (GA), and Particle Swarm Optimization (PSO) were employed as the optimization algorithms to find the optimal SVM parameters and to establish the best model. The results reveal that intermediate-level data fusion significantly improves the detection model, as the accuracy, sensitivity, and specificity of the optimal SVM model reach 100%.

[Md Abdullah Al Noman;2025] proposed a rapid approach for authenticating honey and detecting multiple adulterants using visible, ultraviolet, and near infrared spectroscopy (200–900 nm). PCA was used to extract informative spectral regions. The Random Forest (RF); Support Vector Machine (SVM), and Soft Independent Modeling of Class Analogy (SIMCA) were used for predicting. RF model achieved the highest performance with classification accuracies between 99-100 % depending on adulterant type. Principal Component Analysis (PCA) was used to extract informative spectral regions, and the following ranges were established as significant: 350–450 nm for botanical classification; 200–600 nm for adulteration detection; and 596–605 nm for differentiating adulterants.

[Hategan] developed models for the prediction of honey botanical origin and harvesting year based on isotope and elemental fingerprints. For each classification criterion, significant focus was dedicated to the data preprocessing phase to enhance the models. The application of Artificial Neural Networks (ANNs) is highlighted as efficient, as the obtained classification performances, including accuracy scores, are higher than 86% during the test set.

In this background, the purpose of the present study is to find a rapid, accurate, and effective single analytical method for analyzing honey samples for heavy metal contamination. NIR spectroscopic sensing techniques combined with regression analysis were used for predicting the presence of heavy metals - arsenic (As), cadmium (Cd), lead (Pb), iron (Fe), Hg, and pH, amount of pH, reducing sugars, sweet disaccharide, and water content.

2. Materials and Methods

2.1. Honey Samples Collection and Spectral Data Acquisition

Honey samples from 14 regions of Bulgaria formed a set of 29 samples included in the study. The following regions of production took part in the study: Stara Zagora, Borovo, Pamukchii, Pustrene (Stara Zagora region), Ruse, Brestovica, Yuper, Modjereto, Novo selo, Ivanovo (Ruse region), Razgrad, Haskovo, Turgovishte, Sliven. Thirteen types of honey floral origin were selected: Rapeseed; Sunflower; Linden; Multiflower; Sunflower; Rapeseed & Amorpha; Acacia & Mana; Draka & Pustren; Acacia; Lavender; Mana; Multiflower; Acacia & Rapeseed.

Samples are preliminarily evaluated in the certificated laboratory by chemical analysis for the determination the content of the following metals: arsenic (As); cadmium (Cd); lead (Pb); iron (Fe), and components: amount of pH; reducing sugars; sweet disaccharide, and water content shown in Table 1.

The principal methodology for chemical analysis of honey probes includes burning samples in an analyzer furnace, separation of mercury vapor, transformation into the form of amalgam, heating and release of free mercury, and determination of components by measuring the absorption of radiation with a wavelength of 253.64 nm. AMA 254 Mercury Analyzer (Altec – Czech Republic) was used for the analysis, an oxygen cylinder, and nickel boats for burning samples. As a software tool, AMA/ama.exe analyzer program was used. Our preliminary analyses (Penchev et.al., 2024) include a detailed presentation of laboratory results for all listed chemical parameters.

Preparation of samples for spectral data acquisition includes liquefying (decrystallization) in advance. To keep honey characteristics unaffected, melting at a temperature below 45ºC was carried out. Еach sample consisted of 10 ml decrystallized organic honey. A petri dish with a 60 mm diameter was used to collect the probes.

Figure 1. Honey samples 7, 17, 29.

Figure 1. Honey samples 7, 17, 29.

The spectral characteristics of the copper samples were obtained in daylight, placed in cuvettes, and using a special attachment for placing the probe on the spectrophotometer, as shown in Figure 2.

The spectral data of 29 honey samples are shown in Figure 3.

The obtained spectral data of the copper samples served as the basis for subsequent extraction of informative characteristics of the quality indicators. Data processing was performed using the MATLAB software.

2.2. Chemometric Regression Methods Used for Quality Evaluation of Honey Samples

In this study, two chemometric regression techniques, Partial Least Squares Regression (PLSR) and Principal Components Regression (PCR), were used in Matlab to determine the quality parameters of 29 honey samples. The selection of an appropriate regression analysis method was based on a review of studies conducted by other authors working on similar tasks [17,18,19].

PLSR and PCR are both methods to model a response variable when there are a large number of predictor variables, and those predictors are highly correlated or even collinear. Both methods construct new predictor variables, known as components, as linear combinations of the original predictor variables, but they construct those components in different ways.

PCR is a regression approach that first uses principal component analysis on the training data to reduce its dimensionality; then a regression (e.g., linear) is applied and trained on the resulting samples. In PCA, the transformation is completely unsupervised, meaning that no information about the targets is used. As a result, PCR may perform poorly on some data sets where the target is highly correlated with data with low variance. In effect, PCA dimensionality reduction projects the data into a lower-dimensional space where the variance of the projected data is highly maximized along each axis. Although they have the greatest predictive power on the target, the directions with lower variance will be discarded and the final regression will not be able to use them.

PLSP simultaneously performs data transformation and modeling. This method also applies dimensionality reduction to the samples before applying a linear regressor to the transformed data. The main difference with PCR is that in PLS, the transformation is supervised.

2.3. Evaluation Metrics

Currently, a wide range of metrics are used in classification tasks to evaluate the performance of regression models, which allow for numerical evaluation of the performance of the models. In this research coefficient of determination is used for models performance evaluation.

The coefficient of determination (R²) measures how well a statistical model predicts an outcome. The outcome is represented by the model’s dependent variable.

The lowest possible value of R² is 0 and the highest possible value is 1. The better model has closer R² to 1.

If the R² is 0, the linear regression model doesn’t allow to predict exam scores any better than simply estimating that everyone has an average exam score. If the R² is between 0 and 1, the model allows you to partially predict exam scores. If the R² is 1, the model allows you to predict precise exam score.

The coefficient R² is a measure of goodness of fit. It is the proportion of variance in the dependent variable that is explained by the model. It is calculate using the regression outputs, presented in Equation (1).

$$R^2=1-{\displaystyle \frac{\mathrm{R}\mathrm{S}\mathrm{S}}{\mathrm{T}\mathrm{S}\mathrm{S}}},$$

(1)

Where RSS is sum of squared residuals;

TSS – total sum of squares.

3. Results and Discussion

The results of the experimental study are divided into two parts: one concerns the assessment of the presence of heavy metals, and the other - of the physicochemical indicators.

3.1. Determination the Content of the Heavy Metals

To form a simpler model, a smaller number of components could be used. The percentage of variance explained in the response variable as a function of the number of components, presented in Figure 4 for the presence of heavy metals, is used to select the number of components.

It details the number of underlying dimensions on which most of the variance is observed. The first component makes up from 0 to 38% of the explained variance. This is relatively low as compared to other datasets, but no matter of concern. It simply indicates that a major share (100%–38%=62%) of observations distributes across more than one dimension.

In the study, we need to reduce the dimensionality of the data and retain at least 90% of the variance of the original data. Then, we need to include at least 3 components to reach at least 90% for arsenic, lead and iron; for cadmium - at least 5 components. From 7 components upwards, 100% variance is reached for all four types of metals.

A PCR model with two principal components was created. A principal component analysis of the spectral data was performed. In this case, PCR is a linear regression of the response variable on these two components.

The PCR results are plotted against the original spectral data, with the regression coefficients for the original, uncentered variables transformed.

Figure 5 shows the fit versus observed response for the PLSR and PCR fits.

The two PLS components are much better predictors of the observed y, but they explain slightly less variance in the observed X than the first two principal components used in the PCR.

Estimated mean square prediction error as a function of the number of components for arsenic (As) a), cadmium (Cd) – b), lead (Pb) – c), iron (Fe) – d) is shown in Figure 6.

The number of components can be chosen to minimize the expected error in predicting the response from future observations on the values of the quality parameters - the presence of heavy metals. When using a larger number of components, higher prediction accuracy is achieved, but this is often associated with greater computing power. Cross-validation is often used as a more statistically sound method for choosing the number of components in PLSR or PCR. It avoids overfitting of data by not reusing the same data for both model fitting and prediction error estimation. Thus, the prediction error estimate is not optimistically biased downward. In PLS and PCR, there is an option to estimate the mean squared prediction error (MSEP) through cross-validation.

The mean square error curve of the prediction of the As values for PLSR shows that two components do the best possible job, while PCR needs eight components to obtain similar prediction accuracy. When predicting the Cd values for PLSR and PCR, one component is sufficient to obtain good prediction accuracy. When predicting the lead values for PLSR shows that four components do the best possible job, while PCR needs six components to obtain good prediction accuracy. When predicting the iron values for PLSR and PCR, two components are sufficient to obtain good prediction accuracy.

The effectiveness of the two chemometric methods PLSR and PCR when two components are used for quantitative assessment of the values of heavy metals present in honey samples was assessed by the coefficient of determination indicator, the values of which are presented in Table 2.

The estimated mean squared error of prediction shows that PLSR has a better predictive ability for the presence of heavy metals in honey. The $R^2$values from both regressions confirm this.

3.2. Determination the Content of the Physicochemical Indicators

For the remaining physicochemical indicators - pH value; reducing sugars; sweet disaccharide and water content, a similar analysis was performed as for the presence of heavy metals.

Estimated mean square prediction error as a function of the number of components for amount of pH - a), reducing sugars – b), sweet disaccharide – c), water content – d) is shown in Figure 7.

The effectiveness of the two chemometric methods PLSR and PCR when two components are used for quantitative assessment of the values of physicochemical indicators present in honey samples was assessed by the coefficient of determination indicator, the values of which are presented in Table 3.

The estimated root mean square error of prediction shows that PLSR has a similar predictive ability to PCR for the presence of physicochemical indicators in honey. The $R^2$values from both regressions confirm this.

4. Conclusions

In this article, data from the near infrared spectrum are used to create regression models for assessing and predicting the main quality parameters of copper - the presence of heavy metals and basic physicochemical indicators. The proposed approach for determining the quality parameters of copper using spectral data and regression methods is suitable for the analysis of copper samples and could be easily implemented and used in industrial conditions for assessing the quality of copper. The achieved recognition accuracy meets the requirements in this area.

A comparative analysis of two main chemometric methods for refractometry analysis was made, and experiments were conducted with 29 samples of Bulgarian copper from different regions in the country.

The best results for predicting the main quality parameters of honey are shown by PLSR with values for the presence of heavy metals on $R^2$ from 0.4881 for cadmium to 0.7981 for lead; as well as for the evaluation of the main physicochemical parameters from 0.4061 for reducing sugar to 0.5693 for Sweet disaccharide.

The regression equations obtained from the relationship between the main quality parameters of honey and its spectral characteristics, together with the regression method, serve as the basis for the development of a compact tool for the evaluation of honey quality indicators.