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, подключаемых командами:
library(aod)
library(ggplot2)
Пример. Исследователь ищет, как переменные GRE (Graduate Record Exam scores — оценки во время обучения в вузе), GPA (grade point average — средний балл) и престиж вуза влияют на поступление в аспирантуру. Т. о. искомая переменная, поступил или не поступил (admit/don’t admit), является бинарной.
В описанном ниже анализе мы воспользовались сгенерированными гипотетическими данными, которые можно получить, используя
R:
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
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. 38-40).
Дискриминантный анализ двух групп. Многомерный метод для бинарной зависимой переменной.
Hotelling’s T2. Зависимая переменная 0/1 преобразуется в группирующую переменную, а предикторы — в зависимые переменные. Это позволяет выполнить общий тест значимости, но не даёт индивидуальных коэффициентов для каждой переменной, соответственно остаётся неясным вклад каждого такого «предиктора» в эффект прочих «предикторов».
Использование логистической модели
Код, приведённый ниже, выполняет логистическое моделирование с использованием функции glm() function, но сначала мы превратим переменную rank в фактор, чтобы показать, что она должна трактоваться как категориальная переменная:
mydata$rank <- factor(mydata$rank)
mylogit <- glm(admit ~ gre + gpa + rank, data = mydata, family = "binomial")
Поскольку мы дали нашей модели имя mylogit, R не выдал результаты регресии на экран. Это, однако, легко сделать командой summary():
summary(mylogit)
## 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(), чтобы получить доверительные интервалы рассчитанных коэффициентов.
confint(mylogit)
## 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 (т. е., 4-е и 5-е слагаемое в модели). Для контраста мы умножаем один из них на 1, а другой на −1. Прочие слагаемые модели не участвуют в тесте, поэтому мы умножаем их на 0. Вторая линия кода использует равенство L=l для того, чтобы сказать R, что мы решили выполнить тест с использованием вектора l.
l <- cbind(0, 0, 0, 1, -1, 0)
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():
exp(coef(mylogit))
## (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 на уровне их средних. Создадим и посмотрим на таблицу данных:
newdata1 <- with(mydata, data.frame(gre = mean(gre), gpa = mean(gpa), rank = factor(1:4)))
newdata1
## 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.
newdata1$rankP <- predict(mylogit, newdata = newdata1, type = "response")
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 <- with(mydata, data.frame(gre = rep(seq(from = 200, to = 800, length.out = 100),
4), gpa = mean(gpa), rank = factor(rep(1:4, each = 100))))
Код, генерирующий предсказанные вероятности, аналогичен приведённому выше.
newdata2 <- cbind(newdata2, predict(mylogit, newdata = newdata2, type = "link",
se = TRUE))
newdata2 <- within(newdata2, {
PredictedProb <- plogis(fit)
LL <- plogis(fit - (1.96 * se.fit))
UL <- plogis(fit + (1.96 * se.fit))
})
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)
## [1] 41.46
А разница степеней свободы может быть получена так:
with(mylogit, df.null - df.residual)
## [1] 5
Наконец, p-значение вычисляется так:
with(mylogit, pchisq(null.deviance - deviance, df.null - df.residual, lower.tail = FALSE))
## [1] 7.578e-08
Хи-квадрат 41.46 с 5 степенями свободы и ассоциированным p-значением, меньшим 0.001, показывают нам, что наша модель, в целом, описывает данные значительно лучше, чем нулевая модель. Этот тест иногда называют «likelihood ratio test».
Примечания
Hosmer, D. & Lemeshow, S. (2000). Applied Logistic Regression (Second Edition). New York: John Wiley & Sons, Inc.
Long, J. Scott (1997). Regression Models for Categorical and Limited Dependent Variables. Thousand Oaks, CA: Sage Publications.
Читать далее
