Проекционные методы в линейном регрессионном анализе: РГК/ПЛС Андрей Юрьевич Богомолов Российское хемометрическое общество European Molecular Biology Laboratory (EMBL) «Введение в анализ многомерных данных» (школа WSC-5), 16 февраля 2006, Самара

Андрей Юрьевич Богомолов Российское хемометрическое общество European Molecular Biology Laboratory (EMBL) «Введение в анализ многомерных данных» (школа WSC-5), 16 февраля 2006, Самара Методы многомерной калибровки

Тема лекции Многомерная калибровка Multivariate Calibration Анализ многомерных данных (Хемометрика) Multivariate Data Analysis (Chemometrics)

К вопросу о русской терминологии родной язык хемометрики - английский терминология за 30 лет устоялась: статьи, учебники, книги, конференции привычные аббревиатуры: PCA, PCR, PLS, SIMCA, RMSEP, etc. - не нуждаются в расшифровке русская терминология создается сейчас нужен ли перевод? – да! например: scores and loadings (!?) нужно время, чтобы русские термины вошли в обиход в настоящей лекции - параллельная терминология

Калибровка или градуировка? в русском языке – два сходных термина: « КАЛИБРОВКА (средств измерений) – совокупность операций, выполняемых с целью определения и подтверждения действительных значений метрологических характеристик и (или) пригодности к применению средств измерений…» « ГРАДУИРОВКА – метрологическая операция, при помощи которой устанавливается значение меры или делениям шкалы измерительного прибора придаются значения...» на английский оба переводятся как calibration «градуировка» – официальный термин в лекции будет использоваться некорректный термин «калибровка»

Регрессия & Калибровка Regression is an approach for relating two sets of variables to each other Kim Esbensen Calibration is a process of constructing a mathematical model to relate the output of an instrument to properties of samples Kenneth Beebe Калибровка ~ Регрессия

Регрессионный анализ XY ? зависимая переменная независимая переменная линейная регрессия Y = XB + E МГК (PCA) – моделирование (X) регрессия – моделирование (X,Y)

Спектральные данные Спектры (X) Концентрации (Y) [C 1 ] … [C q ]

Для чего нужна калибровка? замена прямого измерения интересующего свойства, измерением другого, коррелирующего с первым такая потребность возникает если прямое измерение интересующего свойства нежелательно: дорого трудоемко занимает много времени этически нежелательно эксперимент невозможен, и т. п. в подавляющем числе практических ситуаций такая замена оправдана!

Примеры из различных областей ХИМИЯ: калибровка – инструмент 1 количественного анализа БИОЛОГИЯ: непосредственный анализ может быть губителен для живых существ МЕДИЦИНА: неинвазивный анализ, например, определение сахара в крови спектроскопически (ближний ИК) ПСИХОЛОГИЯ: анализ личности может потребовать длительных наблюдений, желательно использовать косвенные данные СОЦИОЛОГИЯ и ФИНАНСЫ: предсказание может быть основано только на исторических данных

Одномерная калибровка: один компонент univariate calibration y = x × r =

двухкомпонентная смесь r= Одномерная калибровка: многокомпонентная смесь компоненты r=0.9988

Многомерная калибровка y=xb+eY=XB+E

Преимущества многомерной калибровки возможность анализировать несколько компонентов одновременно выигрыш в точности от усреднения при использования «избыточных», в т.ч. сильно коррелирующих измерений (спектры) возможность диагностики «плохих» образцов в процессе предсказания «парадигматический сдвиг» в подходах к решению проблем с появлением ПЛС регрессии (PLS-R) спектроскопия ближнего ИК стала одним из наиболее популярных методов анализа

Калибровка и предсказание Калибровка (Calibration) Предсказание (Prediction)

Классические и инверсные методы Два основных подхода в многомерной калибровке: Классический МНК (Classical Least Squares, CLS) основан на прямом решении уравнения Бугера- Ламберта-Бера A = Cε | X = Yε Инверсный МНК (Inverse Least Squares, ILS) решают уравнение вида С = Ab | Y = Xb В настоящей лекции – только ILS

Множественная линейная регрессия (МЛР) Multiple Linear Regression (MLR) Решение: b = (X T X) -1 X T y n - число объектов (спектров) p - число переменных (длин волн) y=b 0 + b 1 x 1 + b 2 x 2 +…+b p x p +e

Недостатки МЛР МЛР может не сработать, если: высока коллинеарность в X (спектры) неустойчивое решение для коллинеарных даных обусловлено преобразованием (X T X) -1 X T высокий уровень шума, ошибки в X переменных больше, чем образцов (типично для спектральных данных) есть линейная зависимость между переменными внутри X визуальная интерпретация МЛР-моделей затруднительна

Пример спектральных данных: полиароматические углеводороды λ, нм ε, M -1 см C × ε + E = D R.S.D. (E) = C = C 0 + E c R.S.D. (E c ) = 5% (C max ) nКомпонент[C n ], M 12-ацетофенантрен ацетиламинофенантрен ацетиламинофенантрен

Полиароматические углеводороды: обучающий и тестовый наборы Обучающий набор Training Set Тестовый набор Validation Set «simdata»

МЛР-калибровка (Simdata) [C 1 ][C 2 ] [C 3 ] точность МЛР-модели для [С 3 ] (3-го компонента смеси ПАУ) неудовлетворительна

МГК (Principle Component Analysis) - преобразование: X = TP T + E счета T (scores) и нагрузки P (loadings) определяют пространство клавных компонент T ортогональны и содержит проекции данных на ГК Метод Главных Компонент (МГК) - оружие против коллинеарности = + X (исходные) r c a rr с E (шум) X (оценка) P (нагрузки) T (счета) c T можно использовать вместо X для анализа (!)

Концепция PCA «на пальцах» X=A(522 nm) Y=A(644 nm) Z=A(714 nm) X=A(430 nm) Y=A(550 nm) Z=A(750 nm)

МГК + МЛР = РГК! (PCA + MLR = PCR) МГК-счета (PCA scores) T можно использовать вместо X для построения МЛР-модели (MLR): MLR: y=Xb+e | b=[XX T ] -1 X T y | y new =X new b (I) PCR: y=Tb+e | b=[TT T ] -1 T T y | y new =T new b (II) Метод называется: регрессия на главные компоненты, РГК (Principal Component Regression, PCR)

Схема РГК (PCR) – подробнее PCA: MLR: n - объектов p - переменных a - главных компонент

Интерпретация РГК-модели интерпретация модели служит для изучения внутренней структуры данных: группы выбросы связь между X и Y инструменты диагностики МГК (PCA) работают в РГК (PCR): график счетов (scores) график нагрузок (loadings) график счетов и нагрузок вместе (bi-plot) график остатков (residuals) инструменты диагностики РГК: совместный график нагрузок X и Y

Строим РГК-модель (Simdata)

Строим РГК-модель (simdata)

Проверка (валидация) модели проверка (validation) модели служит для: определения размерности модели (числа ГК) оценки предсказательной способности модели проверка модели производится с помощью тестовых данных: того же диапазона и того же качества что обущающие данные (та же генеральная выборка) достаточно представительные или кросс-валидации (cross-validation) полная (leave-one-out, LOO) сегментная (например, Venetian blind)

Среднеквадратичная ошибка предсказания (RMSEP) RMSEС = Root Mean Square Error of Calibration RMSEP = Root Mean Square Error of Prediction минимум на кривой RMSEP – основной индикатор числа ГК RMSEP – оценка точности в единицах измерения (!) RMSEP используется для сравнения моделей

Число компонент: почему минимум на кривой RMSEP? включенная ошибка остаточная информация

Оценка числа компонент в РГК правильный выбор числа главных компонент (principle components, PC) - ключевая проблема многомерной калибровки модель с недостаточным числом ГК (underfitting) не использует всей полезной информации из данных модель с избыточным числом ГК (overfitting) начинает моделировать шум (ошибку) найти оптимальную размерность помогают тестовые данные (validation set)

Число компонент: РГК - simdata

Оценка числа ГК в РГК: особенности число главных компонент (размерность модели) определяется в РГК (PCR) нуждами калибровки, и не обязательно совпадает с результатом МГК (PCA) Особенности: в РГК есть RMSEP активно используется тестовые данные (test set) минимум на кривой RMSEP - основной индикатор числа ГК для спектральных данных показательной может быть форма X-нагрузок (X-loadings) решение всегда за экспертом!

Несовершенства РГК РГК (PCR) – мощный метод многомерной калибровки имеет безусловные преимущества перед MLR однако, не вполне оптимизирован для калибровки пространство ГК не учитывает структуры Y и связи между X и Y можно ли учесть эту связь при построении проекционной модели? да, это делает PLS!

Факторные пространства уравнение PCA имеет универсальный смысл: X = TP T + E преобразование называется факторной компрессией, проекцией данных на факторное пространство (factor space) парные вектора в T и P называются факторами (factors) главные компоненты – важный пример факторного пространства, но не единственный факторное пространство можно оптимизировать для решения конкретной задачи ГК (PC) оптимальны для исследования структуры X как оптимизировать пространство для калибровки?

PLS – мощная альтернатива PCR Метод проекции на латентные структуры (ПЛС) и ПЛС-регрессия (ПЛС-Р) PLS = Partial Least Squares -> = Projection on Latent Structures ПЛС-пространство создается при участии двух переменных X и Y одновременно критерий – моделирование той структуры (информации) в X, которая коррелирует с Y например, спектральные полосы (X), которые отвечают за концентрацию компонента(ов), заданные в Y, получат в подели больший вес метод ПЛС оптимизирован для регрессионного анализа

ПЛС-регрессия: схематическое представление участвуют обе матрицы X и Y факторы рассчитываются по очереди – алгоритм NIPALS => 2 набора счетов (scores) T, U и нагрузок (loadings) P, Q плюс матрица W взвешенных нагрузок (loading-weights) итерационное улучшение модели, чтобы максимизировать cov(T,U) Предсказание: Ŷ = T new B t Ŷ = X new B B = W(P T W) -1 Q T T X W P U Y Q X = TP T + E x Y = UQ T + E y [1] S. Wold, H. Martens, H. Wold, Lecture Notes Math. 973 (1983) 286–293

Две разновидности ПЛС: ПЛС1 и ПЛС2 существуют две популярных разновидности ПЛС: ПЛС1 (PLS1) и ПЛС2 (PLS2) ПЛС1 модель строится для единственной переменной y (свойства), например, для концентрации одного компонента смеси если нужна калибровка по нескольким свойствам, строится несколько независимых моделей ПЛС2 рассчитывается для нескольких свойств одновременно расчетные алгоритмы методов отличаются соответственно

Основы алгоритма ПЛС ПЛС-декомпозиция производится алгоримом NIPALS NIPALS = Non-linear Iterative Partial Least Squares факторы находятся по очереди, один за другим, расчет всех факторов (как в SVD) не обязателен итерационная замена векторов u f -> t f и u f -> t f для нахождения текущего фактора f - алгоритмическая основа ПЛС2 алгоритм работает до выполнения критерия сходимости ознакомимся с принципиальной схемой, начиная с более общего ПЛС2

NIPALS алгоритм для ПЛС2 0.0.ufuf выбор начального приближения u 1.1.w f = X f T u f /|X f T u f |расчет нормализованного вектора взвешeнных нагрузок w 2.2.t f = X f T w f расчет вектора весов t 3.3.q f = Y f T t f /|Y f T t f |расчет нормализованного вектора нагрузок q 4.4.u f = Y f T q f расчет вектора счетов u 5.5.|t f.new - t f.old |< lim?проверка сходимости: да -> go to 1. 6.p f = X f T t f /t f T t f расчет вектора весов p 7.7.b f = u f T t f /t f T t f расчет внутненнего коэффициента регрессии b 8.8.X f+1 = X f+1 - t f T p f Y f+1 = Y f+1 - b f t f q f T расчет остатка X и Y 9.f = f + 1f = f + 1Переход к следующему фактору

NIPALS алгоритм для ПЛС1 1.1.w f = X f T y f /|X f T y f |расчет нормализованного вектора взвешeнных нагрузок w 2.2.t f = X f T w f расчет вектора весов t 3.3.q f = y f T t f /|t f T t f |расчет нагрузки q (скаляр) фактора f 4.4.p f = X f T t f /t f T t f расчет вектора весов p 5.5.X f+1 = X f+1 - t f T p f y f+1 = y f+1 - q f t f расчет остатка X и y 6.6.f = f + 1f = f + 1переход к следующему фактору

NIPALS алгоритм для ПЛС1 1.1.w f = X f T y f /|X f T y f |расчет нормализованного вектора взвешeнных нагрузок w 2.2.t f = X f T w f расчет вектора весов t 3.3.q f = y f T t f /|t f T t f |расчет нагрузки q (скаляр) фактора f 4.4.p f = X f T t f /t f T t f расчет вектора весов p 5.5.X f+1 = X f+1 - t f T p f y f+1 = y f+1 - q f t f расчет остатка X и y 6.6.f = f + 1f = f + 1переход к следующему фактору

ПЛС1 и ПЛС2 ПЛС1 моделирует только одну переменную y «за раз» ПЛС2 позволяет моделировать любую комбинацию переменных Y без их разделения – совместно он кажется более подходящим при калибровке нескольких свойств… однако, ПЛС1 дает по отдельной модели на каждое из интересующих свойств, возможно, с различным числом факторов не будет ли набор независимых моделей всегда лучшим решением? однозначного ответа нет… сравним методы на практике!

λ, нм ε, M -1 см C ε + E = D R.S.D. (E) = C = C 0 + E c R.S.D. (E c ) = 5% (C max ) nКомпонент[C n ], M 12-ацетофенантрен ацетиламинофенантрен ацетиламинофенантрен Строим ПЛС2-модель (Simdata)

Интерпретация модели служит для изучения внутренней структуры данных группы выбросы взаимовсвязи Сходство с РГК (PCR): X-счета и нагрузки (scores & loadings) Особенности: график t – u : метод обнаружения выбросов (outliers) графики нагрузок w – w : карта переменных cравнение двух X-нагрузок p – w : насколько Y повлияла на декомпозицию X график w – q Интерпретация ПЛС-моделей

Интерпретация моделей: ПЛС2 против РГК PLS2 PCR

Интерпретация моделей: ПЛС1 против ПЛС2

Интерпретация ПЛС-моделей: связь X и Y (Simdata)

Интерпретация ПЛС-модели: выбросы (Octane)

Проверка регрессионных моделей Проверка (validation) модели преследует две основные цели: Определение оптимального числа компонент Меньше факторов чем в РГК Минимум RMSEP Оценка предсказательной способности модели: График предсказанние относительно измерения (predicted vs measured) RMSEP

Проверка регрессионных моделей: simdata – ПЛС1

Сравнение моделей: Simdata МЛР MLR РГК PCR ПЛС1-Р PLS1-R ПЛС2-Р PLS2-R [C 1 ] [C 2 ] [C 3 ] Сравнение моделей калибровки трехкомпонентной смеси ПАУ (simdata) вывод: модели РГК, ПЛС1-Р, ПЛС2-Р примерно одинково хороши для калибровки этих данных (без осложнений) результаты МЛР значительно хуже, для [C 3 ] - неудовлетворительные

Сравнение методов калибровки МЛР (MLR) плохо пригоден для спектроскопических данных РГК (PCR) имеет недостатки, но хорошо работает при отсутствии осложнений ПЛС регрессия (PLS-R) является лучшим решением для большинства практических задач PLS1 или PLS2? Как выбрать метод? – пробовать! Как сравнивать разные модели? RMSEP

Линейная регрессия и нелинейность X: 100x351 A = εC – 0.1C 2 … Y: 100x1 r=0.999

Предсказание: диагностика соответствия новых образцов не все проблемы заканчиваются с построением калибровочной модели! возможность выявления образцов, не соответствующих данной регрессионной модели является одним из преимуществ проекционного подхода Deviation - эмпирический параметр, характеризующий меру соответствия нового образца калибровочной модели рассмотрим наш пример…

Диагностика предсказания (Simdata)

Диагностика предсказания: ПЛС1 - Simdata [C 1 ] = 0 – 1 M[C 2 ] = 0 – 0.5 M [C 3 ] = 0 – 0.05 M #[C 1 ][C 2 ][C 3 ] U U U U U U60.5

Правила построения «хорошей» калибровки правильно приготовить (собрать) образцы визуально изучить данные, если необходимо, применить предварительную обработку данных (pre-processing) если необходимо применить шкалирование/ взвешивание (scaling/weighting) интерпретировать модель, изучить структуру данных, выявить и удалить возможные выбросы тщательно оценить размерность модели, диагностировать модель диагностировать предсказание

План семинара Пример 1. Концентрационная калибровка трехкомпонентной смеси ПАУ по спектрам в УФ- видимой области (искусственные данные). общие навыки калибровки, интерпретации и диагностики модели, предсказания на «идеальных» данных Пример 2. Определение октанового числа топлива по спектрам ближнего ИК. калибровка на реальных данных, обнаружение и удаление выбросов Пример 3. Качество пшеницы (факультативно). самостоятельное построение калибровки, MSC, выбор переменных

Рекомендуемая литература Richard Kramer Chemometric Tchniques for Quantitative Analysis * Kim H. Esbensen Multivariate Data Analysis - in Practice ** Kenneth R. Beebee et al. Chemometrics: a Practical Guide ** Harald Martens, Tormod Naes Multivariate Calibration ** Richard G. Brereton Chemometrics: Data Analysis for the Laboratory and Chemical Plant *** Edmund R. Malinowski Factor Analysis in Chemistry ****

Пример 1: Калибровка смеси ПАУ Файл Simdata Цель: выработка навыков калибровки с программой Unscrambler изучить наборы данные: обучающий, тестовый, «unknown» - в таблице, как серии спектров построить калибровки: РГК, ПЛС2 - сравнить модели построить ПЛС1 для каждого из 3-х компонентов, определить размерность моделей изучить графики scores, loadings, T-U, predicted vs measured, RMSEP, Variance для [С1] - [С3] с разным количеством факторов предсказать «неизвестные» образцы

Пример 2: Определение октанового числа бензина стр. 139, файл Octane Цель: работа с реальными данными, диагностика и устранение выбросов преимущественно по книге: построить калибровку ПЛС1, диагностировать определить выбросы, удалить, обновить модель проверить модель различными способами, включая тестовый набор построить РГК, сравнить модели предсказать «неизвестные» образцы

Пример 3: Качество пшеницы стр. 150, файл Wheat Цель: самостоятельное построение калибровочной модели построение моделей ПЛС1/2, сравнение моделей определение и удаление выбросов применение MSC попробовать удаление переменных для улучшения модели