Классификация рациона коров на основе спектров среднего инфракрасного диапазона молока: соревнование по анализу данных
на «Международном семинаре ⟦K0⟧
Спектроскопия и хемометрика 2022”

Данные были проанализированы с использованием различных стратегий моделирования, с акцентом как на методы, учитывающие упорядоченность длин волн, так и на методы, которые этого не делают. Все анализы были проведены преимущественно с использованием библиотек Python pandas, sklearn, sktime и matplotbib [see Pedregosa et al.,, 2011, and references therein]: код доступен по адресу https://github.com/mlgig/vistamilk_diet_challenge⟧.

В качестве первого шага были вычислены некоторые описательные статистики и удалены выбросы в соответствии как с рекомендациями, данными перед началом соревнования, так и с визуальным осмотром данных. На следующем этапе размеченный набор данных был разделен в соответствии со стратегией 3-кратной перекрестной проверки (3CV). Таким образом, лучшая модель была выбрана на основе точности перекрестной проверки, а затем обучена на полном обучающем наборе и использована для выполнения прогнозирования на предоставленном неразмеченном тестовом наборе..

Для прогнозирования рациона были рассмотрены следующие стратегии классификации:

• ⟦• Табличные модели: каждый образец рассматривается как вектор неупорядоченных признаков. В частности, были протестированы Ridge Classifier и Linear Discriminant Analysis (LDA). Далее эти методы были объединены как со стратегиями отбора признаков, так и со случайными полиномиальными преобразованиями признаков. Последний подход, основанный на генерации новых полиномиальных переменных из исходных, был направлен на проверку того, улучшают ли нелинейные взаимодействия точность классификации. В частности, представлен новый подход, целью которого является диверсификация полиномиальных признаков при сохранении низких вычислительных требований..
• ⟦• Модели глубоких нейронных сетей: было протестировано семейство подходов на основе глубоких нейронных сетей, как полносвязных, так и сверточных. Эта стратегия неявно генерирует сложные взаимодействия признаков, заложенные в архитектуре сети.

Обратите внимание, что ранее полученные результаты [Frizzarin et al., 2021a, ] показывают, что табличные методы работают довольно хорошо с данными спектроскопии. Более того, следуя рекомендациям в Frizzarin et al., 2021b [], стратегии отбора признаков были совмещены с информацией о наличии областей воды в спектрах. Кроме того, были протестированы современные алгоритмы классификации временных рядов, такие как ROCKET [Dempster et al.,, 2020], MiniROCKET [Dempster et al.,, 2021], MrSQM [Nguyen and Ifrim,, 2021; Nguyen and Ifrim,, 2022] и FreshPrince [Middlehurst and Bagnall,, 2022]. Наконец, были применены ансамблевые методы, направленные на объединение моделей временных рядов и табличных моделей для совмещения их предсказаний и сильных сторон. Тем не менее, эти подходы уступили по эффективности вышеупомянутым, поэтому соответствующие результаты не приводятся в следующих разделах..

Все этапы обработки и реализация алгоритма были выполнены с использованием MATLAB [MATLAB,, 2018]. После импорта набора данных в табличной форме выбросы были определены как наблюдения, имеющие более трех масштабированных медианных абсолютных отклонений от медианы набора данных. Классификация была выполнена с использованием набора алгоритмов, таких как Support Vector Machine (SVM), K-Nearest Neighbors (KNN) и Linear Discriminant Analysis (LDA). Настройка гиперпараметров и оценка точности классификации были выполнены посредством 5-кратной перекрестной проверки..

Наилучшие результаты были получены с использованием LDA, которая смогла отличить молоко коров на пастбищном вскармливании от TMR с точностью 95%, в то время как различение травы и клевера происходило с точностью 68%. Рисунок 3 позволяет визуализировать границы классов путем построения проекций спектров в латентном пространстве, охватываемом двумя дискриминантными функциями. На рисунке видна четкая граница между классами стойлового и пастбищного содержания, в то время как между классами GRS и CLV наблюдается значительное перекрытие. Следовательно, извлеченные компоненты затем рассматривались как входные данные для линейной модели SVM для улучшения классификации между классами пастбищного содержания. Комбинация двух классификаторов (LDA + SVM), представляющая собой двухэтапный подход, значительно улучшила общую точность классификации (87,1%), а также точность классификации между классами, как показано в Таблице 5.

Настоящая работа была выполнена независимо тремя участниками группы после проведения общего предварительного анализа спектральных данных. Результаты прогнозирования на тестовой выборке, предоставленной для хемометрического конкурса, были затем сопоставлены для оценки согласованности между тремя различными примененными статистическими подходами..

Использовался стандартный конвейер машинного обучения, состоящий из выбора признаков (т.е. длин волн) и классификации, при этом из исходного набора данных не удалялись выбросы. Методы снижения размерности, такие как метод главных компонент (PCA) и метод независимых компонент (ICA), а также расширенная мультипликативная коррекция рассеяния (EMSC) и подход к аугментации данных были протестированы для улучшения результатов классификации [Bjerrum et al.,, 2017]. EMSC представляет собой метод предварительной обработки, который устраняет мультипликативные эффекты, потенциально вызванные физическими явлениями, такими как рассеяние света, что обычно наблюдается в спектроскопии отражения, тем самым упрощая моделирование химических эффектов. С другой стороны, схема аугментации данных увеличивает набор данных в десять раз путем девятикратного добавления случайных вариаций смещения, умножения и наклона к каждому образцу..

Впоследствии был использован ряд различных классификаторов, которые ранее успешно применялись к данным инфракрасной спектроскопии. В частности, рассматривались следующие модели: K-nearest Neighbour [Balabin and Safieva,, 2011, K-NN;], Random Forest [Chen et al.,, 2021, RF;], Support Vector Classification [Ji-yong et al.,, 2013, SVC;], Multilayer-perceptron [Balabin and Safieva,, 2008, MLP;], Linear Discriminant Analysis [Khuwijitjaru et al.,, 2020, LDA;], Decision Tree Classification [Geronimo et al.,, 2019], Nu-Support Vector Classification [Terouzi et al.,, 2013, NuSVC;], AdaBoost Classification [Wu et al.,, 2017], Gradient Boosting Classification [Munera et al.,, 2021], Gaussian Naive Bayes [Bhati and Bhattacharya,, 2020] и Quadratic Discriminant Analysis [Oravec et al.,, 2019, QDA;]. Другие исследованные методы прогнозирования относились к группе методов deep Learning (DL), в частности к одномерным (1D) Convolutional Neural Network (CNN). 1D CNN использует шесть одномерных сверточных слоев, а также ряд слоев max pooling, batch normalization и dropout. За каждым слоем 1D CNN следует слой max-pooling и batch normalization. One-dimensional CNN применялась только к исходным спектрам с целью сохранения последовательности данных, что не требуется для PCA и ICA.

Перед проведением анализов набор данных был разделен на обучающую выборку (80% данных) для обучения различных моделей и валидационную выборку (оставшиеся 20% данных), используемую для оптимизации гиперпараметров и определения наилучшего метода для финального тестирования..

Первоначальный эксперимент был проведен на всех классификаторах без использования аугментации данных или отбора признаков. Это было выполнено для изучения того, какой метод классификации показывает лучшие результаты при работе с необработанными спектральными данными. На рисунке 5 представлены результаты, полученные на начальном этапе с использованием разбиения 80/20 для обучения/валидации для различных классификаторов. Все полученные результаты представляют собой средние значения, рассчитанные по трем прогнозам обучения и валидации для каждой модели. LDA показал наилучшие результаты с точностью 76%, в то время как MLP и SVC продемонстрировали одни из худших показателей с точностью около 33%.

На втором этапе классификаторы тестировались в сочетании с PCA, ICA (использовались методы PCA и FastICA из scikit-learn [Pedregosa et al.,, 2011]) или аугментацией данных. Использование PCA и ICA изменяло данные путем уменьшения размерности, в то время как аугментация данных увеличивала количество образцов. Для аугментации данных использовалась функция data augment из Bjerrum et al., [2017]. Это увеличило количество обучающих образцов с 3 244 до 19 464. На данном этапе рассматривалось только подмножество ранее протестированных моделей, основываясь на их эффективности на предыдущем шаге. На рисунке 6 представлены результаты каждого классификатора с каждым методом предварительной обработки (базовый, ICA, PCA, аугментация данных (Aug)). Исходя из этих результатов, было отмечено, что LDA после аугментации данных достиг наивысшей точности в 82,7%. Наибольшее улучшение прогнозов наблюдалось при использовании MLP после ICA (улучшение на 41%). Затем был проведен дополнительный эксперимент только с использованием модели LDA. Это было сделано для того, чтобы показать важность областей внутри спектров, и было протестировано несколько различных диапазонов длин волн. Таким образом, на рисунке 7 показаны результаты LDA при удалении различных спектральных областей..

При использовании аугментации данных наблюдалось общее повышение точности по сравнению с базовым подходом, за единственным исключением в виде CNN. Что касается выбора длин волн, заметного увеличения точности при фокусировке на конкретной области спектра не наблюдалось. Тем не менее, большая часть релевантной информации находилась в диапазоне от 925 см-1 до 1597 см-1, и наблюдалось небольшое увеличение точности прогнозирования примерно на 1% при использовании диапазонов от 925 до 1585 см-1 и от 1717 до 2103 см-1 по сравнению с полным набором длин волн..

Для подготовки набора данных к прогностическому анализу была проведена предварительная обработка. Согласно указаниям организаторов конкурса, выбросы в спектрах были удалены, в результате чего набор данных составил 3243 спектра пропускания, охватывающих 1060 длин волн. Впоследствии спектры были преобразованы в значения поглощения путем взятия $\log_{10}$ от величины, обратной значениям пропускания. Кроме того, следуя Frizzarin et al., 2021b [], подмножество из 534 длин волн, которые находились вне областей с высоким уровнем шума, связанных с водой, было определено как значимое для прогнозирования рациона коровы, хотя области, связанные с водой, на данном этапе анализа исключены не были..

Для обеспечения надежной оценки набор данных был разделен на обучающую и валидационную выборки. В данном случае валидационная выборка была сформирована для контроля влияния эффекта партии, который может искажать оценки при прогнозировании на новых данных [Soneson et al.,, 2014]. Анализ набора данных показал, что строки были упорядочены таким образом, что содержали несколько последовательных наблюдений для каждого рациона. Следовательно, было сделано допущение, что каждый набор последовательных наблюдений рациона относится к одной партии. Таким образом, было идентифицировано 90 партий, по 30 для каждого рациона. Кроме того, данные собирались в течение трех лет [Frizzarin et al., 2021b, ], поэтому было принято, что первые 30 партий были собраны в первый год исследования, следующие 30 — во второй год, а последние 30 — в третий. Исходя из этих допущений, валидационная выборка состояла из 996 спектров из 30 партий, собранных в третий год исследования, что включало десять партий для каждого рациона, в то время как модели обучались на оставшихся 2247 спектрах. Обучающие данные были случайным образом разделены на $V=10$ фолдов, при этом каждый фолд включал две партии для каждого рациона. Возможный эффект партии от повторных измерений для одной коровы не учитывался..

Для описания прогностической модели, используемой в данном анализе, пусть $\mathcal{D}=\left\{y_{i},\mathbf{x}_{i}\right\}_{i=1}^{N}$ обозначает наблюдаемые данные, где зависимая переменная $y_{i}\in\left\{1,\dots,M\right\}$ представляет рацион $i$-й коровы, а ковариаты $\mathbf{x}_{i}\in\mathbb{R}^{D}$ представляют соответствующий спектр поглощения молока. Заметим, что в данном анализе рассматриваются $M=3$ рационов, $D=1060$ длин волн и $N=3243$ обучающих наблюдений. Целью предлагаемых прогностических моделей является изучение $\mathbb{P}\left(y\mid\mathbf{x}\right)$, то есть вероятности того, что данный образец молока получен от коровы, питавшейся травой, клевером или TMR, при заданном спектре для этого образца..

Первым шагом при построении прогностической модели является определение детерминированной функции отображения $g:\mathbf{x}_{i}\to\mathbf{z}_{i}$ для $\mathbf{z}_{i}\in\mathbb{R}^{D^{\prime}}$ с $D^{\prime}<D$, которая описывает процедуру извлечения признаков. Здесь были рассмотрены два подхода к извлечению признаков. Первый заключается в простом выборе $D^{\prime}=534$ релевантных длин волн, определенных с помощью Frizzarin et al., 2021b [], таким образом, что $\mathbf{z}_{i}$ является $i$-м спектром поглощения после удаления областей воды с высоким уровнем шума и стандартизации каждой длины волны. Второй подход основан на вейвлет-преобразовании — популярном методе обработки сигналов, который может применяться для сжатия данных, сглаживания и многомасштабного анализа [Nason,, 2008], и выполняется в три этапа. После установки значений областей с высоким уровнем шума каждого спектра на 0, пороговое вейвлет-преобразование предоставляет набор вейвлет-коэффициентов. Вектор признаков $\mathbf{z}_{i}$ представляет собой вектор вейвлет-коэффициентов, которые являются ненулевыми по крайней мере для одного из $N$ спектров, в данном случае $D^{\prime}=594$. Пороговое вейвлет-преобразование доступно в R-пакете wavethresh [Nason,, 2016] с использованием наименее симметричного вейвлета Добеши в качестве материнского вейвлета и байесовского подхода к пороговой обработке вейвлет-коэффициентов [Abramovich et al.,, 1998]. Обратите внимание, что установка длин волн в областях с высоким уровнем шума на 0 означает, что вейвлет-преобразование сохраняет спектральное расстояние между длинами волн, гарантируя при этом, что соответствующие вейвлет-коэффициенты равны 0.

Учитывая вектор признаков $\mathbf{z}_{i}=g\left(\mathbf{x}_{i}\right)$, была принята мультиномиальная регрессионная модель для рациона, такая что

для $m=1,\dots M$, где $\bm{\beta}_{m}\in\mathbb{R}^{D^{\prime}}$, неявно предполагая, что $\mathbf{z}_{i}$ включает член пересечения. Пакет glmnet [Friedman et al.,, 2010] эффективно подгоняет эту модель под данные. Для простоты была подогнана модель LASSO, в которой 10-кратная перекрестная проверка на обучающих данных определяет гиперпараметр штрафа.

Наконец, прогностическая эффективность предложенных моделей была сравнена путем анализа их логарифмической функции потерь на валидационном наборе данных. То есть, для валидационного набора данных и модели $\mathcal{M}_{j}$ для $\mathbf{z}_{i}=g\left(\mathbf{x}_{i}\right)$, логарифмическая функция потерь определяется как

где $N^{\prime}$ — количество наблюдений в валидационном наборе, $\mathbb{I}\left(\mathcal{A}\right)$ — обычная индикаторная функция, которая равна 1, когда $\mathcal{A}$ истинно, и 0 в противном случае, а $\mathbb{P}\left(y_{i}=m\mid\mathbf{z}_{i},\mathcal{M}_{j}\right)$ — вероятность при $\mathcal{M}_{j}$ того, что $y_{i}=m$ при условии $\mathbf{z}_{i}$. Логарифмическая функция потерь является надлежащим правилом оценки для анализа прогностических моделей [Gneiting and Raftery,, 2007], где меньшие значения являются лучшими, что стимулирует выражение истинного представления о данных. Также довольно просто установить ориентиры для оценки качества прогнозов априори. Например, для любого $M$ средняя логарифмическая функция потерь, равная 0, представляет собой идеальную прогностическую эффективность, в то время как, когда $M=3$, как в рассматриваемом случае, средняя логарифмическая функция потерь, равная $-\ln(1/3)\approx 1.1$, представляет собой «угадывание», при котором мы прогнозируем каждую категорию равномерно случайным образом. Для полноты картины также была оценена точность классификации $\mathcal{M}_{j}$.

Результаты этого анализа представлены в Таблице 7. Первой рассмотренной моделью была мультиномиальная регрессия с LASSO-штрафом для спектров сырого молока в зависимости от рационов питания, в которой области спектра с высоким уровнем шума были исключены, а длины волн стандартизованы. Настраиваемый параметр $\lambda$, управляющий силой штрафа, был выбран путем минимизации мультиномиальной девианты (статистического показателя, пропорционального среднему значению log-loss) с помощью 10-кратной перекрестной проверки. Значение log-loss этой модели на обучающем наборе составляло 0,57, что соответствует точности классификации рациона в 77%. Более детальный анализ прогнозов показал, что при объединении CLV и GRS в одну категорию (кормление на пастбище) можно было прогнозировать TMR с точностью 94%. При попытке предсказать, получала ли корова корм CLV, при условии ее пастбищного содержания, была достигнута точность 72%. Прогностические характеристики на валидационном наборе оказались значительно хуже: общий показатель log-loss составил 0,82, что соответствует точности 58%. Модель прогнозировала TMR с точностью 88%. Однако для коров, которые, как известно, содержались на пастбище, она предсказывала CLV с точностью 49%.

Во второй модели рассматривалась мультиномиальная регрессия ненулевых пороговых коэффициентов вейвлет-преобразования спектров молока в зависимости от рационов. Как и выше, модель подгонялась путем максимизации пенализированного логарифмического правдоподобия с использованием 10-кратной перекрестной проверки для настройки $\lambda$. Для этой модели логарифмические потери на обучающем наборе были равны 0.74, что соответствует точности 69%, хотя она предсказывала TMR с точностью 88%. Для коров пастбищного содержания она предсказывала CLV с точностью 68%. Как и в случае с первой моделью, производительность снизилась для валидационного набора. Логарифмические потери составили 0.88, а точность TMR — 79%. Учитывая, что корова находилась на пастбищном содержании, точность CLV составила 47%. Эти результаты суммированы в Таблице 7.

Полученные результаты наглядно показали, что спектры молока содержат сигнал, позволяющий отличать коров на пастбищном вскармливании от TMR, однако различить CLV и GRS оказалось затруднительно. Тем не менее, прогностическая эффективность на валидационном наборе данных оказалась значительно хуже, чем на обучающем, что указывает на то, что принятые модели не обеспечивали надежных прогнозов вне выборки. Без тщательного учета потенциальных искажающих факторов батч-эффекта в исследованных спектрах мы, вероятно, переоценим производительность наших моделей вне выборки. Сбор данных от большего числа коров в течение более длительного периода времени должен решить эту проблему и позволить разработать более надежные модели.

Наконец, не было обнаружено никаких свидетельств того, что спектры после вейвлет-преобразования предоставили полезную информацию о рационе коров. Тем не менее, это не значит, что какое-то альтернативное разложение базиса не могло бы улучшить текущие предиктивные модели. Фактически, учитывая больше данных о взаимосвязи между спектрами молока и рационом, разработка моделей, допускающих нелинейные зависимости между длинами волн, может оказаться перспективным направлением для будущих исследований..

В качестве первого шага обучающая выборка была центрирована и масштабирована, и та же трансформация была применена к тестовой выборке. В последующих анализах выбросы не удалялись, в то время как все спектры были преобразованы из коэффициента пропускания в оптическую плотность. Длины волн из спектральных областей с высоким уровнем шума между 1720 и 1592 см-1, между 3698 и 2996 см-1 и более 3818 см-1 были исключены из анализа в соответствии с Frizzarin et al., 2021b [].

Оценка Фишера, представляющая собой отношение межгрупповой дисперсии к внутригрупповой дисперсии рациона, была рассчитана для всех длин волн в обучающей выборке. Для длины волны $j$ оценка Фишера определяется как:

где $j$ обозначает индекс длины волны, $i=1,\dots,n$ обозначает спектры, причем $n$ является количеством спектров в обучающем наборе, $m$ обозначает группу рационов с $M=3$, $\mathbb{I}(y_{i}=m)$ является индикатором спектров группы рационов $i$, $\bar{x}_{\cdot j}$ представляет собой среднее значение длины волны $j$ для всех спектров $(i=1,\dots,n)$, $\bar{x}_{\cdot j}^{(m)}$ представляет собой среднее значение длины волны $j$ в группе рационов $m$. Длина волны с наибольшим критерием Фишера в каждой из отброшенных областей была сохранена в анализе. Длины волн с критерием Фишера ниже 0,002 были исключены из дальнейшего анализа, в результате чего осталось 380 длин волн. Для сравнения алгоритмов и выполнения дальнейшего отбора признаков обучающий набор был случайным образом разделен в соотношении 75/25 на обучающую и тестовую выборки с стратификацией по рациону. Генетический алгоритм [Holland,, 1992], реализованный в библиотеке genalg [Willighagen and Ballings,, 2022], использовался в качестве метода стохастического поиска для нахождения оптимального подмножества входных длин волн для классификации. Осуществлялось представление особей в популяции GA в виде бинарных строк, обозначающих длины волн, которые должны быть включены или исключены для прогнозирования. Целевая функция была задана как средняя точность по результатам десяти кросс-валидированных подгонок линейного дискриминантного анализа (LDA) обучающего подмножества. GA выполнялся в течение 200 итераций при размере популяции 200. На рисунке 8 показаны оптическая плотность спектров и соответствующие им значения критерия Фишера, причем точки обозначают длины волн, выбранные с помощью GA.

Наилучшая конфигурация из финальной популяции GA включала 70 длин волн. Эти длины волн использовались в качестве входных данных для следующих алгоритмов классификации:

• ⟦• Линейный дискриминантный анализ (LDA), библиотека MASS [Venables and Ripley,, 2002]⟧;
• • Дискриминантный анализ на основе частичных наименьших квадратов (PLS-DA) [Mevik et al.,, 2020];
• ⟦• Метод наименьшего абсолютного сужения и выбора [Tibshirani,, 1996, LASSO;], библиотека glmnet [Friedman et al.,, 2010]⟧;
• ⟦• Elastic net [Zou and Hastie,, 2005, EN;], библиотека glmnet;
• ⟦• Random Forest [Breiman,, 2001, RF;], библиотека ranger [Wright and Ziegler,, 2017]⟧;
• ⟦• Метод опорных векторов [Vapnik,, 1998], библиотека kernlab [Karatzoglou et al.,, 2004]⟧;
• • Байесовский классификатор на основе ядерной проекции [Domijan and Wilson,, 2011, BKPC], библиотека BKPC [Domijan,, 2018].

Все анализы были выполнены с использованием R [R Core Team,, 2020], код доступен в репозитории github https://github.com/domijan/KD_Vistamilk2022.

Обучающая выборка была случайным образом разделена на десять дополнительных обучающих/тестовых наборов равного размера, стратифицированных по рациону. Средняя точность и стандартное отклонение по десяти случайным разбиениям для всех алгоритмов классификации приведены в таблице 8. Наилучшие показатели продемонстрировал LDA со средней точностью 77,4%. Общая точность PLS-DA и EN составила 76,9% и 76,5% соответственно. Алгоритмы были настроены с использованием дополнительной кросс-валидации обучающих наборов. Для BKPC и SVM наилучшие результаты были получены с линейным ядром. Прогнозы LDA были представлены на конкурс. Более того, генетический алгоритм позволил отобрать значительно меньшее подмножество длин волн без потери эффективности классификации..