Пошаговое построение логистической регрессии в Python
Логистическая регрессия — это алгоритм классификации машинного обучения, используемый для прогнозирования вероятности категориальной зависимой переменной. В логистической регрессии зависимая переменная является бинарной переменной, содержащей данные, закодированные как 1 (да, успех и т.п.) или 0 (нет, провал и т.п.). Другими словами, модель логистической регрессии предсказывает P(Y=1) как функцию X.
Условия логистической регрессии
- Бинарная логистическая регрессия требует, чтобы зависимая переменная также была бинарной.
- Для бинарной регрессии фактор уровня 1 зависимой переменной должен представлять желаемый вывод.
- Использоваться должны только значимые переменные.
- Независимые переменные должны быть независимы друг от друга. Это значит, что модель должна иметь малую мультиколлинеарность или не иметь её вовсе.
- Независимые переменные связаны линейно с логарифмическими коэффициентами.
- Логистическая регрессия требует больших размеров выборки.
Держа в уме все перечисленные условия, давайте взглянем на наш набор данных.
Данные
Набор данных взят с репозитория машинного обучения UCI и относится к прямым маркетинговым кампаниям (телефонный обзвон) португальского банковского учреждения. Цель классификации в прогнозировании успеха подписки клиента (1/0) на срочный депозит (переменная y). Загрузить этот набор данных можно здесь.
Эт и данные предоставляют информацию о клиентах банка, которая включает 41,188 записей и 21 поле.
Вводные переменные
- age — возраст (число);
- job — вид должности (категории: “admin” (администратор), “blue-collar” (рабочий), “entrepreneur” (мелкий предприниматель), “housemaid” (горничная), “management” (руководитель), “retired” (на пенсии), “self-employed” (самозанятый), “services” (сфера услуг), “student” (учащийся), “technician” (техник), “unemployed” (не трудоустроен), “unknown” (неизвестно));
- marital — семейное положение (категории: “divorced” (разведён), “married” (замужем/женат), “single” (холост/не замужем), “unknown” (неизвестно));
- education — образование (категории: “basic.4y”, “basic.6y”, “basic.9y” (базовое), “high.school” (высшая школа), “illiterate” (без образования), “professional.course” (профессиональные курсы), “university.degree” (университетская степень), “unknown” (неизвестно));
- default — имеет ли просроченные кредиты (категории: “no” (нет), “yes” (да), “unknown” (неизвестно));
- housing — имеет ли жилищный кредит (категории: “no” (нет), “yes” (да), “unknown” (неизвестно));
- loan — имеет ли личный кредит (категории: “no” (нет), “yes” (да), “unknown” (неизвестно));
- contact — вид связи (категории: “cellular” (мобильный), “telephone” (стационарный));
- month — месяц последнего обращения (категории: “jan” (январь), “feb” (февраль), “mar” (март), …, “nov” (ноябрь), “dec” (декабрь));
- day_of_week — день недели последнего обращения (категории: “mon” (понедельник), “tue” (вторник), “wed” (среда), “thu” (четверг), “fri” (пятница));
- duration — продолжительность последнего обращения в секундах (число). Важно: этот атрибут оказывает сильное влияние на вывод (например, если продолжительность=0, тогда y=’no’). Продолжительность не известна до момента совершения звонка, а по его завершению y будет, очевидно, известна. Следовательно этот вводный параметр должен включаться только в целях эталонного тестирования, для получения же реалистичной модели прогноза его следует исключать.
- campaign — число обращений, установленных в процессе этой кампании и для этого клиента (представлено числом и включает последнее обращение);
- pdays — число дней, прошедших с момента последнего обращения к клиенту во время предыдущей кампании (число; 999 означает, что ранее обращений не было);
- previous — число обращений, совершённых до этой кампании (число);
- poutcome — итоги предыдущей маркетинговой кампании (категории: “failure” (провал), “nonexistent” (несуществующий), “success” (успех));
- emp.var.rate — коэффициент изменения занятости (число);
- cons.price.idx — индекс потребительских цен (число);
- cons.conf.idx — индекс потребительского доверия (число);
- euribor3m — 3-х месячная европейская межбанковская ставка (число);
- nr.employed — количество сотрудников (число).
Прогнозируемая переменная (желаемая цель):
y —подписался ли клиент на срочный вклад (двоично: “1” означает “Да”, “0” означает “Нет”).
Колонка образования в наборе данных имеет очень много категорий, и нам нужно сократить их для оптимизации моделирования. В этой колонке представлены следующие категории:
Давайте сгруппируем “basic.4y”, “basic.9y” и “basic.6y” и назовём их “basic” (базовое).
Источник
Логистическая регрессия и ROC-анализ — математический аппарат
Математический аппарат и назначение бинарной логистической регрессии — популярного инструмента для решения задач регрессии и классификации. ROC-анализ тесно связан с бинарной логистической регрессией и применяется для оценки качества моделей: позволяет выбрать аналитику модель с наилучшей прогностической силой, проанализировать чувствительность и специфичность моделей, подобрать порог отсечения.
Введение
Логистическая регрессия — полезный классический инструмент для решения задачи регрессии и классификации. ROC-анализ — аппарат для анализа качества моделей. Оба алгоритма активно используются для построения моделей в медицине и проведения клинических исследований.
Логистическая регрессия получила распространение в скоринге для расчета рейтинга заемщиков и управления кредитными рисками. Поэтому, несмотря на свое «происхождение» из статистики, логистическую регрессию и ROC-анализ почти всегда можно увидеть в наборе Data Mining алгоритмов.
Логистическая регрессия
Логистическая регрессия — это разновидность множественной регрессии, общее назначение которой состоит в анализе связи между несколькими независимыми переменными (называемыми также регрессорами или предикторами) и зависимой переменной. Бинарная логистическая регрессия применяется в случае, когда зависимая переменная является бинарной (т.е. может принимать только два значения). С помощью логистической регрессии можно оценивать вероятность того, что событие наступит для конкретного испытуемого (больной/здоровый, возврат кредита/дефолт и т.д.).
Все регрессионные модели могут быть записаны в виде формулы:
y = F (x_1,\, x_2, \,\dots, \, x_n)
В множественной линейной регрессии предполагается, что зависимая переменная является линейной функцией независимых переменных, т.е.:
Можно ли ее использовать для задачи оценки вероятности исхода события? Да, можно, вычислив стандартные коэффициенты регрессии. Например, если рассматривается исход по займу, задается переменная y со значениями 1 и 0, где 1 означает, что соответствующий заемщик расплатился по кредиту, а 0, что имел место дефолт.
Однако здесь возникает проблема: множественная регрессия не «знает», что переменная отклика бинарна по своей природе. Это неизбежно приведет к модели с предсказываемыми значениями большими 1 и меньшими 0. Но такие значения вообще не допустимы для первоначальной задачи. Таким образом, множественная регрессия просто игнорирует ограничения на диапазон значений для y .
Для решения проблемы задача регрессии может быть сформулирована иначе: вместо предсказания бинарной переменной, мы предсказываем непрерывную переменную со значениями на отрезке [0,1] при любых значениях независимых переменных. Это достигается применением следующего регрессионного уравнения (логит-преобразование):
где P — вероятность того, что произойдет интересующее событие e — основание натуральных логарифмов 2,71…; y — стандартное уравнение регрессии.
Зависимость, связывающая вероятность события и величину y , показана на следующем графике (рис. 1):
Рис. 1 — Логистическая кривая
Поясним необходимость преобразования. Предположим, что мы рассуждаем о нашей зависимой переменной в терминах основной вероятности P , лежащей между 0 и 1. Тогда преобразуем эту вероятность P :
P’ = \log_e \Bigl(\frac
<1-P>\Bigr)
Это преобразование обычно называют логистическим или логит-преобразованием. Теоретически P’ может принимать любое значение. Поскольку логистическое преобразование решает проблему об ограничении на 0-1 границы для первоначальной зависимой переменной (вероятности), то эти преобразованные значения можно использовать в обычном линейном регрессионном уравнении. А именно, если произвести логистическое преобразование обеих частей описанного выше уравнения, мы получим стандартную модель линейной регрессии.
Существует несколько способов нахождения коэффициентов логистической регрессии. На практике часто используют метод максимального правдоподобия. Он применяется в статистике для получения оценок параметров генеральной совокупности по данным выборки. Основу метода составляет функция правдоподобия (likehood function), выражающая плотность вероятности (вероятность) совместного появления результатов выборки
L\,(Y_1,\,Y_2,\,\dots,\,Y_k;\,\theta) = p\,(Y_1;\, \theta)\cdot\dots\cdotp\,p\,(Y_k;\,\theta)
Согласно методу максимального правдоподобия в качестве оценки неизвестного параметра принимается такое значение \theta=\theta(Y_1,…,Y_k) , которое максимизирует функцию L .
Нахождение оценки упрощается, если максимизировать не саму функцию L , а натуральный логарифм ln(L) , поскольку максимум обеих функций достигается при одном и том же значении \theta :
L\,*\,(Y;\,\theta) = \ln\,(L\,(Y;\,\theta)\,) \rightarrow \max
В случае бинарной независимой переменной, которую мы имеем в логистической регрессии, выкладки можно продолжить следующим образом. Обозначим через P_i вероятность появления единицы: P_i=Prob(Y_i=1) . Эта вероятность будет зависеть от X_iW , где X_i — строка матрицы регрессоров, W — вектор коэффициентов регрессии:
Логарифмическая функция правдоподобия равна:
где I_0 , I_1 — множества наблюдений, для которых Y_i=0 и Y_i=1 соответственно.
Можно показать, что градиент g и гессиан H функции правдоподобия равны:
g = \sum_i (Y_i\,-\,P_i)\,X_i
H=-\sum_i P_i\,(1\,-\,P_i)\,X_i^T\,X_i\,\leq 0
Гессиан всюду отрицательно определенный, поэтому логарифмическая функция правдоподобия всюду вогнута. Для поиска максимума можно использовать метод Ньютона, который здесь будет всегда сходиться (выполнено условие сходимости метода):
Логистическую регрессию можно представить в виде однослойной нейронной сети с сигмоидальной функцией активации, веса которой есть коэффициенты логистической регрессии, а вес поляризации — константа регрессионного уравнения (рис. 2).
Рис. 2 — Представление логистической регрессии в виде нейронной сети
Однослойная нейронная сеть может успешно решить лишь задачу линейной сепарации. Поэтому возможности по моделированию нелинейных зависимостей у логистической регрессии отсутствуют. Однако для оценки качества модели логистической регрессии существует эффективный инструмент ROC-анализа, что является несомненным ее преимуществом.
Для расчета коэффициентов логистической регрессии можно применять любые градиентные методы: метод сопряженных градиентов, методы переменной метрики и другие.
ROC-анализ
ROC-кривая (Receiver Operator Characteristic) — кривая, которая наиболее часто используется для представления результатов бинарной классификации в машинном обучении. Название пришло из систем обработки сигналов. Поскольку классов два, один из них называется классом с положительными исходами, второй — с отрицательными исходами. ROC-кривая показывает зависимость количества верно классифицированных положительных примеров от количества неверно классифицированных отрицательных примеров.
В терминологии ROC-анализа первые называются истинно положительным, вторые — ложно отрицательным множеством. При этом предполагается, что у классификатора имеется некоторый параметр, варьируя который, мы будем получать то или иное разбиение на два класса. Этот параметр часто называют порогом, или точкой отсечения (cut-off value). В зависимости от него будут получаться различные величины ошибок I и II рода.
В логистической регрессии порог отсечения изменяется от 0 до 1 — это и есть расчетное значение уравнения регрессии. Будем называть его рейтингом.
Для понимания сути ошибок I и II рода рассмотрим четырехпольную таблицу сопряженности (confusion matrix), которая строится на основе результатов классификации моделью и фактической (объективной) принадлежностью примеров к классам.
Модель | Фактически положительно | Фактически отрицательно |
---|---|---|
Положительно | TP | FP |
Отрицательно | FN | TN |
- TP (True Positives) — верно классифицированные положительные примеры (так называемые истинно положительные случаи).
- TN (True Negatives) — верно классифицированные отрицательные примеры (истинно отрицательные случаи).
- FN (False Negatives) — положительные примеры, классифицированные как отрицательные (ошибка I рода). Это так называемый «ложный пропуск» — когда интересующее нас событие ошибочно не обнаруживается (ложно отрицательные примеры).
- FP (False Positives) — отрицательные примеры, классифицированные как положительные (ошибка II рода). Это ложное обнаружение, т.к. при отсутствии события ошибочно выносится решение о его присутствии (ложно положительные случаи).
Что является положительным событием, а что — отрицательным, зависит от конкретной задачи. Например, если мы прогнозируем вероятность наличия заболевания, то положительным исходом будет класс «Больной пациент», отрицательным — «Здоровый пациент». И наоборот, если мы хотим определить вероятность того, что человек здоров, то положительным исходом будет класс «Здоровый пациент», и так далее.
При анализе чаще оперируют не абсолютными показателями, а относительными — долями (rates), выраженными в процентах:
- Доля истинно положительных примеров (True Positives Rate): TPR = \frac
\,\cdot\,100 \,\% - Доля ложно положительных примеров (False Positives Rate): FPR = \frac
\,\cdot\,100 \,\%
Введем еще два определения: чувствительность и специфичность модели. Ими определяется объективная ценность любого бинарного классификатора.
Чувствительность (Sensitivity) — это и есть доля истинно положительных случаев:
Специфичность (Specificity) — доля истинно отрицательных случаев, которые были правильно идентифицированы моделью:
Заметим, что FPR=100-Sp
Попытаемся разобраться в этих определениях.
Модель с высокой чувствительностью часто дает истинный результат при наличии положительного исхода (обнаруживает положительные примеры). Наоборот, модель с высокой специфичностью чаще дает истинный результат при наличии отрицательного исхода (обнаруживает отрицательные примеры). Если рассуждать в терминах медицины — задачи диагностики заболевания, где модель классификации пациентов на больных и здоровых называется диагностическим тестом, то получится следующее:
- Чувствительный диагностический тест проявляется в гипердиагностике — максимальном предотвращении пропуска больных.
- Специфичный диагностический тест диагностирует только доподлинно больных. Это важно в случае, когда, например, лечение больного связано с серьезными побочными эффектами и гипердиагностика пациентов не желательна.
ROC-кривая получается следующим образом:
Для каждого значения порога отсечения, которое меняется от 0 до 1 с шагом d_x (например, 0,01) рассчитываются значения чувствительности Se и специфичности Sp . В качестве альтернативы порогом может являться каждое последующее значение примера в выборке.
Строится график зависимости: по оси Y откладывается чувствительность Se , по оси X — FPR=100-Sp — доля ложно положительных случаев.
Канонический алгоритм построения ROC-кривой
Входы: L — множество примеров f[i] — рейтинг, полученный моделью, или вероятность того, что i -й пример имеет положительный исход; min и max — минимальное и максимальное значения, возвращаемые f ; d_x — шаг; P и N — количество положительных и отрицательных примеров соответственно.
- t=min
- повторять
- FP=TP=0
- для всех примеров i принадлежит L <
- если f[i]>=t тогда // этот пример находится за порогом
- если i положительный пример тогда
- иначе // это отрицательный пример
- >
- Se=TP/P*100
- point=FP/N // расчет (100 минус Sp )
- Добавить точку (point, Se) в ROC-кривую
- t=t+d_x
- пока (t>max)
В результате вырисовывается некоторая кривая (рис. 3).
Рис. 3 — ROC-кривая
График часто дополняют прямой y=x .
Заметим, что имеется более экономичный способ расчета точек ROC-кривой, чем тот, который приводился выше, т.к. его вычислительная сложность нелинейная и равна O(n^2) : для каждого порога необходимо «пробегать» по записям и каждый раз рассчитывать TP и FP . Если же двигаться вниз по набору данных, отсортированному по убыванию выходного поля классификатора (рейтингу), то можно за один проход вычислить значения всех точек ROC-кривой, последовательно обновляя значения TP и FP .
Для идеального классификатора график ROC-кривой проходит через верхний левый угол, где доля истинно положительных случаев составляет 100% или 1,0 (идеальная чувствительность), а доля ложно положительных примеров равна нулю. Поэтому чем ближе кривая к верхнему левому углу, тем выше предсказательная способность модели. Наоборот, чем меньше изгиб кривой и чем ближе она расположена к диагональной прямой, тем менее эффективна модель. Диагональная линия соответствует «бесполезному» классификатору, т.е. полной неразличимости двух классов.
При визуальной оценке ROC-кривых расположение их относительно друг друга указывает на их сравнительную эффективность. Кривая, расположенная выше и левее, свидетельствует о большей предсказательной способности модели. Так, на рис. 4 две ROC-кривые совмещены на одном графике. Видно, что модель «A» лучше.
Рис. 4 — Сравнение ROC-кривых
Визуальное сравнение кривых ROC не всегда позволяет выявить наиболее эффективную модель. Своеобразным методом сравнения ROC-кривых является оценка площади под кривыми. Теоретически она изменяется от 0 до 1,0, но, поскольку модель всегда характеризуются кривой, расположенной выше положительной диагонали, то обычно говорят об изменениях от 0,5 («бесполезный» классификатор) до 1,0 («идеальная» модель).
Эта оценка может быть получена непосредственно вычислением площади под многогранником, ограниченным справа и снизу осями координат и слева вверху — экспериментально полученными точками (рис. 5). Численный показатель площади под кривой называется AUC (Area Under Curve). Вычислить его можно, например, с помощью численного метода трапеций:
AUC = \int f(x)\,dx = \sum_i \Bigl[ \frac
Рис. 5 — Площадь под ROC-кривой
С большими допущениями можно считать, что чем больше показатель AUC , тем лучшей прогностической силой обладает модель. Однако следует знать, что:
- показатель AUC предназначен скорее для сравнительного анализа нескольких моделей;
- AUC не содержит никакой информации о чувствительности и специфичности модели.
В литературе иногда приводится следующая экспертная шкала для значений AUC , по которой можно судить о качестве модели:
Интервал AUC | Качество модели |
---|---|
0,9-1,0 | Отличное |
0,8-0,9 | Очень хорошее |
0,7-0,8 | Хорошее |
0,6-0,7 | Среднее |
0,5-0,6 | Неудовлетворительное |
Идеальная модель обладает 100% чувствительностью и специфичностью. Однако на практике добиться этого невозможно, более того, невозможно одновременно повысить и чувствительность, и специфичность модели. Компромисс находится с помощью порога отсечения, т.к. пороговое значение влияет на соотношение Se и Sp . Можно говорить о задаче нахождения оптимального порога отсечения (optimal cut-off value).
Порог отсечения нужен для того, чтобы применять модель на практике: относить новые примеры к одному из двух классов. Для определения оптимального порога нужно задать критерий его определения, т.к. в разных задачах присутствует своя оптимальная стратегия. Критериями выбора порога отсечения могут выступать:
- Требование минимальной величины чувствительности (специфичности) модели. Например, нужно обеспечить чувствительность теста не менее 80%. В этом случае оптимальным порогом будет максимальная специфичность (чувствительность), которая достигается при 80% (или значение, близкое к нему «справа» из-за дискретности ряда) чувствительности (специфичности).
- Требование максимальной суммарной чувствительности и специфичности модели, т.е. Cutt\underline<\,\,\,>off_o = \max_k (Se_k\,+\,Sp_k)
- Требование баланса между чувствительностью и специфичностью, т.е. когда Se \approx Sp : Cutt\underline<\,\,\,>off_o = \min_k \,\bigl |Se_k\,-\,Sp_k \bigr |
Второе значение порога обычно предлагается пользователю по умолчанию. В третьем случае порог есть точка пересечения двух кривых, когда по оси X откладывается порог отсечения, а по оси Y — чувствительность или специфичность модели (рис. 6).
Рис. 6 — «Точка баланса» между чувствительностью и специфичностью
Существуют и другие подходы, когда ошибкам I и II рода назначается вес, который интерпретируется как цена ошибок. Но здесь встает проблема определения этих весов, что само по себе является сложной, а часто не разрешимой задачей.
Источник
Применение логистической регрессии: задача о программистах
Подобные инструменты доступны в пакете Углубленные методы анализа — модули Нелинейное оценивание и Обобщенные линейные/нелинейные модели. Все, что нужно сделать — это задать переменные и выбрать способ оценивания целевой функции.
Приведем пример такого анализа. Основная часть данных для данного примера взята из работы Neter, Wasserman, Kutner (1985). Однако отметим, что они использовали для подгонки линейную регрессионную модель.
Предположим, что вы хотите проверить, правда ли, что стаж работы помогает программистам в написании сложных программ, если на написание отпущен ограниченный промежуток времени. Для исследования были выбраны двадцать пять программистов с различным стажем работы (выраженным в месяцах). Их попросили написать сложную компьютерную программу за определенный промежуток времени.
Бинарная переменная отклика принимала значение 1, если программист справился с поставленной задачей, и 0, если нет.
Эти исходные данные выглядят следующим образом:
Шаг 1. Визуализация
Первым шагом для любого анализа является осознание структуры представленных данных. У нас есть таблица с двумя переменными. Для начала посмотрим, как распределен стаж работы кандидатов — построим гистограмму для переменной EXPERENCE.
Действие 1. Выделите переменную EXPERENCE и правым кликом вызовите контекстное меню. В этом меню выберете пункт Графики блоковых данных -> гистограммы: все столбцы.
Будет отображена гистограмма для переменной EXPERENCE. Она выглядит следующим образом:
Мы видим, что опыт работы для программистов распределен довольно равномерно. Представлены как опытные, так и неопытные кандидаты и их примерно одинаковое число.
Насколько эффективно программисты справлялись с заданием? Построим диаграмму рассеяния.
Действие 2. Выберете опцию Графика -> Диаграммы рассеяния. Будет отображена панель задания параметров диаграммы рассеяния.
Действие 3. Выберете переменные для построения диаграммы рассеяния. Для этого нажмите на кнопку Переменные и задайте EXPERENCE — как переменную по оси х и SUCCESS — как переменную оси y.
Оставьте остальные параметры по умолчанию. Нажмите на кнопку ОК — будет отображена диаграмма рассеяния. Для данного случая диаграмма примет такой вид:
На диаграмме рассеяния выделяются два облака точек. Одно — вблизи программистов с небольшим опытом и проваливших задание, второе — вблизи программистов с обширным опытом и выполнивших задание. Гипотеза подтверждается прямо на графике!
Шаг 2. Задание модели.
Теперь строго подтвердим наши догадки. Построим логистическую регрессию.
Действие 1. Выберете пункт меню Анализ -> Углубленные методы анализа -> Нелинейное оценивание. В появившемся окне выберете Логит регрессия. Стартовая панель модуля выглядит следующим образом:
Действие 2. Выберем переменную SUCCESS как зависимую и EXPERNCE как независимую. Для этого нажмите на кнопку Переменные.
Программа автоматически выберет коды зависимой переменной.
Шаг 3. Задание метода оценивания.
После нажатия на кнопку ОК на стартовой панели будет отображен диалог определения оценивания модели. Здесь вы можете выбрать метод оценивания, уточнить критерий сходимости, начальные значения и т.д. Вы можете также выбрать вычисление (с использованием метода конечных разностей) асимптотических стандартных ошибок оценок параметров. Панель оценивания модели выглядит следующим образом.
Действие 1. На вкладке Дополнительно выберете Метод оценивания — Квази- Ньютоновский. Установите опцию Асимптотические стандартные ошибки на Вкл.
Действие 2. Нажмите на кнопку ОК, чтобы начать вычисления. При этом будут отображаться результаты итераций. Если процесс сойдется, переходите на следующий шаг. Если нет, выполните следующее действие.
Действие 3. На вкладке Дополнительно попробуйте изменить метод оценивания или же Начальные значения. Методов оценки много и скорее всего итерационный процесс сойдется.
Когда процесс сошелся, можно переходить к Шагу 4.
Шаг 4. Просмотр результатов.
После проведения вычислений будет отображена панель диалога просмотра результатов. Здесь собрана вся информация, касающаяся построенной модели и результатов оценивания. Для данного примера окно выглядит следующим образом.
Действие 1. На панели диалога отображения результатов содержится р-уровень гипотезы. Если этот р-уровень менее 5%, то модель значима.
В данном случае р-уровень гипотезы оказался ниже 5% — значение статистики хи-квадрат для разницы между текущей моделью и моделью, содержащей лишь свободный член, высоко значимо. Поэтому можно заключить, что стаж работы влияет на успехи программиста в выполнении поставленной задачи. Результаты работы собраны в виде нескольких таблиц.
Шаг 5. Интерпретация результатов
Действие 1. Выберем опцию Параметры и стандартные ошибки. Рассмотрим таблицу, в которой содержатся данные об оценках регрессионных коэффициентов. В таблице результатов ниже оба параметра имеют уровень значимости p<.05.
В принципе, оценки параметров могут быть проинтерпретированы, как и в случае стандартной линейной регрессионной модели, т.е. в терминах свободного члена (Const.B) и углового коэффициента (EXPERNCE). По существу, результаты исследования показывают, что продолжительность имеющегося стажа существенно влияет на успешное проведение порученной работы по программированию.
Однако, оцениваемые параметры относятся к предсказанию логит-преобразования (вычисляемого как log[p/(1-p)]), а не самой вероятности (p), определяющей возможность успеха или неудачи. Логит преобразование принимает значения от минус до плюс бесконечности, когда значения вероятности p пробегают отрезок от 0 до 1.
Действие 2. На вкладке Быстрый выберите опцию Наблюдаемые, предсказанные и значения остатков. Напомним, что регрессионная модель логит гарантирует, что предсказанные значения всегда будут находиться внутри отрезка [0,1]. Поэтому вы можете рассматривать полученные значения как вероятности. Например, предсказанная вероятность успеха для второго программиста (Henry) равна (.84).
Шаг 6. Оценка качества модели
Действие 1. На вкладке Дополнительно нажмите на кнопку Классификация. Будет отображена таблица с результатами классификации.
Оценить качество построенной модели можно, если оценить параметр Отношение несогласия. Выведем на экран таблицу с числом наблюдений, которые были правильно и неправильно классифицированы в соответствии с полученной моделью.
Все наблюдения с предсказанными значениями (вероятностью) меньше или равными .5 классифицируются как неудача — Failure, остальные, с предсказываемыми значениями больше .5, классифицируются как успех — Success. Отношение несогласия вычисляется как отношение произведения чисел правильно расклассифицированных наблюдений к произведению чисел неправильно расклассифицированных. Отношение несогласия больше 1 показывает, что построенная классификация лучше, чем, если бы мы просто провели классификацию наугад.
Однако следует помнить, что наша классификация была подобрана так, чтобы максимизировать вероятность успеха для уже полученных данных, которым соответствовал успех. Поэтому не следует заранее рассчитывать на хорошую классификацию, если вы в будущем примените нашу модель к новым наблюдениям (как уже говорилось, логит регрессионная модель сильно подвержена излишней подгонке).
Осталось отметить, что при использовании логит регрессионной модели необходимо дополнять исследования другими методами, например, деревьями классификации.
Источник
Блог странного учёного
NB: Этот материал представляет собой сокращённый перевод публикации R Data Analysis Examples: Logit Regression, http://www.ats.ucla.edu/stat/r/dae/logit.htm. Большая часть описанных манипуляций давно автоматизирована в пакете epicalc (см. напр. http://donbas-socproject.blogspot.com/2011/03/blog-post_25.html), однако настоящее описание позволяет взглянуть «под капот» и разобраться в логике построения вывода.
Логистическая регрессия, называемая также логит-моделью, используется для моделирования двоичных зависимых переменных. В этом методе логарифм шансов наступления исследуемого события представляется линейной комбинацией переменных-предикторов.
Для анализа мы дополнительно воспользуемся пакетами aod и ggplot2 , подключаемых командами:
Пример. Исследователь ищет, как переменные GRE (Graduate Record Exam scores — оценки во время обучения в вузе), GPA (grade point average — средний балл) и престиж вуза влияют на поступление в аспирантуру. Т. о. искомая переменная, поступил или не поступил (admit/don’t admit), является бинарной.
В описанном ниже анализе мы воспользовались сгенерированными гипотетическими данными, которые можно получить, используя R :
mydata
head(mydata)## посмотрим на первые строки данных
## admit gre gpa rank
## 1 0 380 3.61 3
## 2 1 660 3.67 3
## 3 1 800 4.00 1
## 4 1 640 3.19 4
## 5 0 520 2.93 4
## 6 1 760 3.00 2
Этот масив содержит бинарную зависимую переменную admit, а также три предиктора: gre, gpa and rank. Мы воспользуемся переменными gre и gpa как непрерывными. Переменная rank принимает значения от 1 до 4. У вузов с рангом 1 самый высокий престиж, а у вузов с рангом 4 — самый низкий.
Методы анализа, которые вы также можете использовать
Пробит-регрессия даёт результаты, похожие на логистическую. Выбор пробит или логит-моделирования в значительной степени зависит от личных предпочтений исследователя.
OLS-регрессия, используемая с бинарной зависимой переменной, известна как линейно-вероятностная модель (linear probability model) и может применяться для описания условных вероятностей. Тем не менее, ошибки (т. н. остатки — residuals), даваемые этим методом, не соответствуют предположению о гомоскедастичности (равной дисперсии) и нормальном распределении, заложенном в OLS-регрессии, что приводит к неправильным значениям стандартных ошибок и ложным результатам тестирования гипотез (подробнее см. Long, 1997, p.
Дискриминантный анализ двух групп. Многомерный метод для бинарной зависимой переменной.
Hotelling’s T2. Зависимая переменная 0/1 преобразуется в группирующую переменную, а предикторы — в зависимые переменные. Это позволяет выполнить общий тест значимости, но не даёт индивидуальных коэффициентов для каждой переменной, соответственно остаётся неясным вклад каждого такого «предиктора» в эффект прочих «предикторов».
Использование логистической модели
Код, приведённый ниже, выполняет логистическое моделирование с использованием функции glm() function, но сначала мы превратим переменную rank в фактор, чтобы показать, что она должна трактоваться как категориальная переменная:
Поскольку мы дали нашей модели имя mylogit , R не выдал результаты регресии на экран. Это, однако, легко сделать командой summary() :
## Call:
## glm(formula = admit
gre + gpa + rank, family = «binomial»,
## data = mydata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.627 -0.866 -0.639 1.149 2.079
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.98998 1.13995 -3.50 0.00047 ***
## gre 0.00226 0.00109 2.07 0.03847 *
## gpa 0.80404 0.33182 2.42 0.01539 *
## rank2 -0.67544 0.31649 -2.13 0.03283 *
## rank3 -1.34020 0.34531 -3.88 0.00010 ***
## rank4 -1.55146 0.41783 -3.71 0.00020 ***
## —
## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ‘ 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 499.98 on 399 degrees of freedom
## Residual deviance: 458.52 on 394 degrees of freedom
## AIC: 470.5
##
## Number of Fisher Scoring iterations: 4
В приведённых результатах вначале мы видим саму модель (секция «Call») и заданные опции.
Затем видим отклонение остатков (deviance residuals), измеряющий подгонку модели. Ниже мы обсудим, как пользоваться этой строкой, чтобы оценить соответствие модели.
Следующая часть результатов показывает коэффициенты, их стандартные ошибки, z-статистики (называемые иногда z-статистиками Вальда — коэффициенты, делённые на стандартные ошибки) и соответствующие p-значения. Очевидно, что переменные gre и gpa суть статистически значимы, точно также как и три уровня переменной rank . Коэффициенты логистической регресии указывают, насколько изменится вероятность искомого события при увеличении соответствующего предиктора на одну единицу.
Так, возрастание на одну единицу gre увеличивает логарифм шансов поступить в аспирантуру (по сравнению с непоступлением) на 0.002. Возрастание на одну единицу gpa , увеличивает логарифм шансов поступления на 0.804.
Уровни переменной rank несколько отличаются интерпретацией. Напр., у выпускников вуза с рангом 2 (по сравнению с вузом первого ранга), меньше логарифм шансов быть в аспирантуре на 0.675.
Под таблицей коэффициентов находятся служебные данные. Чуть позже мы покажем пример, как можно использовать эти цифры для оценки качества модели.
Мы можем воспользоваться функцией confint() , чтобы получить доверительные интервалы рассчитанных коэффициентов.
## 2.5 % 97.5 %
## (Intercept) -6.2716202 -1.792547
## gre 0.0001376 0.004436
## gpa 0.1602959 1.464143
## rank2 -1.3008888 -0.056746
## rank3 -2.0276713 -0.670372
## rank4 -2.4000265 -0.753543
Мы можем протестировать общее влияние переменной rank функцией wald.test() из библиотеки aod . Порядок, в котором коєффициенты размещены в таблице коэффициентов, повторяет порядок слагаемых в модели. Это важно, поскольку функция wald.test() применяется к коэффициентам согласно этому порядку (опция Terms указывает R , какие именно слагаемые из модели должны быть протестированы, в данном случае это 4, 5 и 6 суть три слагаемых, относящихся к уровням переменной rank ).
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), Terms = 4:6)
## Chi-squared test:
## X2 = 20.9, df = 3, P(> X2) = 0.00011
Значение теста хи-квадрат 20.9 с тремя степенями свободы ассоциировано с p-значением 0.00011, показывающим, что общее влияние переменной rank является статистически значимым.
Мы также можем протестировать другие гипотезы о различиях в коэффициентах разных уровней переменной rank . Напр., ниже мы тестируем предположение, что коэффициент для rank=2 равен коэффициенту для rank=3. Первая строка приведённого кода создаёт вектор l, описывающий тот тест, который мы намерены выполнить. В нашем случае мы тестируем разницу слагаемых для rank=2 и rank=3 (т. е., и слагаемое в модели). Для контраста мы умножаем один из них на 1, а другой на −1. Прочие слагаемые модели не участвуют в тесте, поэтому мы умножаем их на 0. Вторая линия кода использует равенство L=l для того, чтобы сказать R, что мы решили выполнить тест с использованием вектора l.
l
wald.test(b = coef(mylogit), Sigma = vcov(mylogit), L = l)
## Chi-squared test:
## X2 = 5.5, df = 1, P(> X2) = 0.019
Значение теста хи-квадрат 5.5 с одной степенью свободы ассоциировано с p-значением 0.019, показывающим, что разница между коэффициентами для rank=2 и rank=3 статистически значима.
Вы также можете экспоненцировать коэффициенты и интерпретировать их как отношения шансов (OR). R сделает это, если вы укажете ему, что откуда взять коэффициенты — из coef(mylogit) . Аналогично вы можете получить и доверительные интервалы отношений шансов. Чтобы вставить это всё в одну таблицу, мы используем команду cbind() :
## (Intercept) gre gpa rank2 rank3 rank4
## 0.0185 1.0023 2.2345 0.5089 0.2618 0.2119
exp(cbind(OR = coef(mylogit), confint(mylogit)))
## OR 2.5 % 97.5 %
## (Intercept) 0.0185 0.001889 0.1665
## gre 1.0023 1.000138 1.0044
## gpa 2.2345 1.173858 4.3238
## rank2 0.5089 0.272290 0.9448
## rank3 0.2618 0.131642 0.5115
## rank4 0.2119 0.090716 0.4707
Теперь мы можем сказать, что при возрастании на одну единицу переменной gpa отношение шансов поступить в аспирантуру (в сравнении с непоступлением) возрастает в 2.23 раза.
Вы можете использовать предсказанные вероятности, чтобы лучше понять эту модель. Предсказанные вероятности можно рассчитать как для категориальных, так и для непрерывных предикторов. Для того, чтобы создать предсказываемые вероятности, нам сначала нужно создать новую таблицу с данными, в которой будут значения независимых переменных такие, для которых нужно произвести предсказание.
Мы начнём с расчёта презсказанных вероятностей поступления для каждого значения престижности вуза (переменная rank ), удерживая переменные gre и gpa на уровне их средних. Создадим и посмотрим на таблицу данных:
## gre gpa rank
## 1 587.7 3.39 1
## 2 587.7 3.39 2
## 3 587.7 3.39 3
## 4 587.7 3.39 4
Заметьте, что в новой таблице имена переменных должны быть такими же, как и в вашей регрессионной модели (напр., среднее для gre должно быть названо gre). Теперь, имея таблицу, мы можем поручить R создать пресказанные вероятности. Первая линия кода очень компактна — значения rankP должны быть предсказаны [функция predict() ] на основе результатов анализа, содержащихся в mylogit и значений предикторов из таблицы newdata1 .
## gre gpa rank rankP
## 1 587.7 3.39 1 0.5166
## 2 587.7 3.39 2 0.3523
## 3 587.7 3.39 3 0.2186
## 4 587.7 3.39 4 0.1847
Из приведённых результатов следует, что предсказанные вероятности быть принятым в аспирантуру суть 0.52 для студентов из наиболее престижных вузов (rank=1), но 0.18 для студентов из наименее престижных (rank=4) при одинаковых средних показателях успеваемости gre и gpa. Аналогично можно варьировать значения переменных gre и rank . Если мы намерены нарисовать эти зависимости, мы должны создать 100 значений gre в интервале 200 и 800 для каждого значения переменной rank (т. е., 1, 2, 3 и 4):
Код, генерирующий предсказанные вероятности, аналогичен приведённому выше.
newdata2
newdata2
head(newdata2)
## gre gpa rank fit se.fit residual.scale UL LL
## 1 200.0 3.39 1 -0.8115 0.5148 1 0.5492 0.1394
## 2 206.1 3.39 1 -0.7978 0.5091 1 0.5499 0.1424
## 3 212.1 3.39 1 -0.7840 0.5034 1 0.5505 0.1454
## 4 218.2 3.39 1 -0.7703 0.4978 1 0.5512 0.1485
## 5 224.2 3.39 1 -0.7566 0.4922 1 0.5519 0.1517
## 6 230.3 3.39 1 -0.7429 0.4866 1 0.5525 0.1549
## PredictedProb
## 1 0.3076
## 2 0.3105
## 3 0.3134
## 4 0.3164
## 5 0.3194
## 6 0.3224
Рисунок полученной модели делается так:
ggplot(newdata2, aes(x = gre, y = PredictedProb)) + geom_ribbon(aes(ymin = LL,
ymax = UL, fill = rank), alpha = 0.2) + geom_line(aes(colour = rank), size = 1)
Мы можем также захотеть узнать, насколько хорошо наша модель описывает реальные данные. Это может быть особенно полезно, если мы сравниваем конкурирующие модели. Результаты, выводимые на экран командой summary(mylogit) , включают показатели подгонки (они выводятся сразу под коэффициентами), в частности разброс остатков нулевой и данной моделей, а также AIC. Один из параметров подгонки — значимость модели в целом. Этот тест показывает, действительно ли модель с предикторами описывает реальные данные значимо лучше, чем модель без предикторов (т. н. нулевая модель). Тестовая статистика — это разница между отклонениями остатков в моделях с предикторами и без них. Эта статистика описывается распределением хи-квадрат с числом степеней свободы равном разнице в степенях свободы тестируемой и нулевой моделью. Чтобы вычислить эту величину, мы можем воспользоваться командой:
with(mylogit, null.deviance — deviance)
А разница степеней свободы может быть получена так:
with(mylogit, df.null — df.residual)
Наконец, p-значение вычисляется так:
with(mylogit, pchisq(null.deviance — deviance, df.null — df.residual, lower.tail = FALSE))
Хи-квадрат 41.46 с 5 степенями свободы и ассоциированным p-значением, меньшим 0.001, показывают нам, что наша модель, в целом, описывает данные значительно лучше, чем нулевая модель. Этот тест иногда называют «likelihood ratio test».
Источник