28 лют. 2013 р.

Снова о логистической регрессии

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.

Читать далее

24 лют. 2013 р.

Небольшой алгоритм для импорта данных из SPSS в R

При работе со сторонними данными, набранными в SPSS, у меня возникала проблема: спсс-овские метки переменных не экспортировались в R, что затрудняло работу с получившимся массивом (переменные в этом случае назывались примерно так v1, v2 etc). Решение нашлось в блоге http://strengejacke.wordpress.com/2013/02/22/migrating-from-spss-to-r-rstats/:

my.data=read.spss('data.sav',to.data.frame=FALSE, use.value.labels=T)

my.data.table=as.data.frame(my.data)

my.data.lab=attr(my.data,'variable.labels')

names(my.data.table)=my.data.lab

Чтобы нарисовать картинку с нужным заглавием, можно воспользоваться созданной вспомогательной переменной my.data.lab (напр., для переменной в 11 столбце):

plot(my.data.table[,11], main=my.data.lab[11])


Читать далее

Кребс и Авсоний

Одна из работ покойного литературо­веда и переводчик­а Михаила ГАСПАРОВА повествует­ о ныне забытом поэте времени упадка Римской империи — АВСОНИИ. Основная мысль: АВСОНИЙ, сам того не подозревая­ (ведь это только мы знаем, что жил он в «эпоху упадка», а для него Рим был жив...), как бы наводил порядок в наследии своих гениальных­ предшестве­нников, готовя его к архивному хранению; он виртуозно пользовалс­я всей уходящей техникой работы со словом, изящно соединяя в образе землю и небо, растения и животных, культуру и природу, каждый раз выстраивая­ мысленную перспектив­у универсума­.

Сегодня мне показалось­, что мой любимец последнего­ времени — Иоганн Людвиг КРЕБС — является таким вот Авсонием эпохи немецкого барокко: он столь же малоизвест­ен, столь же композитор­ски виртуозен,­ он пользуется­ всеми приёмами своих учителей, так же широк дыханием, универсале­н мыслью, он почти ничего не открывает нового, живя в эпоху «упадка контрапунк­тического искусства». Он наводит порядок, чтобы наследие баховского­ искусства передать потомкам.

Ирония судьбы: после КРЕБСА наступает эпоха, возродивша­я АВСОНИЯ — «est in arundineis modulatio musica ripis»...

Читать далее

17 лют. 2013 р.

Шесть маленьких прелюдий

Как снимают клипы на Баха? По моим наблюдениям: либо какая-то статичная фотография, либо музыкант за инструментом. А вот тут милая девушка, играя то, что играют дети в школе, отошла от этого канона и преподнесла свою улыбку — это заразительно и красиво!


Читать далее

14 лют. 2013 р.

Слабая женщина

Сегодня перед манипуляционным кабинетом старушка-далеко-за-70 витийствовала. Надо сказать, что у неё на шее был шов, залитый зелёнкой. Так вот:

«Захожу я, девоньки (от автора — внимательно слушавшие её бабки-ровестницы), в свой подьезд, а из лифта мне кричат: долго ещё ждать? Ну я, дура старая, забежала, нажала свой 8 этаж, повернулась лицом к дверям и еду»…

Прервусь. На этом месте я было понадеялся услышать, что стоявший в лифте за её спиной плюгавенький мужичок лет 30… ну… вы понимаете: стал грязно домогаться. Сюжет оказался занимательно непредсказуем:

…«приехали мы, я вышла, держу в руке кошелёк, и тут он как накинется, приставил нож к горлу, молчи, говорит, сука, а сам другой рукой кошелёк вырывает. Я женщина слабая, прижала кошелёк к груди и как заору! Выскочил мой внук, а он такой — высокий, тот побежал, внук его в два прыжка настигнул — ну прямо тигр, придавил к полу, так и держал, пока милиция приехала, человек двенадцать. Потом нас следователь допрашивал, потом прокурор приехал. 6 часов мы были в участке! А потом мне в челюстно-лицевой порез зашивали — 5 швов наложили! Вот теперь хожу и сдаю кровь на столбняк!»

Слабая женщина! Слона на скаку остановит! И хобот ему оторвёт!

Читать далее