По долгу службы мне пришлось разбираться с логистической регрессией в R. Не могу сказать, что я с ней особо хорошо разобрался. но, что выросло — то выросло. Я написал небольшой тьюториал по этой штуке, может вам он тоже будет полезен. Буду рад комментариям, замечаниям, предложениям, ибо статистика наверное все-таки не мой конек.

Логистическая регрессия — регрессия для которой зависимой переменной является логарифм шансов получить то или иное значение бинарной переменной. В моем случае, меня интересовала связь между тремя предикторами и одной зависимой переменной. Зависимой переменной была вероятность выбора стимула в задаче вынужденного выбора (chosen), предикторами — частота предъявления этого стимула на предыдущем этапе (target_freq, 1/5), тип выбора (task_type, выбор А/выбор B, суть неважна),  и сложность стимула (complexity, метрическая переменная, 5..30).

Итак, для начала смотрим на данные:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
> summary.formula(as.numeric(chosen)~cut2(complexity,g=4)+task_type,data=mydata,method="cross")
 
 mean by cut2(complexity, g = 4), task_type 
 
+------------------+
|N                 |
|as.numeric(chosen)|
+------------------+
+-----------------------+--------+--------+--------+
|cut2(complexity, g = 4)| Type A | Type B |   ALL  |
+-----------------------+--------+--------+--------+
|        [ 5,13)        |1064    |1048    |2112    |
|                       |1.554511|1.472328|1.513731|
+-----------------------+--------+--------+--------+
|        [13,16)        | 749    | 727    |1476    |
|                       |1.524700|1.524072|1.524390|
+-----------------------+--------+--------+--------+
|        [16,21)        |1001    |1021    |2022    |
|                       |1.469530|1.514202|1.492087|
+-----------------------+--------+--------+--------+
|        [21,30]        | 746    | 712    |1458    |
|                       |1.379357|1.453652|1.415638|
+-----------------------+--------+--------+--------+
|        ALL            |3560    |3508    |7068    |
|                       |1.487640|1.491448|1.489530|
+-----------------------+--------+--------+--------+
 
> xtabs(~ chosen+task_type, mydata)
      task_type
chosen Type A Type B
     0   1824   1784
     1   1736   1724
 
> xtabs(~ chosen+target_freq, mydata)
      target_freq
chosen    1    5
     0 1803 1805
     1 1727 1733

chosen — фактор с двумя уровнями (не выбран, выбран), при использовании as.numeric он принимает значения 1, 2. Соответственно, 1.55 в верхней таблице означает 55% вероятность выбора стимула. Единственное, что заметно — снижение вероятности выбора в зависимости от сложности стимула для задач типа A. Особо отсутствия данных нигде вроде нет, действуем дальше.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
> mylogit<- glm(chosen~task_type*complexity*target_freq, mydata, family=binomial(link="logit"), na.action=na.pass)
> summary(mylogit)
 
Call:
glm(formula = chosen ~ task_type * complexity * target_freq,
    family = binomial(link = "logit"), data = mydata, na.action = na.pass)
 
Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.4602  -1.1618  -0.9218   1.1887   1.5665  
 
Coefficients:
                                         Estimate Std. Error z value Pr(&gt;|z|)
(Intercept)                              0.630081   0.150501   4.187 2.83e-05 ***
task_typeType B                         -0.557237   0.210878  -2.642  0.00823 **
complexity                              -0.043474   0.008897  -4.887 1.03e-06 ***
target_freq5                             0.318621   0.210413   1.514  0.12996
task_typeType B:complexity               0.037877   0.012443   3.044  0.00234 **
task_typeType B:target_freq5            -0.522420   0.296817  -1.760  0.07840 .
complexity:target_freq5                 -0.017481   0.012434  -1.406  0.15975
task_typeType B:complexity:target_freq5  0.028035   0.017514   1.601  0.10944
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1 
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 9795.2  on 7067  degrees of freedom
Residual deviance: 9718.5  on 7060  degrees of freedom
AIC: 9734.5
 
Number of Fisher Scoring iterations: 4

R мило сообщает нам:

1) какую модель мы проверяем

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

3) собственно регрессионные коэффициенты, их стандартные ошибки, z-статистику и уровни значимости

4) отклонение для нулевой модели (модели без предикторов) и остаточное отклонение (я не совсем уверен, что на русском это будет называться именно так) и индекс пригодности модели (AIC)

Из всего этого нас больше всего интересуют собственно значения коэффициентов, их уровни значимости и то, что было в пункте 4. Коэффициенты можно интерпретировать следующим образом. При изменении значения предиктора на единицу, значение логарифма отношения шансов зависимой переменной изменится на величину коэффициента. Т.е. при увеличении сложности на 1, значение логарифма отношения шансов снизится на 0.04. Это не совсем удобно, поэтому имеет смысл перевести логарифмы в собственно отношение шансов. Делается это так:

1
2
3
4
5
> exp(mylogit$coefficients)
 (Intercept)                         task_typeType B                              complexity                            target_freq5
 1.8777618                               0.5727895                               0.9574576                               1.3752295
 task_typeType B:complexity            task_typeType B:target_freq5                 complexity:target_freq5 task_typeType B:complexity:target_freq5
 1.0386030                               0.5930838                               0.9826705                               1.0284316

Вполне логично, взяли экспоненты от логарифма, получили значение шансов. Теперь можно сказать, что при увеличении complexity на 1, отношение шансов того, что наш стимул будет выбран, увеличится в 0.95 раз, т.е. по сути снизится. Т.е. было, допустим 3 к 2 (1.5), станет 1.425.  В общем. тоже не шибко информативно, но уже лучше.

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

1
2
> mylogit$null.deviance - mylogit$deviance
[1] 76.77171

Количество степеней свободы:

1
2
> mylogit$df.null - mylogit$df.residual
[1] 7

И уровень значимости:

1
2
> dchisq(mylogit$null.deviance-mylogit$deviance, mylogit$df.null-mylogit$df.residual)
[1] 2.931273e-14

Таким образом можно сказать, что наша модель значимо лучше чем нулевая модель.

Взаимодействие переменных target_freq и task_type проще всего посмотреть так:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
> summary.formula(as.numeric(chosen)~target_freq+task_type,data=mydata,method="cross")
 
 mean by target_freq, task_type 
 
+------------------+
|N                 |
|as.numeric(chosen)|
+------------------+
+-----------+--------+--------+--------+
|target_freq| Type A | Type B |   ALL  |
+-----------+--------+--------+--------+
|    1      |1775    |1755    |3530    |
|           |1.482817|1.495726|1.489235|
+-----------+--------+--------+--------+
|    5      |1785    |1753    |3538    |
|           |1.492437|1.487165|1.489825|
+-----------+--------+--------+--------+
|    ALL    |3560    |3508    |7068    |
|           |1.487640|1.491448|1.489530|
+-----------+--------+--------+--------+

Данные glm говорят нам о том, что при task_type=Type B стимулы, вероятность выбора снижается, если стимулы предъявлены 5 раз, что мы и видим в этой таблице. Но различия незначительны и значимость лишь на уровне тенденции, так что я сочту это артефактом.

Чтобы понять, что все-таки происходит при изменении сложности стимула, полезно посмотреть на графики. Для этого используем пакет ggplot2.

1
2
> qplot(complexity, as.numeric(chosen)-1, data=mydata, facets=task_type~.,fun.y="mean",stat="summary")+
  stat_smooth(method="glm",family="binomial", formula = y ~ ns(x, 2)) + theme_bw()

Первая часть команды говорит, что нас интересует график средних (fun.y=»mean», stat=»summary»), где переменная x — complexity, y — численное значение (chosen)-1, при этом график должен быть разбит на уровни по task_type. Вторая часть — stat_smooth — говорит, что нам нужна линия аппроксимации для биномиальной модели с использованием формулы (y~ns(x,2)). Эту формулу я взял из описания stat_smooth, насколько понимаю, она делает аппроксимацию нелинейной. Наконец, theme_bw() слегка меняет цветовую палитру графика.

Вот, что получается:

Графики логистической регрессии в зависимости от типа задач (щелкните для увеличения)

В одной задаче имеет место почти линейное убывание вероятности выбора в зависимости от сложности, в другой зависимости нет.

Примерно вот так все, буду рад комментариям, замечаниям, предложениям.