7.5 Модель мультиномиального логита

Логистическая регрессия обычно используется как бинарный классификатор для выборок с альтернативным откликом. Однако этот метод можно обобщить и на случай с несколькими классами. В качестве моделируемого отклика \(Y\) могут использоваться номинальные или порядковые переменные, и в обоих случаях мы сталкиваемся с многомерным биномиальным распределением.

Рассмотрим сначала модель логита (multinomial logit) для отклика с несколькими классами. Пусть \(J\) обозначает число классов для случайной независимой величины \(Y\) (требование независимости всех \(J\) альтернатив тут играет важное значение). Многомерное распределение вероятности наблюдения каждого из ее классов может быть представлено вектором \(\{ \pi_1, \dots \pi_J\}\). Тогда для каждого из \(n\) независимых наблюдений оценка \(\pi_j\) будет соответствовать вероятности того, насколько предпочтительно присвоение произвольному тестируемому объекту метки класса \(j\).

Отметим, что для оценки всех вероятностей достаточно построить \(J - 1\) моделей, поскольку недостающее значение легко получить из условий нормировки (сумма \(\pi_j\) всегда равна 1). Поэтому любой из классов можно принять за базовый (примем для определенности, что это класс с номером 1), и тогда совокупность логит-моделей (baseline-category logit) будет иметь вид:

\[ \log(\pi_j / \pi_1) = \alpha_j + \beta_j \mathbf{x}, \, \, j=2, \dots, J.\]

Коэффициенты системы логит-уравнений оцениваются совместно путем максимизации общего критерия правдоподобия, т.е. другая форма анализа, основанная, например, на трех последовательно построенных моделях регрессии локально для каждого класса, может дать совсем иные результаты (Agresti, 2007).

Для анализа категориальных данных в R существуют мощные специализированные пакеты VGAM и mlogit. Мы же воспользуемся более простой в освоении функцией multinom() из пакета nnet, которая осуществляет оценку коэффициентов системы логит-моделей с использованием алгоритмов построения искусственных нейронных сетей (см. раздел 7.6). Для данных по ирисам Фишера получаем:

library(nnet)
MN.iris <- multinom(Species ~. , data = iris)

## # weights:  18 (10 variable)
## initial  value 164.791843 
## iter  10 value 16.177348
## iter  20 value 7.111438
## iter  30 value 6.182999
## iter  40 value 5.984028
## iter  50 value 5.961278
## iter  60 value 5.954900
## iter  70 value 5.951851
## iter  80 value 5.950343
## iter  90 value 5.949904
## iter 100 value 5.949867
## final  value 5.949867 
## stopped after 100 iterations

summary(MN.iris)

## Call:
## multinom(formula = Species ~ ., data = iris)
## 
## Coefficients:
##            (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor    18.69037    -5.458424   -8.707401     14.24477   -3.097684
## virginica    -23.83628    -7.923634  -15.370769     23.65978   15.135301
## 
## Std. Errors:
##            (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor    34.97116     89.89215    157.0415     60.19170    45.48852
## virginica     35.76649     89.91153    157.1196     60.46753    45.93406
## 
## Residual Deviance: 11.89973 
## AIC: 31.89973

Анализируя коэффициенты, можно заметить, например, что ширина чашелистика Sepal.Width уменьшается при переходе от setosa к virginica, причем логарифм отношения шансов \(\log(\pi_{vir} / \pi_{set})\) уменьшается на величину \(\beta_{23} = -15.3\) при увеличении ширины на 1 см.

Метод fitted(), применяемый к объекту multinom, возвращает матрицу апостериорных вероятностей принадлежности каждого наблюдения к соответствующим классам. Предсказание класса осуществляется с учетом максимума такой вероятности:

Probs <- fitted(MN.iris)

pred <- apply(Probs, 1, function(x) colnames(Probs)[which(x == max(x))])
(table(Факт = iris$Species, Прогноз = pred))

##             Прогноз
## Факт         setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         49         1
##   virginica       0          1        49

Acc <- mean(pred == iris$Species)
paste("Точность=", round(100*Acc, 2), "%", sep = "")

## [1] "Точность=98.67%"

Поскольку функции, выполняющие подгонку моделей мультиноминального логита, обычно не проводят анализа статистической значимости коэффициентов, рассчитаем \(p\)-значения самостоятельно с использованием теста Вальда и \(t\)-статистики:

z <- summary(MN.iris)$coefficients/summary(MN.iris)$standard.errors
# p-значения на основе теста Вальда 
(1 - pnorm(abs(z), 0, 1))*2

##            (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor   0.5930295    0.9515807   0.9557828    0.8129231   0.9457075
## virginica    0.5051288    0.9297757   0.9220685    0.6955898   0.7417773

# p-значения на основе  t-статистики
pt(z, df = nrow(iris) - 5, lower = FALSE) 

##            (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
## versicolor   0.2969240    0.5241679   0.5220705    0.4066286   0.5270993
## virginica    0.7469061    0.5350513   0.5388982    0.3480821   0.3711264

И мы с недоумением видим, что, несмотря на приличную точность предсказаний модели, все отношения \(\log(\pi_k / \pi_1)\) согласно вычисленным \(р\)-значениям статистически не отличаются от 0, т.е. все шансы равны. Ответ видится в резковатой реплике проф. Б. Рипли (B. Ripley) при дискуссии в Интернете: “Попробуйте почитать наши книги. Там объясняется, почему тесты Вальда выполнять нельзя, и что асимптотическая теория может здесь дико вводить в заблуждение”. Читатели, желающие убедиться в правоте сказанного могут пересчитать \(р\)-значения коэффициентов с использованием теста множителей Лагранжа или бутстрепа, в частности, функцией cluster.im.mlogit() из пакета clusterSEs.