Численно устойчивый вероятностный классификатор логистической регрессии Текст научной статьи по специальности «Математика»

АПВПМ-2019

ЧИСЛЕННО УСТОЙЧИВЫЙ ВЕРОЯТНОСТНЫЙ КЛАССИФИКАТОР ЛОГИСТИЧЕСКОЙ РЕГРЕССИИ

В, Л, Лукинов

Институт вычислительной математики и математической геофизики СО РАН, 630090, Новосибирск Новосибирский государственный университет, 630090, Новосибирск Сибирский государственный университет телекоммуникаций и информатики, 630102, Новосибирск

УДК 519.254

DOI: 10.24411/9999-016А-2019-10047

При обработке реальных данных методом логистической регрессии требуется определить и устранить негативное влияние на качество оценок коэффициентов модели, связанное с неполнотой данных, существованием "выбросов" в данных и особых ковариантов, коллинеарностью, критериями отбора ковариантов в многофакторную модель. В данной работе предложен и реализован алгоритм, позволяющий контролировать и устранить перечисленные выше негативные факторы. Другой задачей является разработка параллельных алгоритмов расчета коэффициентов логистической регрессии и характеристик качества. При этом необходимо решить оптимизационную задачу, реализующую критерий максимального правдоподобия. Известный подход, основанный на итерационном методе Ныотона^Рафсона, является численно неустойчивым и может привести к неправильному нахождению коэффициентов регрессии [9]. В данной работе предлагается новый параллельный численно-устойчивый итерационный алгоритм решения задачи минимизации на основе случайного поиска. Проведено сравнение с базовым алгоритмом glm из языка статистической обработки данных R [8]. Ключевые слова: задача классификации, бинарная логистическая регрессия.

Введение

Одним из наиболее распространенных и широко используемых в статистическом анализе больших данных является разработанный в прошлом веке метод логистической регрессии [1]. В биологии и медицине логистическая регрессия применяется в самых разных областях: выявлении и исследовании степени влияния предикторов заболеваний и послеоперационных осложнений, фундаментальном для ретроспективных групповых сравнительных исследований методе Propensity Score Matching, методах автоматического и полуавтоматического распознавания медицинских изображений, методах статистического кластерного анализа [2-5].

1 Задача бинарной логистической регрессии

Модель логистической регрессии строится на обучающей выборке, определяемой парами вида где

/

Xj =

V

является набором значений г-го объекта (значения ковариант), бинарные переменные yi определяют принадлежность одному из двух классов, например, yi == — 1 — для первого класса и уj == 1 — для второго класса, т — 1 — количество ковариант (признаков) у каждого объекта, п — число наблюдений.

Работа выполнена при частичной финансовой поддержке Российского фонда фундаментальных исследований (код проекта 17-01-00698 А) и Института вычислительной математики и математической геофизики СО РАН (гос. задание 0315-2016-0002).

ISBN 978-5-901548-42-4

(Xj,^), г = 1, 2 ,...,п,

(1)

Хц \ Xi2

1

Xi2 \ Xim J

( xf \ ( 1 X\2

X

x2

V xn /

1 ^22 V 1 X„2

^2m

(2)

x

x

nm

Конечной целью является построение апостериорной функции принадлежности произвольного объекта х классу у == —1, называемой функцией логистической регрессии, имеющей вид [1]

р(х)

1 + exp(—bтх) '

и принимающей значения в интервале (0; 1). Принадлежность объекта х классу у == 1 определяется через дополнение вероятности как 1 — р(х). Объекты с полученным значением р(х) > 0.5 относят к первому классу.

( х)

ции на априорных данных можно получить методами ROC анализа [6].

Известно [7], что поиск функции логистической регресии можно представить в виде задачи минимизации функционала

Q(b) = ¿In (l + e-y'bTx>) ^ mrnm, (3)

• =1

где bT = (bu b2, ..., bm) — искомый вектор коэффициентов регрессии, представляющий границу классификации объектов в виде гиперплоскости

т

Ь\ + biXi = 0

i=2

2 Итерационный метод поиска минимума

По причине отсутствия аналитического решения задачи (3), для численного поиска коэффициентов регрес-b

том числе и метод соответствия обобщенных линейных моделей glm базового пакета stats в языке R [8]) использует метод Ныотона-Рафсона. Далее приводится описание основных шагов этого алгоритма. В качестве нулевого приближения берется "наивное" решение задачи классификации методом многомерной линейной регрессии:

b(o) = (XTX)-1 (XTY) . (4)

На к-ом итерационном шаге происходит уточнение вектора коэффициентов b(k):

b(k) = b(k-1) - kk(Q"(b(k-1))jQ'(b(k-1)), (5)

где Q'(b(k-1)) — градиент функционала Q в точке b(k-1); Q"(b(k-1)) — матрица вторых производных функционала Q в точке b(k-1) ; hk — шаг итерации, который может, в тривиальном случае, быть единичным, "оптимальным" фиксированным или изменяться на каждой итерации [7].

Для исследования свойств алгоритма (5) выразим градиент функционала Q'(b) и матрицу вторых производных Q"(b). Обозначим значение так называемой сигмоидной функции как а = ai = ffi(y¿bTXj) = 1/(1 + exp(—уjbTXj)). Заметим, что для производной сигмоидной функции справедливо следующее соотношение а'(b) = ai(b)(1 — ai(b)). Выразим первые частные производные функционала Q(b):

) = — (1 — ai)ViXk, к = 1,. .. ,т (6)

^ г=1

()

d2Q(b) д ^ .А .

дъ дъ = — дг- 2^(1 — ai)yiXk = }_^ai(1 —ai)xjxk j = 1,...,m, к = 1,...,m (7)

J k i i=1 i=1

Положим:

W„x„ = diag j уЯЛ—О-) I — диагональная матрица весов наблюдений; Я-пхт = W„x„X„xm — матрица взвешенных ковариант;

X = (1 — а^/а^, у = (гу)г=\..,п — взвешенные отклики наблюдений.

Тогда произведение обратной матрицы вторых производных на вектор градиента из (5) в новых обозначениях можно записать в виде:

(я "(Ь))-1 а '(Ь) = — (х^ 2х) -1ХТУ/ у = — (хТ х)-1 ХТ у (8)

3 Причины неустойчивости и медленной сходимости

Существует несколько причин возможной численной неустойчивости и медленной сходимости метода Ныо-тона-Рафсона (4), (5). В случае медленной сходимости, итерации в программе могут прерваться из-за принудительного ограничения количества итераций и не достичь минимума функционала с заданной точностью. Рассмотрим некоторые причины.

Коррелируемость ковариант. Наихудшим случаем является линейная зависимость признаков, тогда мат-

Т

рица Ц "(Ь) = XX является вырожденной и не может быть обратимой. В статистических расчетах па реальных данных такая ситуация редко возникает при формальных расчетах без исследования свойств наблюдений. Намного чаще возможен случай коррелируемых или так называемые мультиколлинеарных ковариант, при котором матрица ф''(Ь) обратима, но близка к вырожденной. Столбцы такой матрицы приближенно линейно зависимы в пределах малой точности, а сама матрица является плохо обусловленной. Причиной плохой обусловленности является наличие собственных значений близких к нулю. Действительно, число обусловленности матрицы есть

№ ») = ^,

^тг п

где ^тах и \тт соответственно максимальное и минимальное собственное значение матрицы. Известно [10], что обращение плохо обусловленной матрицы численно неустойчиво и умножение обратной матрицы на вектор может увеличивать относительную погрешность в ^ раз.

Из (4) следует, что мультиколлинеарность может являться причиной плохого начального приближения, что в свою очередь влияет на число итераций.

Кроме того, плохая обусловленность матрицы ф''(Ь) повышает дисперсии коэффициентов логистической регрессии, что приводит к неверной оценки вероятности нулевой гипотезы с помощью Е-статнстпк и статистик Вальда.

Резюмируя, можно сказать, что коррелируемость ковариант является причиной:

— неправильной оценки параметров Ь в следствии численной неустойчивости;

— неверного вычисления уровня значимости р для коэффициентов регрессии;

— увеличения вычислительных итераций вплоть до принудительной остановки подпрограммы д1т.

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

()

пых минимумов [11]. В этом случае возможен случай, когда, попав в точку локального минимума, алгоритм завершает свою работу, предоставляя значения коэффициентов Ь для этой точки, вместо глобального минимума.

4 Исходный метод построения однофакторных и многофакторных моделей логистической регрессии

Во многих клинических исследованиях логистическая регрессия применяется для выявления отдельных предикторов послеоперационных осложнений (однофакторные модели) и мультипликативных предикторов (многофакторные модели) [3,12,13]. Краткая схема метода:

1. построение для каждой коварианты однофакторной модели логистической регресси;

2. отбор в многофакторную модель ковариант с достигнутым уровнем значимости коэффициентов в од-нофакторных моделях меньше фиксированного порога значений;

3. при необходимости — устранение части ковариант из отобранных для устранения мультиколлипеарпо-сти;

4. построение многофакторной модели логистической регрессии;

5. итеративное построение оптимальной по информационному критерию Акаике модели многофакторной логистической регрессии одним из алгоритмов прямого, обратного или смешанного шага;

6. верификация построенных моделей на использованных для обучения данных.

В таблице 2 представлен пример оформления результатов расчетов. Включение возраста в многофакторную модель привело к значительному увеличению достигнутого уровня значимости, что является признаком неустойчивого счета в следствии коррелируемости возраста и стеноза.

Таблица 1: Логистические регрессии для предикторов осложнений

Однофакторные модели Полная многоф. модель Оптим. многоф. модель

Коварианты ОШ [95 % ДИ] Р ОШ [95 % ДИ] Р ОШ [95% ДИ] Р

Возраст 0.96 [0.897; 1.02] 0.162 0.991 [0.92; 1.07] 0.812

хнмк 1.5 [0.959; 2.48] 0.074 1.51 [0.902; 2.51] 0.117 1.62 [0.977; 2.68] 0.062

ПК 1 [1; 1-02] 0.003 1.01 [1; 1.02] 0.008 1.01 [1; 1.02] 0.003

Стеноз 0.94 [0.889; 0.998] 0.042 0.969[0.915; 1.02] 0.266

5 Численный метод

Представим модифицированный устойчивый параллельный метод, выполняющий расчеты для исходного метода:

варианты, при этом сохранаяем полученные значения свободного коэффициента Ъ\, и коэффициента коварианты Ь2; с полученной оценкой стандартного отклонения коэффициента ¿2;

2. отбираем для включения в полную модель коварианты с достигнутым уровнем значимости меньше порогового уровня;

3. при необходимости — удаляем часть ковариант из отобранных для устранения мультиколлинеарности;

4. устраняем коварианты с процентом пропущенных значений больше заданного порогового значения;

5. формируем вектор начального приближения из коэффициентов Ь2 и среднего арифметического свободных коэффициентов Ь1;

6. на каждом к-ом шаге итерации следим, чтобы значение коэффициента Ьк* не превосходило область значений Ь2 ± с * й2, где с — заранее фиксированный пороговый уровень; при достижении порогового уровня делаем случайный отскок от граничного значения внутрь области;

7. используя предложенную выше модификацию строим оптимальную многофакторную модель.

6 Результаты

/ \ ( 1 \

Для численной проверки генерировались наборы значений (хг,уг), хг = хг2 = I хг2 К г = 1,...,п.

\ ха ) \ х{3 )

Объем выборки п = 40, первая часть выборки моделируется для первого класса уг = — 1, г = 1,..., 20, для второго класса уг = 1, г = 21,..., 40. Признаки моделируются как независимые гауссовские X1 = X2 = (Х|, Х3) с одинаковой дисперсией равной 1.5 и математическими ожиданиями = 6, ЕХ3 = 3,

ЕХ| = —6, ЕX2 = —3. Всего проведено 1000 моделирований.

Таблица 2: Сравнение методов логистической регрессии

Новый алгоритм Алгоритм glm

E Оценка bi b2 b3 bi b2 b3

X1 =6 !3j 0.00000 -0.89443 -0.44721 0.00000 -0.89443 -0.44721

X1 =3 bj -0.00063 -0,89737 -0.44578 -0.09851 -0.80379 -0.39406

*22 = -6 SD(bj 0.05367 0.00859 0.02984 0.24789 0.14301 0.23874

X32 = -3 Abj -0.00082 0.00506 0.00677 0.08537 0,08189 0.07853

В таблице 2 используются следующие обозначения: fij — аналитические рассчитанные величины коэффициентов логистической регрессии; bj — математические средние оценок коэффициентов логистической регрессии; SD(bj — средние квадратические значения ошибок оценок коэффициентов; Abj — средние значения ошибок оценок коэффициентов. Нетрудно заметить, что разработанный алгоритм выигрывает по сравнению с glm.

Заключение

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

Список литературы

fl] Сох D.R., "The regression analysis of binary sequences (with discussion)"// J Roy Stat Soc B. 20 (2): 215-242, 1958.

[2] Walker S.H., Duncan D.B. "Estimation of the probability of an event as a function of several independent variables" // Biometrika. 54 (1/2): 167-178, 1967 doi: 10.2307/2333860. JSTOR 2333860.

[3] Truett J., Cornfield J., Kannel W, "A multivariate analysis of the risk of coronary heart disease in Framingham" // Journal of Chronic Diseases. 20 (7): 511-24, 1967 doi:10.1016/0021-9681(67)90082-3

[4] Rosenbaum P.R., Rubin D.B. "The Central Role of the Propensity Score in Observational Studies for Causal Effects"// Biometrika. 70 (1): 41-55, 1983 doi:10.1093/biomet/70.1.41

[5] Venkatesan R., Meng J.E., "A novel progressive learning technique for multi-class classification" // Neurocomputing. 207: 310-321, 2016. doi:10.1016/j.neucom.2016.05.006

[6] Hanley, James A., McNeil, Barbara J., "The Meaning and Use of the Area under a Receiver Operating Characteristic (ROC) Curve"// Radiology. 143 (1): 29-36, 1982. doi:10.1148/radiology.l43.1.7063747

[7] Cleveland W. S.Robust locally weighted regression and smoothing scatter plots //Journal of the American Statistical Association. 1979. Vol. 74, no. 368. Pp. 829-836.

[8] R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

[9] H. П. Васильев, А. А. Егоров, "Опыт расчета параметров логистической регрессии методом Ныото-на-Рафсона для оценки зимостойкости растений" // Матем. биология и биоинформ., 6:2 (2011), 190-199

[10] Бахвалов Н. С., Жидков Н. П., Кобельков Г. М., Численные методы, — 8-е изд., 90 ЭЛ. — М.:БИНОМ. Лаб. знаний, 2015. — 639с.

[11] Тырсин А. Н., Костин К. К. Оценивание логистической регрессии как экстремальная задача // Вестн. Том. гос. ун-та. Управление, вычислительная техника и информатика. 2017. № 40. С. 52-60.

[12] Kareva Yulia Е., Efendiev Vidady U., et al, Risk factors for the return of mitral insufficiency after coronary-artery bypass graft and reconstruction of the mitral valve in patients with ischemic cardiomyopathy // Circulation Pathology and Cardiac Surgery, Vol.23, Issue 2, P. 31-42, 2019 DOI: 10.21688/1681-3472-2019-231-42

[13] Simakova M.A., Ryzhkov A.V., et al, Perspectives of using pulmonary arterial stiffness indicators to evaluate the prognosis of patients with pulmonary arterial hypertension // Terapevticheskii arkhiv, V. 90 Issue 1, P. 86-92. DOI: 10.26442/terarkh201890186-92

Лукинов Виталий Леонидович — к.ф. м.н., ст. науч. сотр. Института вычислительной математики и математической геофизики СО РАН;

e-mail: vitaliy.lukinov@sscc.ru. Дата поступления — 30 апреля 2019 г.