Логистическая регрессия в R: glm() с family = binomial

Когда исход бинарный

Линейная регрессия предсказывает число. Но многие вопросы, достойные моделирования, бинарны: уйдёт ли клиент, выздоровеет ли пациент, кликнут ли по письму. Подгонка прямой к исходу 0/1 ломается немедленно — прямая бодро предсказывает вероятности −0.3 или 1.4, что является бессмыслицей.

Логистическая регрессия исправляет это, моделируя вероятность исхода через преобразование логарифма шансов (логит): log(p / (1 − p)) = свободный член + наклон × x. Шкала логарифма шансов пробегает всю числовую прямую, поэтому линейное уравнение ложится на неё естественно, а обратное отображение сжимает каждый прогноз в интервал (0, 1) по знакомой S-образной кривой. Плата за трюк: коэффициенты живут на шкале логарифма шансов, и вся игра в чтение логистической регрессии состоит в переводе их обратно во что-то понятное людям.

Подгонка: glm() с family = binomial

glm() (обобщённая линейная модель) — старший брат lm(); family = binomial выбирает логистическую регрессию. В mtcars переменная am фиксирует тип коробки передач (1 = механическая, 0 = автоматическая) — склонны ли экономичные машины быть механическими?

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
summary(fit)

Две вещи до чтения вывода. Во-первых, family = binomial не является необязательным: опустите его — и glm() молча подгонит обычный метод наименьших квадратов. Во-вторых, отклик должен быть бинарным: 0/1, логический или фактор с двумя уровнями (R моделирует вероятность второго уровня).

Теперь сводка, блок за блоком:

• Coefficients — Estimate при mpg составляет около 0.31, и это наклон в логарифмах шансов: каждая дополнительная миля на галлон добавляет 0.31 к логарифму шансов быть механической. Положительный знак означает «повышает вероятность», отрицательный — «понижает»; кроме знака, ничья интуиция на этой шкале не работает, поэтому и существует следующий раздел.
• z value и Pr(>|z|) — та же логика, что и у t-критериев регрессии (Estimate ÷ Std. Error, затем p-значение для вопроса «может ли это быть нулём?»), просто с нормальным приближением, отсюда z вместо t. Здесь p ≈ 0.011: связь между расходом и типом коробки вряд ли является шумом.
• Null deviance против Residual deviance — девианс это мера плохости подгонки в мире glm (меньше = лучше). Null deviance (43.2 при 31 df) относится к модели только со свободным членом; residual deviance (29.7 при 30 df) — к вашей. Падение примерно на 13.6 за одну степень свободы предиктора является glm-аналогом фразы «R-квадрат вырос».
• AIC — оценка для сравнения моделей, балансирующая качество подгонки и сложность; побеждает меньшая. Сама по себе бессмысленна, полезна при сравнении кандидатов на одних и тех же данных.

От логарифмов шансов к отношениям шансов: exp(coef())

Экспоненцирование переводит коэффициенты с аддитивной шкалы логарифма шансов на мультипликативную шкалу шансов:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

print(exp(coef(fit)))
print(exp(confint(fit)))

exp(0.307) ≈ 1.36, и вот честный шаблон фразы, который стоит запомнить: «каждая дополнительная миля на галлон умножает шансы механической коробки примерно на 1.36». Отношение шансов выше 1 повышает шансы, ниже 1 понижает, ровно 1 означает отсутствие эффекта, — поэтому вердикт доверительного интервала для отношений шансов формулируется как «не содержит ли интервал 1?» (а не нуля; ноль был границей ещё на шкале логарифма шансов).

Следите за формулировками: шансы не являются вероятностями. Шансы = p / (1 − p), поэтому вероятности 0.75 соответствуют шансы 3. Умножение шансов на 1.36 — не то же самое, что умножение вероятности на 1.36, и когда исход частый, разрыв велик. Отношение шансов 2 для редкого исхода ведёт себя как «примерно двойной риск»; для исхода с частотой 50% — категорически нет. Никогда не сообщайте отношение шансов формулировками отношения рисков («в 1.36 раза вероятнее»), если только исход не редкий.

Предсказанные вероятности: подвох type = "response"

Самый частый баг логистической регрессии в дикой природе:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)
new_cars <- data.frame(mpg = c(15, 20, 25, 30))

print(predict(fit, newdata = new_cars))                      # log-odds!
print(predict(fit, newdata = new_cars, type = "response"))   # probabilities

Первый вызов возвращает умолчание type = "link" — предсказания на шкале логарифма шансов, вместе с отрицательными значениями. Второй возвращает настоящие вероятности. Если ваши «вероятности» когда-нибудь выходят отрицательными или больше 1, вот почему. Запустите блок: у машины с расходом 15 миль на галлон практически нет шансов оказаться механической, машина с 30 очень вероятно механическая, а S-образная кривая изгибается посередине.

Классификация: порог и таблица ошибок

Вероятности превращаются в предсказанные классы выбором порога — по умолчанию берут 0.5, — а честной оценочной ведомостью служит таблица «предсказанное против фактического»:

fit <- glm(am ~ mpg, data = mtcars, family = binomial)

probs <- predict(fit, type = "response")   # fitted probabilities, one per car
predicted <- probs > 0.5

print(table(actual = mtcars$am, predicted = predicted))

Ячейки на диагонали — верные решения; две внедиагональные ячейки — две разные ошибки (предсказать механическую для автоматической и наоборот). Одна только общая доля верных ответов может сильно польстить модели: если 95% клиентов не уходят, стратегия «предсказывать, что не уходит никто» даёт 95% при нуле пойманных уходов, — так что всегда смотрите на оба типа ошибок. И 0.5 — соглашение, а не закон: когда две ошибки стоят по-разному, сдвигайте порог соответственно.

Одна оговорка о честности: эта таблица оценивает модель на тех же данных, на которых она подгонялась, что ей льстит. Настоящая оценка откладывает данные, которых модель никогда не видела.

Несколько предикторов

В точности как в lm() — добавляйте члены через +, и каждая интерпретация обзаводится оговоркой «при прочих равных»:

fit2 <- glm(am ~ mpg + qsec, data = mtcars, family = binomial)
summary(fit2)
print(exp(coef(fit2)))

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

Предостережения

• Полное разделение. Если предиктор идеально разделяет исход (каждая машина выше некоторого расхода механическая, каждая ниже — автоматическая), коэффициент максимального правдоподобия хочет стать бесконечным. R предупреждает — glm.fit: fitted probabilities numerically 0 or 1 occurred — и сообщает огромные коэффициенты с абсурдными стандартными ошибками. Не пускайте такие числа в дело: упростите модель, соберите больше данных или используйте штрафной метод (пакеты brglm2 или logistf).
• Достаточно событий. Связывающим ограничением является число более редкого исхода, а не общее число строк. Старое правило большого пальца требует порядка 10–15 событий на предиктор; примеры на 32 машинах здесь нужны для объяснения механики, а не как образец публикуемых размеров выборки.
• Отношения шансов не являются отношениями рисков, когда исход частый, — разобрано выше и повторено потому, что рецензенты это заметят, даже если вы нет.

Что вы уносите с собой

• Бинарный исход → glm(y ~ x, data = df, family = binomial); никогда не забывайте family.
• Сырые коэффициенты — это логарифмы шансов; exp(coef(fit)) даёт отношения шансов, и нулевым значением для их интервалов является 1.
• Шаблон фразы: «увеличение x на единицу умножает шансы исхода на exp(b)».
• predict(..., type = "response") для вероятностей — умолчание возвращает логарифмы шансов, путаница номер один.
• Классифицируйте по порогу и оценивайте по таблице ошибок; одна лишь доля верных ответов может лгать.
• Следите за предупреждениями о разделении, считайте свои события и не наряжайте отношения шансов в отношения рисков.

Дальше: механика, стоящая за каждым интервалом, который вы уже видели, — доверительные интервалы через t.test(), confint() и prop.test().