CC BY
487
Известия Тульского государственного университета Естественные науки. 2012. Вып. 1. С. 153-162 = ИНФОРМАТИКА
УДК 519.252
Многоклассовая логистическая регрессия для прогноза вероятности наступления инфаркта *
А. П. Мотренко, В. В. Стрижов
Аннотация. В работе описан алгоритм классификации четырех групп пациентов: перенесших инфаркт, имеющих
предрасположенность к инфаркту и здоровых пациентов двух типов. Признаками для определения состояния пациента служат измерения концентрации белков в крови. Работа посвящена прогнозу вероятности принадлежности пациента к одному из нескольких неупорядоченных классов. Решается задача оценки параметров функции регрессии и выбора признаков при многоклассовой классификации. Классификация выполняется по всевозможным парам групп.
Ключевые слова: логистическая регрессия, многоклассовая
классификация, выбор признаков, прогноз предрасположенности к инфаркту.
Введение
Заболевания сердечно-сосудистой системы могут протекать, не проявляясь клинически. Тем не менее, обнаружение нарушений, связанных с работой сердца, по косвенным признакам вполне возможно [1, 2]. В данной работе в качестве признаков (биомаркеров) используются концентрации белков и их соединений, абсорбированные на поверхности кровяных телец. Разделение пациентов на группы по состоянию здоровья приводит к задаче многоклассового прогнозирования. Эта задача сведена к задаче двуклассовой классификации; используется подход «каждая группа против каждой». В этом случае рассматриваются все возможные пары групп пациентов и решается задача вида «к какой из двух данных групп пациент принадлежит с большей вероятностью?». Данный подход принят в связи с относительно небольшим объемом выборки, на которой проводился вычислительный эксперимент.
* Работа выполнена при финансовой поддержке РФФИ (проект № 10-07-00422).
Для каждой пары групп решается задача логистической регрессии [3], в основе которой лежит предположение о биномиальном распределении независимой переменной, и оцениваются параметры функции регрессии [4,
5].
Предполагается, что число измеряемых признаков избыточно; требуется отыскать оптимальный набор признаков, эффективно разделяющий классы. Признаки в логистической регрессии как правило выбираются с помощью шаговой регрессии [6, 7]. В данной работе используется полный перебор, т.к. он дает экспертам гарантию, что рассмотрены все возможные сочетания признаков при выборе модели. При этом экспертами вводились ограничения на сложность модели. Задача выбора признаков поставлена с использованием площади под ИОС-кривой [8] в качестве внешней функции ошибки.
Задача классификации сопряжена с оценкой минимального объема выборки, достаточного для проведения классификации. Для этого используются метод доверительных интервалов [9], метод скользящего контроля [10], сравнение предполагаемых распределений на различных подвыборках [11].
При проведении вычислительного эксперимента и прогноза вероятности наступления инфаркта были использованы данные [12], предложенные специалистами парижской лаборатории анализа крови «Иммуноклин».
1. Задача классификации и оценка параметров
Дана выборка Б = {(хг ,уг)} ,г = 1,...,т, состоящая из т объектов (пациентов), каждый из которых описывается п признаками (биомаркерами) Хг £ М” и принадлежит одному из двух классов Уг £ {0,1}. Рассмотрим задачу логистической регрессии. Предполагается, что вектор ответов у = [у1,... ,ут]Т — бернуллиевский случайный вектор с независимыми компонентами уг ~ В(рг, 1 — рг) и плотностью
т
Р(уИ = П РУ (1 — Рг)-. (1)
г=1
Определим функцию ошибки следующим образом:
т
Е(ш) = — 1пр(у\ш) = —^2 Уг 1пРг + (1 — У г) 1п (1 — Рг). (2)
г=1
Другими словами, функция ошибки есть логарифм плотности, или функции правдоподобия, со знаком минус. Требуется оценить вектор параметров , доставляющий минимум функции ошибки:
= агг шш Е(^).
(3)
Вероятность принадлежности объекта к одному из двух классов определим
как 1
р = -(—~) = а(хТ w) = (4)
1 + ехр(-хТ w)
Для оценки параметров, воспользовавшись тождеством
йа(в) йд
вычислим градиент функции Е^):
= а(1 — а),
m
VE (w) = — ^ (y»(1 — ai) — (1 — yi)ai) x» = ^(ai — y» )x = XT (a — y), i=1 i=1
lT
где вектор а = [о1,...,от]Т и матрица X = [хТх^] состоит из
векторов-описаний объектов.
Оценка параметров осуществляется по схеме Ньютона-Рафсона. Введем обозначение X — диагональная матрица с элементами Ец = оі(1 — оі), г = = 1,..., т. В качестве начального приближения w = [аді,..., адга]Т вектора возьмем
m
wj = 2_j Уі(1 — Уi), j = 1,---,п.
i=1
Оценка параметров w^+1 логистической регрессии (4) на k + 1-м шаге итеративного приближения имеет вид
wfc+i = wfc — (XT £X)-1XT (a — y) = (XT £X)-1XT S(Xwfc — £-1(a — y)).
(5)
Процедура оценки параметров повторяется, пока норма разности ||wfc+1 — w^у2 не станет достаточно мала.
Алгоритм классификации имеет вид:
a(x) = sign(a(x, w) — a0), (6)
где a0 — задаваемое в (8) пороговое значение функции регрессии (4).
Вычисления качества прогноза. Для контроля за качеством прогноза множество индексов объектов I разбивается случайным образом на два подмножества, I = LUT, обучающее и тестовое. Параметры w оцениваются на подвыборке Dl, а качество прогноза вычисляется на подвыборке Dt. В данной работе для оценки качества прогноза и для выбора признаков используется один из двух функционалов: площадь AUC под кривой ROC и функционал
Q = (1 — TPR)2 + FPR2,
где
1m
m
i=1
TPR = — V[a(xi) = 1][yi = 1]
m ^
есть доля объектов выборки, правильно классифицированных в пользу данного класса, и
1 m
FPR = — ^[a(xi) = 1][yi = 0]
— . л
i=1
есть доля ошибочно классифицированных в пользу данного класса объектов выборки. Здесь используется обозначение индикаторной функции:
[У = 1] = Г’ У = 1; (7)
W J \0, y = 1. ^
Таким образом, алгоритм тем лучше разделяет классы, чем меньше значение функционала Q или чем больше значение площади AUC под кривой ROC. Отложив на графике для каждого значения связанной переменной £ G [0,1] по оси абсцисс значения FPR(£), а по оси ординат — TPR(£) получим кривую ROC, каждая точка которой соответствует некоторому значению оо Отыскание параметра оо алгоритма классификации. В алгоритме
(6) используется то значение оо, которое соответствует наибольшему расстоянию от отрезка [(0,0),(1,1)], означающего отказ от принятия решения о классификации, до кривой ROC:
сто = arg max ||(TpR(£)’FPR(£)) - . (8)
Последнее выражение включает вычисление значения функционала качества, и как следствие, вычисление выражения (6) и итеративную оценку параметров (5).
2. Выбор признаков в задаче классификации
Введем обозначения A — некоторое подмножество индексов признаков, J = {1’...’П} и A — оптимальный набор индексов. Обозначим Хд множество столбцов-признаков матрицы X, заданное набором A и wa — соответствующие им параметры. Рассмотрим задачу выбора признаков как задачу максимизации:
A = argmaxAUC(A) при условии |A| = const. (9)
В задаче использована площадь под кривой AUC(A) = AUC(X^Wa,<г0,y), значение которой вычислено для набора индексов признаков A, а параметры Wд и о0 получены в результате решения задач (3) и (8).
Набор признаков отыскивается путем полного перебора. Такой подход возможен благодаря сравнительно небольшому количеству признаков в данной задаче и диктуется требованиями экспертов. Запишем выражение для функции регрессии o(xTw) в виде
xT w = aixiiwi + a2Xi2W2 + ... + anXinWn,
где а^ € {0,1} — структурный параметр. Таким образом, алгоритм нахождения признаков сводится к перебору значений элементов вектора [а\,..., ап] структурных параметров
а1 а2 ... ап
~ 0 ... 0“
0 1 ... 0 .
1111
В данном случае не оговаривается необходимость разбиения выборки О на обучающую и тестовую подвыборки, так как эксперты зафиксировали максимальное число признаков при решении задачи: |Д| не должна превышать четырех. Наборы признаков, полученные в результате решения задачи (9), будем называть оптимальным для данной пары классов, а сами признаки — наиболее информативными (рис. 1 и 2).
4-markers sets
Рис. 1. Вариационный ряд значений функционала для всех наборов признаков фиксированной мощности
3. Прогноз при многоклассовой классификации
Пациенты из исследуемой выборки разделены по состоянию здоровья на четыре группы:
Ai — пациенты, уже перенесшие инфаркт,
A3 — пациенты, имеющие предрасположенность к инфаркту,
Bi, B2 — здоровые пациенты двух типов.
При появлении в выборке нового объекта Xm+i, состояние которого необходимо спрогнозировать, выполняется следующая процедура. Для каждой пары групп из шести возможных пар выполняется
2 4 6 8 10 12 14 16 18 20
Рис. 2. Число вхождений каждого из двадцати маркеров в набор «К
лучших»
классификация (6). Используется оптимальный для соответствующей пары набор признаков (9). В каждом случае алгоритм а(хт+і) возвращает решение о принадлежности объекта к одному из двух рассматриваемых классов и вероятность принадлежности рт+і = ^хт+і^. По этим результатам составляется таблица следующего вида:
Аі А3 Ві В2
¿1 — 0 0 1
А3 1 — 1 1
Ві 1 0 — 0
В2 0 0 1 —
Например, третья строка содержит следующие результаты классификации: объект хт+1 скорее принадлежит классу Б\, чем классу Л\, и скорее принадлежит классам Аз и Б2, чем классу Б1. Присвоим классам А1,Аз,Б1,Б2 номера 1, 2, 3, 4 и введем нижние индексы а^(х) € € {0,1}, 1,к € {1,..., 4} для пары классов (I, к). Тогда для рассмотренного примера
а2з(хт+1) = 0.
При прогнозировании объект относится к тому классу, для которого сумма элементов таблицы по строке наибольшая:
4
с1а88(хт+1) = ащ шах V аш (хт+1). (11)
Если эта сумма для двух классов совпала, результатом будет классификация (6), полученная для этих двух классов.
4. Результаты вычислительного эксперимента
Вычислительный эксперимент проводился на данных лаборатории анализа крови «Иммуноклин» [12]. Данные содержат измерения концентрации 20-ти белков и их соединений на поверхности кровяных телец 98-ми пациентов четырех классов; в каждом классе примерно равное число пациентов. В табл. 1 приведен список исследуемых биомаркеров с их порядковыми номерами.
Таблица 1
Список исследуемых биомаркеров
1 2 3 4 5 6 7 8 9 10
к Ь к/м Ь/М к/ы к/о Ь/О к/р Ь/р к/д
11 12 13 14 15 16 17 18 19 20
к/и Ь/И Ь/И/ЯЛ Ь/Т/ЯЛ Ь/Т/ЯО и/У и/ш и/х и/У и/г
Построение процедуры прогнозирования выполнялась следующим образом.
(1) Для множества объектов обучающей выборки, принадлежащих каждой паре классов, рассматривались все наборы {Л} признаков мощностью 4.
(2) Оценивались параметры ,<го алгоритма классификации, соответствующие набору Л.
(3) Вычислялось значение функции качества алгоритма классификации ЛИС на выборке, состоящей из пары классов.
(4) Для каждой пары классов выбраны К наборов признаков, доставляющих большее значение функции качества ЛИС.
Таблица 2
Результаты выбора признаков
Классы Ш2 Л АиС(Д)
со 1 31 14 [2, 11, 19, 20] [2, 13, 19, 20] 0.953 0.953
А - в 55 14 [3, 13, 18, 19] [12, 13, 15, 19] 0.829 0.829
А\ — в 55 14 [5, 15, 17, 19] [б, 12, 15, 19] 0.901 0.901
Аз — Б\ 58 17 [5, 6, 11, 17] [2, 7, 9, 13] 0.814 0.829
Аз — В2 43 17 [2, 3, 5, 9] [2, 3, 9, 19] 0.954 0.957
В1 — В2 67 41 [1, 2, 3, 9] [2, 3, 9, 11] 0.821 0.823
В табл. 2 для каждой пары классов указаны наборы маркеров, доставивших наибольшие значения максимизируемому критерию ЛИС, и сами значения этого критерия. Приведены два набора для каждой пары классов. Для исследования были выбраны пять лучших наборов. Число наборов К задано экспертами, исходя из рис. 1, на котором изображен вариационный ряд (в порядке возрастания) значений ЛИС, полученных при классификации на различных наборах {Л}. Предлагается выбирать такое число К, при котором рост графика меняется еще достаточно сильно (справа налево). В данном случае выбрано значение К = 5 (табл. 3).
Таблица 3
Число вхождений признаков в К оптимальных наборов для каждой пары
классов.
Классы к Ь к/м к/м к/О Ь/О к/р Ь/р к/и
Аг — Аз 0 5 0 0 0 0 0 0 1
Аг — Вг 0 0 1 0 0 0 1 0 0
Аг — В2 0 0 1 1 1 0 0 1 0
Аз — Вг 0 3 1 3 2 2 0 3 1
Аз — В2 0 5 4 1 0 0 0 5 0
Вг — В2 2 5 3 0 0 0 0 3 1
Классы Ь/И Ь/И/ЯЛ Ь/Т/ЯЛ Ь/Т/ЯО и/У и/ш и/х и/У и/г
Аг — Аз 0 1 1 0 1 1 0 5 5
Аг — Вг 2 4 0 2 0 0 4 5 1
Аг — В2 4 0 0 5 0 2 0 5 0
Аз — Вг 0 2 0 0 0 1 1 0 1
Аз — В2 0 1 0 0 0 1 1 2 0
Вг — В2 0 1 0 0 0 0 3 2 0
Согласно полученным результатам, итоговый алгоритм прогнозирования имеет следующий вид.
(1) Для нового объекта выполняется классификация на всех К наборах признаков каждой пары классов.
(2) Строится таблица (10), включающая каждую пару классов К раз.
(3) Решается задача прогнозирования (11) с учетом числа наборов К. Одной из важных практических задач, решаемых в рамках проводимых
исследований, является задача снижения стоимости клинического обследования одного пациента, решаемая путем уменьшения числа измеряемых биомаркеров. Предложено измерять только наиболее информативные биомаркеры, выбранные следующим образом. Для ]-той
к
пары классов найдено множество оптимальных наборов Бз = и Л^,
i=1
] = 1,..., 6. Объединив признаки из всех наборов из колонки «Л» табл. 2,
6
получим множество наиболее информативных признаков У Бз. Для
3=1
каждого признака подсчитано количество его вхождений в это множество.
Гистограмма на рис. 2 показывает, насколько часто каждый признак входит в K лучших наборов, полученных для каждой пары классов. Таблица 3 показывает число вхождений биомаркера в набор наиболее информативных признаков каждой пары групп пациентов.
Заключение
В работе описан алгоритм прогнозирования вероятности наступления инфаркта пациентов при многоклассовом прогнозировании; описан способ оценки параметров и выбора наиболее информативных признаков. Выборка пациентов разбита на четыре группы относительно наличия нарушений в работе сердечно-сосудистой системы. При этом задача многоклассовой классификации сводилась к двуклассовой путем
рассмотрения всевозможных пар групп. Для каждой из таких пар получен оптимальный набор признаков. Получена оценка качества прогнозирования в парах групп на исследуемой выборке.
Список литературы
1. Azuaje F, Devaux Y, Wagner D. Computational biology for cardiovascular biomarker discovery // Brief Bioinform. 2009. V.10. №4. P.367-377.
2. Transcriptomic biomarkers for individual risk assessment in new-onset heart failure / Heidecker [et al.] // Circulation. 2008. V.118. №3. P.238-246.
3. Hosmer D., Lemeshow S. Applied logistic regression. N. Y.: Wiley, 2000. 375 p.
4. Bishop C.M. Pattern recognition and machine learning. Springer, 2006. 738 p.
5. MacKay D.J.C. Information theory, inference, and learning algorithms. Cambridge University Press, 2003. 628 p.
6. Friedman J, Hastie, Tibshirani R. Additve logistic regression: a statistical way of boosting // The Annals of Statistics. 2000. V.28. №2. P.337-407.
7. Madigan D, Rideway G. Discussion of least square regression / В сб. Efron B. [et al.]. Least Angle Regression // The Annals of Statistics. 2004. V.32. №2. P.465-469.
8. Fawcet T. ROC graphs: notes and practical considerations for researchers // HP Laboratories. 2004. 38 p.
9. Реброва О.Ю. Статистический анализ медицинских данных. Применение прикладного пакета Statistica. М.: Медиасфера, 2002. 312 с.
10. Bos S. How to partition examples between cross-validation set and fraining set? / Saitama. Japan: Laboratory for information representation RIKEN. 1995. 4 p.
11. Perez-Cruz F. Kullback-Leibler divergence estimation of continuous distributions // IEEE International Symposium on Information Theory, 2008.
12. Standart flow cytometry analysis of nondental patients. Paris: ImmunoClin laboratory. 2007. 1 p.
Стрижов Вадим Викторович ([email protected], http://strijov.com), к.ф.-м.н., н.с., Вычислительный центр Российской Академии Наук, Москва.
Мотренко Анастасия Петровна ([email protected]), студент, Московский физико-технический институт.
Multiclass logistic regression for cardio-vascular disease
forecasting
A. P. Motrenko, V. V. Strijov
Abstract. The paper describes an algorithm to classify four groups of patients: a cardio-vascular disease group, a cardio-risk group and two types of healthy groups. The blood-cells protein measurements are the description features for an investigated patient. The paper develops an algorithm to forecast a patient’s cardio-vascular disease case as one of four unordered classes. The problem is to estimate the regression parameters and select the most informative features for multi-class classification. During the forecasting all pairs of the classes are considered.
Keywords: logistic regression, multiclass classification, feature selection, cardio-vascular disease forecasting.
Strijov Vadim ([email protected], http://strijov.com), candidate of physical and mathematical sciences, researcher, Computing Center of the Russian Academy of Sciences, Moscow.
Motrenko Anastasiya ([email protected]), student, Moscow Institute of Physics and Technology.
Поступила 01.02.2012
