

\section{Introduction}
Following the interesting results obtained during similar events \citep[see e.g.][and references therein]{pierna2011case,pierna2020applicability}, a chemometric challenge has been held during the inaugural edition of the ``International Workshop on Spectroscopy and Chemometrics'', organized by the Vistamilk SFI Research Centre in April 2021. 

The dataset provided to the participants contained the values for three milk quality traits for different milk samples along with the corresponding mid-infrared spectra. Mid-infrared spectroscopy (MIRS) represents a convenient and non-disruptive way to collect vast amount of data in a relatively cheap and fast way. In recent years, such data have been proven useful to predict several different quantities of interest in the dairy framework; see for example the results in \citet{cecchinato2009mid,mcparland2014mid,bonfatti2017comparison} for coagulation properties of milk, animal energy efficiency, and fatty acids concentration, respectively. Despite these advantages, MIRS data introduce some relevant challenges from a modelling perspective, such as high-dimensionality and the strong correlations among the variables (i.e. the wavenumbers). For this reason several different strategies have been explored in literature and readers may refer to the recent work by \citet{frizzarin_predicting_2021} to see how different machine learning and statistical techniques performed on such data. 

The primary objective of the competition was, therefore, to evaluate the ability of the participants to propose different regression strategies to predict milk quality traits. In particular, it has been required to provide accurate predictions for the traits exploiting only the information contained in the corresponding milk spectra. The variables to be predicted have been chosen in order to introduce some trait specific challenges, as it will be clear in the next section. 

Different participants, or groups of, took part in the competition. In this work, the best six contributions, evaluated in terms of their predictive performances, are presented focusing both on the results obtained and on the subjective choices the participants made in terms of pre-processing and modelling paradigms. The rest of the paper is structured as follows. The data that have been provided and the criteria considered to evaluate the contributions are introduced in Section \ref{sec:sec2}. In Section \ref{sec:sec3} the different modelling strategies devised by the participants are outlined. In Section \ref{sec:sec4}, a critical discussion of the results and of the proposed approaches, together with an overview on the recurring themes discussed during authors' presentations at the workshop, is provided. Lastly, Section \ref{sec:conclusion} concludes the paper with some considerations on the competition. 

\section{Data description and challenge}\label{sec:sec2}
A dataset containing 622 milk samples from 622 cows were collected between August 2013 and August 2014 from 7 different Irish research herds. The samples originated from Holstein-Friesian, Jersey and Norwegian Red cows, as well as their crosses; all cows were fed a predominantly grass-based diet with occasional concentrate and grass silage supplementation. The samples were collected during morning and evening milking at different stages of lactation and different parities. All samples were analyzed by the same MilkoScan FT6000 (Foss Electronic A/S, Hillerød, Denmark) producing 1060 transmittance data points in the mid-infrared light region. 

\begin{figure}[t]
    \centering
    \includegraphics[scale = 0.5]{figs/histograms.pdf}
    \caption{Frequency distributions of the three reference traits for the training set.}
    \label{fig:trait_hist}
\end{figure}

Milk pH, Casein Micelle Size (CMS) and $\kappa$-casein were provided for the analyses. Milk pH of all samples was assessed with a SevenCompact pH meter S220 (Mettler Toledo AG, Greifensee, Switzerland). The casein micelle hydrodynamic diameter was determined using a Zetasizer Nano system (Malvern Instruments Inc., Worcester, UK). Lastly, $\kappa$-casein was determined using reverse-phase high performance liquid chromatography (HPLC) using an adaptation of the method of Visser et al. (1991) and is expressed as grams per liter of milk. 

The original dataset was divided into training and test set. The training set contained 399 observations for $\kappa$-casein, 538 observation for CMS, and 548 observations for pH. Some spectra were therefore associated with only one or two of the reference traits, with the missing trait being considered as a missing value. Some summary statistics for the mentioned traits are reported in Table \ref{table:descriptive_stat} while a graphical representation of their frequency distributions is provided in Figure \ref{fig:trait_hist}. For all the observations the corresponding mid-infrared spectra, with measures for 1060 wavenumbers, were available. A test dataset has been provided to the participants, with 69 spectra with unknown values for the three reference traits. Therefore, the objective of the competition was to predict the test set values for $\kappa$-casein, CMS, and pH values by exploiting the spectral information only. Note that the information about the specific wavenumbers has not been shared with the participants; nonetheless they were aware that the region spanned by the provided spectra goes from 900 to 5000 cm$^{-1}$.
The three traits provided have been carefully selected according to their different characteristics \citep[see][]{frizzarin_predicting_2021}; CMS has been shown to have low predictability while $\kappa$-casein and pH has a medium-low and a medium-high predictability, respectively. Moreover, all these traits are important for milk processing.

The contributions from the participants have been ranked considering the \emph{Root Mean Squared Error} computed on the test dataset (RMSEP)
$$
\text{RMSEP}_{kj} = \sqrt{\frac{\sum_{i=1}^n (y_{ik} - \hat{y}^{(j)}_{ik})^2}{n}}
$$
where $y_{ik}$ and $\hat{y}^{(j)}_{ik}$ represents respectively the $i$-th reference value and the corresponding prediction for the $k$-th trait by the $j$-th participant, while $n$ is the total number of observations in the test set. In order to provide a single measure to rank the performances across all the traits a \emph{relative error} has been considered, that has been defined as 
\begin{eqnarray}\label{eq:rel_err}
\text{RERR}_j &=& \frac{1}{3}\left(\sum_{k=1}^3 \frac{\text{RMSEP}_{kj}}{\text{min}_j \,\, \text{RMSEP}_{kj}} \right)
\end{eqnarray}
with $j$ spanning over the total number of participants. Thus for each trait, each participant’s RMSEP is divided by 
the best RMSEP for that trait and this is summed up over the three traits.  The more the value of (\ref{eq:rel_err}) is close to 1, the better the corresponding method has performed over the available traits. The rationale behind the choice to consider (\ref{eq:rel_err}) to evaluate the results, lies in the necessity to have a relative measure being trait independent and not being influenced by the different error magnitudes and degree of difficulty for the three traits.

\begin{table}[t]
\centering
\caption{Mean, median, standard deviation (SD), minimum and maximum values for the provided traits.}
\begin{tabular}{l|ccccc}
  \hline

  Trait & Mean & Median & SD & Minimum & Maximum \\ 
  \hline
  $\kappa$-casein & 5.92 & 5.71 & 1.83 & 1.36 & 22.23 \\
  CMS & 229.99 & 172.45 & 343.10 & 63.12 & 4063.00 \\
  pH & 6.68 & 6.68 & 0.12 & 5.86 & 7.06 \\
   \hline
\end{tabular}
\label{table:descriptive_stat}
\end{table}

\section{Modelling approaches and results}\label{sec:sec3}

\subsection{Participant 1}\label{sec:sec3_1}
\input{mlgig-tabular}

\subsection{Participant 2}\label{sec:sec3_2}
\input{mlgig-timeseries}

\subsection{Participant 3}\label{sec:sec3_3}
\input{Giovanna}

\subsection{Participant 4}\label{sec:sec3_4}
\input{Ferr_Stef}

\subsection{Participant 5}\label{sec:sec3_5}
\input{WorkshopPaperKD}




\subsection{Participant 6}\label{sec:sec3_7}
First of all, an exploratory data analysis has been performed with Unscrambler X (Camo Analytics) and with the \texttt{R} software. The data were visualized by plotting the raw spectra, as well as the calculated interquartile-range of each wavelength which allowed to identify two highly noisy regions. Several spectral pre-treatments were tested such as Standard Normal Variate, first and second derivative Savitzky-Golay transformations and Multiplicative Scatter Correction (MSC) with the latter one being chosen as the most promising. According to the identification of the high-noise-level regions, in the subsequent analyses two approaches have been taken: in the first case the original 1060 wavenumbers have been considered while, in the latter, the noisy regions were removed leading to datasets with 533 features. In both the cases, CMS values were log-transformed to obtain data following more closely a Gaussian distribution.

\begin{figure}[t]
    \centering
    \includegraphics[width = 16.3cm, height = 7.5cm]{figs/image1AK.pdf} \\
    \vspace*{-1cm}
    \hspace*{-7.85cm} \includegraphics[width = 7.8cm, height = 7.5cm]{figs/image2AK.pdf} 
    \caption{Variables selected (red vertical lines) by S-PLS for the reference traits.}
    \label{fig:Agnieszka1}
\end{figure}


After removing the missing values for each trait, a calibration set was created based on the training data using the Puchwein algorithm implemented in the \texttt{prospectr} package \citep{prospectr} with the remaining samples left out for validation. The Puchwein algorithm selects samples by iteratively eliminating similar samples using the Mahalanobis distance. Even if no outliers have been removed beforehand, they have not been included in the calibration subset by the algorithm.

Afterwards principal component analysis was considered for mean-centered and MSC-treated mean-centered spectra. As an exploratory insight, it has been observed that a higher number of principal components was needed to capture 95\% of variance for MSC-treated spectra.

Several regression models were tested, including Random forest, Support Vector Machines, and sparse PLS (S-PLS). 
The models used to provide the test dataset predictions were the ones with the lowest RMSECV and with the best fit when considering a regression of the predicted values versus the reference ones. For $\kappa$-casein and pH, the best model obtained was S-PLS, implemented in the \texttt{caret} package \citep{caret}, after removal of high-noise-levels (4 and 5 factors, respectively). For CMS the best results were obtained for the S-PLS model trained on the full spectrum after applying MSC (1 latent variable). The \texttt{caret} training included three resampling iterations with 5-folds cross validation. 
Wavenumbers selected by the S-PLS model for each trait were plotted over the spectra in Figure \ref{fig:Agnieszka1} (311, 145, and 266 variables for pH, $\kappa$-casein, and CMS respectively). A noticeable overlap in variables selected has been observed between models predicting pH and CMS.



\subsection{Other approaches}\label{sec:sec3_8}
In the previous sections the best strategies, in terms of the obtained results, have been outlined and discussed. Nonetheless some interesting approaches have been adopted also by other participants to the competition. In particular one of the participants detected the outliers by means of a multiple linear regression, fitted for each trait. Afterwards the outlying observations have been identified using the Bonferroni outlier test \citep[see e.g.][]{fox2015applied} and consequently deleted. Another researcher considered, when building the calibration models, the Interval Partial Least Squares model \citep[iPLS,][]{norgaard2000interval}. This is a graphically oriented method which can be seen as a local and interval-based extension of the standard PLS. As a consequence it might be particularly useful to interpret the results, since it highlights  different contributions from distinct spectral regions. Moreover a modification of the iPLS algorithm, called Backward Interval Partial Least Square \citep[biPLS,][]{norgaard2005multivariate}, has been tested. This approach works similarly to iPLS but it deletes sequentially the poorest performing interval. For an application of these algorithms on near-infrared data readers may refer to \citet{xiaobo2010genetic}.

\iffalse
\subsubsection{Uche}\label{sec:sec3_6}
As a first step, the missing values where removed for each one of the available reference traits. Afterwards Principal Component Analysis was performed on the full training dataset, for all the traits, and the number of principal components accounting for the 99\% of the total variability in the data were retained. For outlier detection multiple linear regression was fitted for each reference trait. 
Then, the outliers in the multiple linear regression response were detected using the Bonferroni outlier test \citep[see e.g.][]{fox2015applied} and consequently deleted. These pre-processing steps led to use 77, 78 and 78 principal components in the subsequent predictive analyses and resulted in a number of spectra retained in the training set equal to 393, 517 and 545 for $\kappa$-casein, CMS and pH respectively.

\iffalse
\begin{table}[h!]
\centering
    \begin{tabular}{ |l|c|c| } 
        \hline
        \multirow{2}{4em}{Response}&\multicolumn{2}{|c|}{Dimension} \\\cline{2-3}
         & No. Obs & PCs \\
        \hline
        $\kappa$ -casein & 393 & 77 \\ 
        CMS & 517 & 78 \\ 
        pH& 545 & 78 \\ 
        \hline
    \end{tabular}
\caption{Number of observations and number of principal components used for the final analysis}
\label{table:dim}
\end{table}
\fi

Four methods were considered for the prediction of the reference traits: multiple linear regression, support vector machine (SVM), random forest, and deep neural network with two layers. The training dataset has been divided randomly 100 times in training and validation set. The validation set contained 20\% of the training samples. 

The performances of the tested methods were compared using the RMSE computed on the validation set with SVM outperforming the other approaches for all the three traits. It has been therefore used to provide predictions for the test dataset. When using SVM a radial basis kernel has been employed while the value for epsilon was set to zero and the cost values used were 29 for $\kappa$-casein and 2 for both CMS and pH. All the hyperparameters of the model were tuned using a grid search.

\subsubsection{Elena approach}

First of all the missing values for all the available traits and the corresponding spectra was removed. Afterwards the training dataset was further split in order to obtain an internal validating set to be used for model accuracy assessment. Some commonly adopted pre-processing techniques, such as Standard Normal Variate and first and second derivative with Savitzky–Golay smoothing, were applied to the spectra but they were not considered for the final analyses as they did not seem to improve the model performances. 

For each reference trait different models were tested. As a first one Partial Least Squares regression (PLSR) has been used, with the optimal number of factors being selected via 10-fold CV. 
Moreover, iPLS and biPLS \textcolor{red}{need citations for both of them} were considered, with the wavenumbers divided into 20 and 30 intervals respectively and with a maximum of 6 factors considered for the latter approach \textcolor{red}{(don't know if iPLS and biPLS are that widespread. Probably they need some words explaining how they work, e.g. I do not know how the intervals are used and why there's a maximum number of factors + I suppose some tuning parameters -e.g. the optimal number of factors and intervals- have to be set. How did you do that? 10-fold cv here too?)}. These procedures have been successfully applied to improve PLSR regression models. \citep{xiaobo2010genetic}.

Models developed using PLSR provided poor predictions for $\kappa$ - casein, CMS, and pH. Slight improvements on the models were made using iPLS as visible in Table \ref{table:ElenaHayes}. iPLS and biPLS were used as variable selection methods to remove uninformative wavelengths from the spectrum to improve the predictive ability of the model. 

\begin{table}[t]
\centering
\caption{RMSECV averaged over the cross-validation folds for the methods tested.\textcolor{red}{(we need results from PLSR and standard deviation of the results over cv-folds)}}
\begin{tabular}{l|ccc}
  \hline
  Method/Trait & $\kappa$-casein  & CMS & pH \\ 
  \hline
   PLSR &  &  &   \\
  iPLS & 2.95 & 435.03 & 0.11 \\
  biPLS & 2.62 & 435.86 & 0.13 \\
   \hline
\end{tabular}
\label{table:ElenaHayes}
\end{table}


\iffalse
\begin{table}[t]
\centering
\caption{Statistics based on iPLS model prediction}
\begin{tabular}{l|ccc}
  \hline
  Trait & Interval no. & RMSE & RMSECV \\ 
  \hline
  $\kappa$-casein & 1 & 0.69 & 2.954 \\
  CMS & 1 & 218.54 & 435.031 \\
  pH & 1 & 0.111 & 0.115 \\
   \hline
\end{tabular}
\label{table:ElenaHayes}
\end{table}

\begin{table}[t]
\centering
\caption{Statistics based on biPLS model prediction}
\begin{tabular}{l|ccc}
  \hline
  Trait & Interval no. & RMSE & RMSECV \\ 
  \hline
  $\kappa$-casein & 10 &0.6812 &  \\
  CMS & 11 & 213.557 & 435.858 \\
  pH & 5 & 0.114 & 0.126 \\
   \hline
\end{tabular}
\label{table:ElenaHayes2}
\end{table}
\fi

In this study, biPLS proved to be the most effective procedure based on the lowest RMSECV and highest $R^{2}$ values in the validation set, and the results are reported in Table \ref{table:ElenaHayes2}. However, CMS was poorly predicted both in the training and validating set. Both $\kappa$ - casein and pH were predicted well in the validating set. The RMSE in the validation set increased significantly for $\kappa$ - casein indicating that the validation set was not well predicted. This information suggested that the model predicted on the training set was based on correlations found in the data but not relevant to the chemical composition or physical structure of the samples, therefore making the model inaccurate for unknown samples. 

\textcolor{red}{in order to be coherent across contributions I'd remove the RMSE for the training, considering only RMSECV. Anyway the last sentence might suggest that the splitting among training and validation was not performed correctly, more than a focus on ``wrong'' correlation. I need the results for PLSR and some more words about how iPLS and biPLS work in order to write the comments to the results and to highlight properly the advantages of using iPLS and biPLS with respect to plsr}
\fi

\section{Discussion}\label{sec:sec4}

As mentioned in Section \ref{sec:sec2}, the contributions have been evaluated considering the RMSEP for each single trait while, for the overall performances, the relative error as defined in (\ref{eq:rel_err}) has been computed. By thoroughly looking at the final results, reported in Table \ref{table:final_results}, some relevant insights can be obtained. First of all it is clear that the performances of different methods are usually trait specific. This sounds coherent with the so called \emph{no free lunch theorem} \citep[see e.g.][]{wolpert1997no} which provides an essential warning about the danger of comparing and ranking different algorithms by studying their performances on a small set of prediction problems. A general conclusion to draw is concerned with the poor predictability of the CMS, thus providing coherent indications with the ones obtained in previous works \citep[see e.g.][]{visentin2015prediction,frizzarin_predicting_2021}. Lastly note that different contributions performed comparably when predicting the test set values. 

A common choice that all the participants had to make was concerned with the possibility to transform the data prior to run their analyses. A first possible transformation consisted in considering absorbance instead of transmittance values. The well known \emph{Beer-Lambert law} seems to suggest that working on the absorbance scale can lead to improvements in the performances, even if sometimes the benefits are minimal. The results in Section \ref{sec:sec3_1}, where the participants tried both the approaches without noting significant improvements, as well as prior experiences of some of the authors, seem to suggest that this transformation is not strictly needed when analyzing MIRS data. 

Other transformations have been applied on the spectra by some of the participants, as a pre-treatment step. Procedures as the Multiplicative Scatter Correction or the Savitzky-Golay derivative smoothing have been tried but, again, without strong ameliorations of the results. Such transformations are usually hugely important when dealing with highly-noisy spectra, often encountered in the near-infrared spectroscopy applications. Nonetheless, in the dairy framework and when MIRS data are available, they seem to have a lower impact.

\begin{table}[t]
\centering
\caption{Root Mean Squared Error computed for the test set (RMSEP) for the reference traits and the overall relative error (RERR) as defined in \ref{eq:rel_err}. In bold, the best results are highlighted for each single trait and for the overall prediction. }
\begin{tabular}{l|ccc|c}
\hline
  Participant & & RMSEP & & RERR \\
   & $\kappa$-casein & CMS & pH &  \\ 
  \hline
  1 & 1.77 & 430.78 & {\bf 0.09} & 1.019 \\
  2 & 1.82 & 431.48 & 0.11 & 1.104 \\
  3 & 1.76 & 428.42 & 0.12 & 1.127 \\
  4 & {\bf 1.68} & 430.91 & {\bf 0.09} & {\bf 1.002} \\
  5 & 1.71 & {\bf 428.27} & 0.11 & 1.080 \\
  6 & 1.69 & 429.01 & 0.11 & 1.076 \\
   \hline
\end{tabular}
\label{table:final_results}
\end{table}

Another subjective decision that might have influenced the results was concerned with the outliers removal process. Outliers detection can be a critical step when working with spectroscopy data as it might be troublesome to discriminate between true extreme values and data collection errors, often due to the instruments used. In some of the contributions the analyses have been run considering all the available observations while some participants decided to follow the same approach considered in \citet{frizzarin_predicting_2021} where the observations have been considered as outliers when they had unusually large values for a reference trait; none of the participants considered procedures identifying outlying observations based on the spectra and therefore on unusual values of the wavenumbers. After the competition, the values for the traits on the test dataset were made publicly available. After re-running the analyses and re-evaluating their approach, a group of participants noted an improvement in the performances on the test set when outlying observations were not removed from the training set. This behavior might act as a partial confirmation of the challenges involved in the definition of the concept of outliers.

The participants were not provided with the specific information about the wavenumbers. Therefore, it was not straightforward to include prior information on the role of different wavenumbers in the modelling strategies. One of the possibly most impactful choices is concerned with the removal of the highly noisy regions related to water absorption processes, which are usually present when analyzing milk samples. In \citet{frizzarin_predicting_2021} these regions have been removed, according to the suggestions in \citet{hewavitharana1997fourier}. In the competition only two contributions considered the presence of the water absorption phenomenon: in fact in Section \ref{sec:sec3_4} this information has been used to devise an adaptive functional approach which provided good results both in terms of performances and in terms of interpretability. On the other hand, in Section \ref{sec:sec3_7}, the interquantile range for each wavenumbers has been used as a clever data-driven approach to identify noisy variables in the data. From previous experiences and from other analyses on similar data, it has been noted that the removal of highly noisy regions can be strongly influential and lead to improvements in the predictive performances. 

In the chemometrics framework, Partial Least Squares based methods still constitute the state of the art, as they are widely used in different applications usually providing good results. Nonetheless this challenge showed that it is worth looking at different modelling tools when analyzing spectroscopic data as they might slightly improve the quality of the obtained predictions or enhance the interpretability. In fact, one of the most interesting outcomes of the competition is related to the heterogeneity of the proposed approaches. Gathering together participants having different backgrounds, spanning from animal science to more statistical or machine learning ones, the contributions provided a wide overview about different possible ways to look at MIRS data. In most of the cases the spectra have been treated as standard vectors, also called \emph{tabular} approaches in Section \ref{sec:sec3_1}, thus disregarding the natural ordering of the wavenumbers. Nonetheless, the contributions in Section \ref{sec:sec3_2} and \ref{sec:sec3_4}, interestingly propose time series and functional approaches, respectively. The latter approach produced generally better predictions with respect to the first one; this might be due to the specific correlation structure usually encountered in MIRS data. In fact, even if neighboring wavenumbers tend to be highly correlated, strong relationships are often seen among spectral regions being far one from the other, since chemical constituents might have different absorption peaks at different spectral regions \citep[see][]{casa2021parsimonious}. Time series models, more often accounting for short term relationships and for more regular dependency patterns, might fail to properly capture this behaviour in spectroscopic data. On the other hand, the functional approach perfomed very well across the different traits, providing the lowest RERR in the challenge. Even if the spectra has already been treated before as continuous functions \citep{alsberg1993representation}, functional strategies are still not widespread in the chemometrics framework. Nonetheless, as pointed out in \citet{saeys2008potential} and as partially confirmed by the results of the competition, functional data analysis can find a fertile ground in infrared spectroscopy applications.  

A last relevant point, which have been discussed and faced both in the competition and during the whole workshop, is concerned with the adoption of methods that automatically perform wavelength selection. Often the participants advocated for the use of simpler methods which provide indications about the relevance of the features, therefore possibly providing more interpretable predictions from a practical standpoint. In fact such methods, jointly with proper collaborations with researchers with subject matter knowledge in the dairy framework, might lead to obtain information having a practical utility at different levels. For example from a commercial perspective, a proper identification of the most influential variables, can lead to the production of new spectrometers being cheaper and faster and possibly collecting data for those specific wavenumbers only. This can lead to a rapid and more pervasive diffusion of these technologies and, therefore, to an increase in the amount of data collected and to new challenges from a modelling standpoint. 

\section{Conclusion}\label{sec:conclusion}
The challenge organized during the first edition of the ``International Workshop on Spectroscopy and Chemometrics'' provided some insightful take home messages to both the participants and the researchers which attended to their presentations. First of all, the results showed how, when dealing with MIRS data, several different models provide good prediction and how their ranking can be strongly trait dependent. This advocates for the consideration of heterogeneous strategies when building calibration models in this framework. Moreover, it has been shown how different and less explored ways to treat spectral data, as for example adopting functional approaches, can lead to good results. Lastly, during the whole workshop, a lot of attention has been put on the necessity of relatively simple, parsimonious and interpretable models which can be more adequate when good predictions have to be complemented by the extraction of new knowledge about the phenomenon under study. 


\section*{Acknowledgements}
This publication has emanated from research conducted with the financial support of Science Foundation Ireland (SFI) and the Department of Agriculture, Food and Marine on behalf of the Government of Ireland under grant number (16/RC/3835), the SFI Insight Research Centre under grant number (SFI/12/RC/2289\_P2) and the SFI Starting Investigator Research Grant ``Infrared spectroscopy analysis of milk as a low-cost solution to identify efficient and profitable dairy cows'' (18/SIRG/5562). Agnieszka Konkolewska is funded from the European Union's Horizon 2020 program (Marie Sklodowska-Curie grant agreement No. 841882).

\bibliographystyle{apalike}

\subsection{Participant Katarina Domijan (KD)}

In the pre-processing steps the data in the training set were centered and scaled, with the same transformation being applied to the test set. Moreover CMS was log transformed, in order to reduce its skewness. Note that no outlier observations were removed and that the full set of 1060 wavenumbers was used for training and prediction. 

Hierarchical clustering with distance matrix based on the absolute value of correlation between the spectra and Ward's clustering criterion was used to find groupings of the spectra. The resulting dendrogram is displayed in Figure \ref{fig:spectradendro}; after a graphical inspection, the transmittance values were assigned to one of the fifteen clusters. For each cluster, the transmittance value with the highest correlation with the response variable was selected. This formed a set of fifteen features that has been used as predictors for one of the regression models described below. 

\begin{figure}[t]
    \centering
    \includegraphics[width=11cm, height = 7.5cm]{figs/spectradendro.png}
    \caption{Dendrogram of the spectra clustered by hierarchical clustering with distance matrix based on the absolute value of correlation between the spectra and Ward's clustering criterion.}
    \label{fig:spectradendro}
\end{figure}

The available training dataset with known trait values was initially split into 50 random splits of training (50\%) and validation (50\%) sets to select the optimal algorithm for regression. The algorithms were tuned using further cross-validation of the training subsets. 
All analysis were done in \texttt{R} and the code is available at \texttt{https://github.com/domijan/vistamilk}.
The following algorithms were used for regression:
\begin{itemize}
\item Lasso from library \texttt{glmnet} \citep{glmnet}.
\item Random forest (RF) from library \texttt{ranger} \citep{ranger} with and without regularization (RF + VS) \citep{regrf}.
\item Linear regression with the first six principal components (lm+PCA) as input features.
\item Linear regression with the first six kernel principal components (lm+kPCA) using Gaussian kernel implemented in library \texttt{kernlab} \citep{kernlab}.
\item Random forest with six kernel principal components (RF+kPCA) as input features.
\item Random forest with a Gaussian kernel in place of the transmittance values (RF+kernel) as input features.
\item Support vector machine (SVM) for regression implemented in library \texttt{e1071} \citep{e1071}.
\item Partial least squares (PLS) with 3, 4 and 5 components implemented in library \texttt{pls} \citep{pls}.
\item Linear regression with fifteen wavelengths (lm+15) selected from hierarchical clustering as described above. 
\item Bayesian additive regression trees \citep[BART,][]{bart} as in \texttt{bartMachine} library \citep{bartMachine}.
\item Post-hoc ensemble model that averaged over test set predictions from all of the models described above.
\end{itemize}


\begin{table}[t]
\centering
\caption{RMSE obtained over 50 random splits of the training set into training and validation for $\kappa$-casein, CMS  and pH. In brackets the standard deviations are reported.}
\begin{tabular}{l|ccc}
  \hline
  Algorithm & $\kappa$-casein & CMS & pH \\ 
  \hline
ensemble & 1.48 (0.21) & 367.5 (31.7) & 0.11 (<0.01) \\ 
lasso & 1.53 (0.21) & 368.3 (31.3) & 0.12  (<0.01) \\ 
lm+15 & 1.52 (0.18) & 366.2 (31.9) & 0.11  (0.01) \\ 
lm+PCA & 1.50 (0.20) & 367.9 (31.4) & 0.12  (<0.01) \\ 
lm+kPCA & 1.57 (0.19) & 368.2 (31.4) & 0.11  (<0.01) \\ 
PLS  & 1.50 (0.19) & 366.0 (31.4) & 0.12  (<0.01) \\ 
BART & 1.62 (0.12) & 368.3 (31.9) & 0.12 (0.01) \\ 
RF & 1.53 (0.18) & 367.9 (31.5) & 0.12  (0.01) \\ 
RF+VS & 1.52 (0.19) & 367.8 (31.6) & 0.12  (0.01) \\ 
RF+kPCA & 1.57 (0.21) & 369.4 (31.6) & 0.12  (<0.01) \\ 
RF+kernel & 1.53 (0.21) & 370.0 (31.6) & 0.12  (<0.01) \\ 
SVM & 1.54 (0.21) & 371.0 (31.3) & 0.11 (<0.01) \\ 
   \hline
\end{tabular}
\label{table:KD}
\end{table}

The predictive ability of the algorithms was evaluated according to the RMSE computed on the validation set. Table \ref{table:KD} displays the average RMSE over the 50 random splits with standard deviation in brackets. All the models considered had comparable performance; the best scoring algorithm for all three traits was the ensemble model, which slightly outperformed the other models. As a consequence it was selected and trained over the entire training set and used to predict the values for the test set. 






