ПЛС-регрессия и ковариация

Агнар Хосдкульдсон, Датский технический университет
Homepage of Chemometrics, Editorial, July 2004
Agnar Höskuldsson
IPL, DTU, Blgd 358, 2800 Kgs.
Lyngby, Denmark

ПЛС-регрессия характеризуется оптимизационным по отношению к исходным данным процессом, при котором ведутся поиски вектора счетов в пространстве X и вектора счетов в пространстве Y, обеспечивающих максимальную ковариацию. При анализе данных, часто возникает необходимость в наглядном представлении результатов такого процесса. В данной статье рассказывается, каким образом это сделать лучше всего. Кроме этого, мы так же проиллюстрируем, как при этом выделить наиболее полезную, с точки зрения задачи моделирования, информацию.

Задача оптимизации

ПЛС-регрессию можно рассматривать как пошаговую процедуру, на каждом этапе которой вычисляется вектор счетов t по следующей формуле:

t = Xw = x1w1+…+xkwk

где, X = Xa — приведенная матрица плана (матрица NxK, X = (x1,…,xk)), соответствующая a-тому шагу процедуры. Неизвестным параметром является вектор взвешенных нагрузок w = (w1,…,wk). И, таким же образом вычисляется вектор счетов u,

u = Yq = y1q1+…+ymqm

где Y = Ya — приведенная матрица откликов (матрица NxM, Y = (y1,…ym)), а неизвестным параметром являются вектор взвешенных нагрузок q = (q1,…,qm). Задача оптимизации в ПЛС-регрессии выглядит как:

maximize(tTu) = maximize(t1u1+…+tnun), при условии |w| = |q| = 1


Рис.1. Иллюстрация решения задачи оптимизации

Используя множитель Лагранжа можно показать, что w задается собственным вектором, принадлежащим наибольшему собственному значению в системе линейных уравнений XTYYTXw = λw, а u вычисляется из уравнения u = YYTXw.

Такой взгляд на процедуру оптимизации подразумевает, что при моделировании матрица X и матрица Y имеют одинаковые векторы взвешенных нагрузок. Задачу оптимизации можно поставить и по-другому, только относительно вектора w:

maximize|q|2 = maximize|YTXw|2, при условии |w| = 1

Отсюда видно, что мы должны найти такие взвешенные нагрузки w, которые доставляют максимум вектору Y-нагрузок — q. Такая интерпретация особенно важна, если результаты задачи моделирования будут в дальнейшем использоваться для построения других более детализированных моделей.

Максимальная ковариация, достигаемая на данном этапе, равна наибольшему собственному значению λ. Это можно записать как λ = tTYYTt = |YTt|2. Полезно после каждой итерации проверять насколько сильно уменьшился этот показатель.

Когда вектор взвешенных нагрузок w найден, вычисляются векторы t = Xw, q = YTt и u = Yq, после чего можно проанализировать полученные результаты детально. Рассчитывается коэффициент регрессии Y на t как dq, где d = 1/(tTt). Перед следующим шагом матрицы X и Y пересчитываются следующим образом:

X <− X − d t pT
Y <− Y − d t qT

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

Рис.2. Слева — график значений лямбда и соответствующего им верхнего доверительного предела, справа — квадратный корень из лямбда и соответствующий ему верхний доверительный предел

Анализ ковариационной матрицы данных

Процедура, описанная выше, рассматривает в качестве основной величины квадрат ковариационной матрицы C = XTYYTX. Максимальное собственное значение, которое можно получить, определяется как: λ = wTCw = tTYYTt. Обычно весьма полезным является отслеживание изменения величины λ на каждом этапе, как показано на рисунке 2. В качестве данных взяты результаты БИК-спектроскопии, которые были предварительно автошкалированы, т.е. центрированы и взвешены, так что дисперсия по каждой переменной стала равняться единице. Всего имеется пять y-переменных. Рисунок показывает, что λ уменьшается достаточно быстро. Обычно когда λ достигает нуля, это верный признак того, что была достигнута подходящая размерность. Таким образом, из левой части рисунка мы можем заключить, что размерность нашей модели должна быть равна 4 или 5. Кроме этого на рисунке так же изображен верхний доверительный интервал для λ. На правом графике изображен корень квадратный из λ и соответствующий ему верхний доверительный интервал. На нем отчетливо видно, что графики пересекаются при размерности равной 6, что является признаком того, что для модели нужно использовать размерность 5. Формула для определения границы верхнего доверительного интервала представлена в приложении. Если мы проведем более подробный анализ, то обнаружим, что на самом деле наиболее подходящая размерность равна четырем, а число переменных — 120.

Кроме этого, полезным бывает анализ диагональных элементов матрицы C,

Ckk = (xkTy1)2 +…+(xkTym)2,  k = 1,…,K
Рис.3. Графики квадрата ковариации для первых шести размерностей

Графики значений этой величины для первых шести размерностей изображены на рисунке 3, из которого можно сделать два основных вывода. Первый заключается в том, что значение ковариации становится пренебрежимо малым для размерностей 5 и 6, из чего следует, что для построения модели достаточно четырех компонент. Второй вывод состоит в том, что большая часть переменных не коррелирует со значениями y, начиная с размерности равной двум. В случае, подобном этому, когда большая часть переменных не вносит никакого вклада в модель необходимо исключить их из анализа. Однако мы не будем рассматривать здесь это более подробно.

Заключение

ПЛС-регрессия основывается на решении задачи оптимизации, отдельной для каждой размерности. При этом весьма полезно отслеживать информацию о том, насколько хорошо был выполнен каждый этап оптимизационного процесса. В данной статье мы показали, как получить такую информацию. Графики на рисунках 2 и 3 весьма полезны для демонстрации "энергии" в данных, которая приводит к искомому результату.

Приложение

Здесь мы коротко покажем, как вычислить верхний 95% доверительный интервал для квадрата ковариации λ. Более подробно об этом можно прочитать в [1]. Предполагается, что теоретически ковариация между X и Y представляет собой нули, Σxy = 0. В [1] показано, что мы можем приближенно записать:

E(λ)  tr(tTt) tr(YTY)/N, и Var(λ)  2 (tTt)2 tr([YTY]2)/N2.

Через tr() обозначается след матрицы, λ = (tTy1)2+…+(tTym)2. Можно показать, что каждый элемент (tTyi) имеет почти нормальное распределение (это следует из центральной предельной теоремы). Используя эти выражения, мы можем получить приближенное распределение для λ. Однако на практике можно считать, что λ имеет почти нормальное распределение со средним значением и вариацией, определенными выше. При нормальном распределении со средним μ и стандартным отклонением σ, 95% верхний доверительный интервал задается выражением μ+σ*1.645. Применяя это к λ, мы получим следующую формулу:

E(λ)  tr(tTt) tr(YTY)/N+1.645 {2(tTt)2 tr([YTY]2)/N2}1/2.

Это выражение и использовалось на графиках, представленных на рисунке 2. Оно вычислялось перед каждой итерацией, т.е. и верхний доверительный интервал вычислялись одновременно. Такой же доверительный интервал можно изобразить и для кривых на рисунке 3.

Литература

1. Prediction Methods in Science and Technology. Vol. 1. Basic Theory. Thor Publishing, Copenhagen, ISBN 87-985941-0-9.