Данный ноутбук является конспектом по курсу «R для теории вероятностей и математической статистики» (РАНХиГС, 2018). Автор ноутбука вот этот парень по имени Филипп. Если у вас для него есть деньги, слава или женщины, он от этого всего не откажется. Ноутбук распространяется на условиях лицензии Creative Commons Attribution-Share Alike 4.0. При использовании обязательно упоминание автора курса и аффилиации. При наличии технической возможности необходимо также указать активную гиперссылку на страницу курса. На ней можно найти другие материалы. Фрагменты кода, включенные в этот notebook, публикуются как общественное достояние.
В прошлый раз мы попытались просечь метод максимального правдоподобия и оценили парочку простых моделей. Сегодня мы продолжим разбираться с ним и оценим свою первую настоящую модель. Она называется логистической регрессией и очень часто, с небольшими видоизменениями, используется на практике.
Что мы делали до этого? Мы:
В этом подходе мы забывали про очень существенную вещь: наблюдаемая нами случайная величина зависит от других случайных величин. Например, цена на квартиру, зависит от её метража, этажности здания, криминагенности района и многих других факторов. Выживет ли человек на Титанике зависит от его пола, возраста, каюты и тп. Так почему бы не учесть все эти факторы при моделировании.
В случае квартир мы должны спрогнозировать цену. Это непрерывная переменная. Задача, состоящая в строительстве модели для непрерывной переменной, называется задачей регрессии. В случае Титаника мы должны ответить на вопрос "Умрёт ли человек?" либо да, либо нет, а также найти вероятность того, что да. Такая задача, с конечным множеством ответов, называется задачей классификации. Именно задачу классификации мы с вами будем решать в этой тетрадке.
У нас в руках есть данные о пассажирах титаника с соревнования на Kaggle (Kaggle -это такая платформа для проведения соревнований по анализу данных. Там много чего интересного есть). Наша цель: научиться предсказывать умрёт ли пассажир титаника по его характеристикам. Попутно неплохо было бы научиться предсказывать вероятность смерти.
Подгружаем все необходимые пакеты. Если каких-то из них у нас нет, устанавливаем.
library('dplyr')
library('ggplot2')
library('maxLik') # пакет для метода макс. правдоподобия
options(repr.plot.width=4, repr.plot.height=3)
Выборку пассажиров титаника надо будет разбить на две части: тренировочную и тестовую. Обычно, на тренировочной выборке учат модель, на тестовой выборке проверяют насколько хорошо она работает. Выборка, на которой мы учимся, не подходит для этих целей, так как знания об этих объектах уже будут лежать в нашей модели, а в реальной жизни будут встречаться и новые, ранее неизвестные объекты. Часть, которую мы отложили, алгоритм ещё не видел. Значит можно использовать её для имитации этих самых новых объектов из реальной жизни и посмотреть какое на них будет качество.
Такой подход помогает избежать переобучения модели. Переобучением называют ситуации, когда модели тупо запоминают тренировочную выборку вместо того, чтобы извлечь из неё какие-то закономерности. Обычно на обучение пускают $70-80\%$ выборки, а на тестирование $30-20\%$.
df = read.csv('/Users/fulyankin/Yandex.Disk.localized/R/R_prob_data/titanic_train.csv', sep=',')
# Можно посмотреть на первые пять строк!
head(df, 5)
PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked |
---|---|---|---|---|---|---|---|---|---|---|---|
1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22 | 1 | 0 | A/5 21171 | 7.2500 | S | |
2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Thayer) | female | 38 | 1 | 0 | PC 17599 | 71.2833 | C85 | C |
3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26 | 0 | 0 | STON/O2. 3101282 | 7.9250 | S | |
4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35 | 1 | 0 | 113803 | 53.1000 | C123 | S |
5 | 0 | 3 | Allen, Mr. William Henry | male | 35 | 0 | 0 | 373450 | 8.0500 | S |
# Можно посмотреть на последние пять строк!
tail(df, 5)
PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
887 | 887 | 0 | 2 | Montvila, Rev. Juozas | male | 27 | 0 | 0 | 211536 | 13.00 | S | |
888 | 888 | 1 | 1 | Graham, Miss. Margaret Edith | female | 19 | 0 | 0 | 112053 | 30.00 | B42 | S |
889 | 889 | 0 | 3 | Johnston, Miss. Catherine Helen "Carrie" | female | NA | 1 | 2 | W./C. 6607 | 23.45 | S | |
890 | 890 | 1 | 1 | Behr, Mr. Karl Howell | male | 26 | 0 | 0 | 111369 | 30.00 | C148 | C |
891 | 891 | 0 | 3 | Dooley, Mr. Patrick | male | 32 | 0 | 0 | 370376 | 7.75 | Q |
# Можно кинуть свой взор на таблицу и посмотреть какой тип у каких переменных
glimpse(df)
Observations: 891 Variables: 12 $ PassengerId <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, ... $ Survived <int> 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0... $ Pclass <int> 3, 1, 3, 1, 3, 3, 1, 3, 3, 2, 3, 1, 3, 3, 3, 2, 3, 2, 3... $ Name <fctr> Braund, Mr. Owen Harris, Cumings, Mrs. John Bradley (F... $ Sex <fctr> male, female, female, female, male, male, male, male, ... $ Age <dbl> 22, 38, 26, 35, 35, NA, 54, 2, 27, 14, 4, 58, 20, 39, 1... $ SibSp <int> 1, 1, 0, 1, 0, 0, 0, 3, 0, 1, 1, 0, 0, 1, 0, 0, 4, 0, 1... $ Parch <int> 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 1, 0, 0, 5, 0, 0, 1, 0, 0... $ Ticket <fctr> A/5 21171, PC 17599, STON/O2. 3101282, 113803, 373450,... $ Fare <dbl> 7.2500, 71.2833, 7.9250, 53.1000, 8.0500, 8.4583, 51.86... $ Cabin <fctr> , C85, , C123, , , E46, , , , G6, C103, , , , , , , , ... $ Embarked <fctr> S, C, S, S, S, Q, S, S, S, C, S, S, S, S, S, S, Q, S, ...
Давайте наконец выясним, что какая переменная обозначает.
Часть из этих переменных являются действительными (возраст, цена билета), часть категориальными (класс, в котором плыл человек, пол человека). В выдаче от glimpse
числовые переменные обозначены как <int>
и <dbl>
. Для категориальных переменных у R есть свой специальный тип, <fctr>
, факторный. Это бывает очень удобно при оценке моделей.
Дело в том, что метод максимального правдоподобия (как в прочем и любые другие методы и модели) понимает только цифры. Напрмер, он не знает что такое male
и female
и не может посчитать от них среднее. Приходится придумывать способы перегнать такие переменные в цифры. В случае пола, например, можно обозначить всех мужчин единичками, а всех женщин нулями. Тогда итоговая модель, получая единичку, будет понимать, что c ней делать, а мы будем понимать, что она делает это с мужчинами. Такая переменная, принимающая два значения: ноль и один, называется дамми-переменной.
В случае более сложных переменных, например порта посадки, возникает желание переобозначить Шербур за единички, Квинстаун за двойки, а Саутгемптон за тройку. Это неправильный подход. Если так сделать, то алгоритмы будут дамать, что третий порт самый лучший, а первый самый плохой. Очень большой вопрос почему мы определили порядок портов именно так. Это запутает модель. Чтобы не путать алгоритм, поступают следующим образом: делают из категориальной переменной много дамми. Вот так:
Если человек плыл из порта $C$, то у него в соответствующем столбце будет стоять единичка, а в остальных нули. По аналогии получится с другими портами.
Приятный бонус: Не нужно бороться с пропусками в данных. Для пропусков будет своя дамми и мы не потеряем информации. Такой способ работы с категориями называется One Hot Encoding. Дословно можно перевести как одно горячее кодирование.
В чём минус одного горячего кодирования? Можно обжечься. Если категорий внутри переменной много, то мы переработаем её в слишком большое число столбцов. В конечном итоге это может привести к тому, что модели не хватит наблюдений, чтобы научиться в этом огромном числе столбцов ориентироваться.
Есть и другие способы обработки категориальных переменных, но мы будем пользоваться этим. Он самый простой. Более того, тип переменной <fctr>
позволяет делать это по дефолту без предобработки. Обратите внимание, что часть категориальных переменных, R воспринял как числовые. Это нужно исправить.
# команда mutate делает из выборки мутанта, изменяя её :)
df <- mutate(df, Pclass = as.factor(Pclass), Survived = as.factor(Survived))
Давайте сразу же разберёмся со всякими техническими тонкостями, связанными с факторными переменными. Этот тип данных нужен, чтобы было удобнее работать с пакетами, которые оценивают модели, а также чтобы было удобно строить рисунки. При этом, при ручной работе с факторными переменными, будут возникать проблемы.
mean(df_train$Survived)
Warning message in mean.default(df_train$Survived): “argument is not numeric or logical: returning NA”
Посмотрим на структуру факторной переменной. Это помогает сделать команда str
. В каждой факторной переменой записаны две штуки. Первая - её значения, в данном случае это ноль и один. Обратите внимание, они записаны как строки, а не как числа. Вторая - уровни, которые этим значениям приписываются. R отталкиваясь от этих уровней сам делает внутри пакетов для оценивания моделей одно горяее кодирование.
str(df_train$Survived)
Factor w/ 2 levels "0","1": 1 2 2 1 2 1 1 2 1 2 ...
Чтобы работать с переменной, как с числом, придётся преобразовать её обратно. Если мы сделаем это командой as.numeric
от столбца, на метсо нулей и единиц подставятся номера категорий (факторов). В данном случае, либо первая, либо вторая. Это плохо, потому что изначально категории были закодированы нулями и единицами.
str(as.numeric(df_train$Survived))
num [1:713] 1 2 2 1 2 1 1 2 1 2 ...
Чтобы избежать этого, придётся сделать переменную сначала тектсовой, потом уже числовой.
str(as.numeric(as.character(df_train$Survived)))
num [1:713] 0 1 1 0 1 0 0 1 0 1 ...
Это была малоприятная техническая деталь, которая обычно взрывает мозг многим ребятам, которые впервые работают в R с таблицами. Поэтому я решил это тут показать.
На самом деле, на этом предобработка переменных не заканчивается, есть огромное количество тонкостей. Большую часть времени аналитики проводят, подготавливая данные. Меньшую часть, оценивая модели. Всем этим тонкостям посвящены курсы по машинному обучению и эконометрике. Мы же всё ещё в курсе матстата, поэтому в этом блокноте сосредоточимся на методе максимального правдоподобия.
# Выбросим неинтересные для нас переменные и посмотрим по интересным основную статистику
variables <- c('Survived', 'Pclass', 'Sex', 'Age', 'SibSp', 'Parch', 'Fare')
df = df[variables]
summary(df[c( 'Age', 'SibSp', 'Parch', 'Fare')])
Age SibSp Parch Fare Min. : 0.42 Min. :0.000 Min. :0.0000 Min. : 0.00 1st Qu.:20.12 1st Qu.:0.000 1st Qu.:0.0000 1st Qu.: 7.91 Median :28.00 Median :0.000 Median :0.0000 Median : 14.45 Mean :29.70 Mean :0.523 Mean :0.3816 Mean : 32.20 3rd Qu.:38.00 3rd Qu.:1.000 3rd Qu.:0.0000 3rd Qu.: 31.00 Max. :80.00 Max. :8.000 Max. :6.0000 Max. :512.33 NA's :177
summary(df[c('Survived', 'Pclass', 'Sex')])
Survived Pclass Sex 0:549 1:216 female:314 1:342 2:184 male :577 3:491
Последняя проблема предобработки, которую мы сегодня исправим заключается в пропусках. Обратите внимание на переменную Age
. В ней $177$ пропусков. Метод максимального правдоподобия оценивает модели на числах и не любит пропуски. Поэтому с ними надо что-то сделать.
Первое решение: выбросить все эти наблюдения. Проблема такого подхода в том, что другие то переменные для этих наблюдений не пропущены. Значит, выбросив строчки, соответствующие пропускам, мы потеряем кучу информации. Полезной информации. Информацию и без всех наших преобразований легко потерять, сложно найти и невозможно забыть, а мы её ещё и выкинуть хотим.
Обычно, если пропусков в переменной не очень много, их заполняют чем-нибудь нейтральным. Например, медианами. Когда алгоритм будет натыкаться на них, он будет думать: "Хммм, ну тут ничего особенного, всё как в среднем по выборке, пойду посмотрю другую переменную". И тогда всё будет хорошо.
Тем не менее, если пропусков в переменной очень много, то лучше выбросить её. Заполнение будет только путать алгоритм. Что лучше, выкинуть строки или заполнить их, во многом зависит от задачи, которую вы решаете. Можно построить несколько разных моделей с разной предобработкой пропусков и посмотреть у какой из них качество будет лучше.
Заполним пропуски в возрасте медианами.
median(df$Age) # не расчитало медиану из-за пропусков :(
median(df$Age, na.rm = TRUE) # включили опцию игнорить пропуски и медиана расчиталась
# сделали срез по всем ячейкам с пропусками и залили их медианой
df[is.na(df$Age),]$Age = median(df$Age, na.rm = TRUE)
median(df$Age) # теперь пропусков нет и всё считается!
Повторюсь, что проблемы, связанные с предобработкой, на этом не заканчиваются. Большую часть времени приходится возиться с данными: работать с пропусками, придумывать новые переменные, обрабатывать неструктурированные переменные (например, текст) и т.п.
Мы останавливаемся на этом только потому, что хотим посмотреть на то, как работает метод максимального правдоподобия на реальной задаче, а не прослушать вводный курс машинки. Разобьём выборку, как и планировали в самом начале, на тестовую и тренировочную.
RowIndex <- sample(1:nrow(df), round(0.8*nrow(df))) # отбираем рандомные индексы для теста и трэйна
RowIndex[1:10]
df_train <- df[RowIndex, ] # тренировчная будет состоять из 80% наблюдений
df_test <- df[-RowIndex, ] # тестовая будет состоять из 20%
dim(df_train)
dim(df_test)
Давайте построим пару визуализаций датасета. Просто, чтобы было. Во-первых давайте посмотрим на распределение выживших и мёртвых. Классы несбалансированы. Примеров с выжившими меньше, чем с мёртвыми. Запомните эту деталь. Она нам понадоибится в самом конце.
ggplot(df_train, aes(x=Survived))+
geom_bar(stat="count", width=0.7, fill="steelblue")
Можно посмотреть на кучу других картинок. Например, вот на такую:
#Построим плотности выживших и невыживших на одном графике.
qplot(data=df_train, x=Age, y=..count.., fill=Survived, geom="density")
# количество закраска
# наблюдений в соотв.
# комп сам выжил или нет
На картинке построены две условные плотности распределения возраста. Одна для выживших, другая для мёртвых. На них видно, что большая часть пожилых погибла. Большая часть младенцев выжила.
qplot(data=df_train,x=Survived,y=Age,geom="violin")
#если поделить пассажиров на выживших и не выживших
#и отдожить по вертикали количественную переменную
#а именно возраст, то мы видим что возраст меняется от
# 0 до 80 и среди тех и других по возрастам распределение
# примерно одинаковое.
Ещё больше картинок можно найти в статье по титанику на хабре.
На первои шаге забудем про все переменные, кроме выживаемости. Вытащим из таблицы только колонку из нулей и единиц, где $1$ означает, что человек выжил, а $0$, что нет. По этой колонке нам хотелось бы оценить вероятность выживания. До сегодняшнего дня мы делали это следущим способом:
$$ L = p \cdot p \cdot \ldots \cdot p \cdot (1-p) \cdot (1-p) \cdot \ldots \cdot (1-p) = p^{\sum_{i=1}^n y_i} \cdot (1-p)^{n- \sum_{i=1}^n y_i}. $$Прологорифмировав, мы получали функцию для оптимизации:
$$ \ln L = \sum_{i=1}^n y_i \ln p + (n- \sum_{i=1}^n y_i) (1-p). $$Пришло время вбить функцию логарифма правдоподобия и посмотреть на пару её значений в разных точках. Мы будем делать это руками. Поэтому давайте вернём переменной Survived
числовой тип.
# вернём пока что числовой тип целевой переменной
df_train = mutate(df_train, Survived = as.numeric(as.character(Survived)))
y_sur = df_train$Survived
lnL <- function(p,y){
ans <- sum(y)*log(p) + (length(y) - sum(y))*(1-p)
return(ans)
}
lnL(0.5, c(1,0,1,0,1))
lnL(0.5, y_sur)
Оптимизируем!
res <- maxLik(lnL, start=0.5, y = y_sur) # какую стартовую точку взять?
summary(res) # 0.5 это логично для оценки вероятности
-------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 4 iterations Return code 1: gradient close to zero Log-Likelihood: 34.16908 1 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] 0.63158 0.03802 16.61 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 --------------------------------------------
mean(y_sur) # всё ок! :) Получили ожидаемый результат.
Для того, чтобы в дальнейшем нам было удобно сопоставлять разные формулы, давайте перепишем уже полученное выше правдоподобие немного иначе и прологарифмируем
$$ L = \prod_{i=1}^{n} p^{y_i} (1-p)^{1-y_i} $$$$ \ln L = \sum_{i=1}^n [y_i \ln p + (1-y_i) \ln(1-p)]. $$Итак, что только что произошло? Мы получили вероятность того, что человек выживет на титанике. Она для каждого из нас равна $0.38$. На интуитивном уровне сразу же понятно, что это очень плохая оценка. Она никак не учитывает в себе, что люди бывают разными.
На картинках выше, мы с вами увидели, что для женщин вероятность выжить намного больше, чем для мужчин. Было бы неплохо оценить для мужчин свою вероятность, а для женщин свою. Мы, зная, что $\hat p_{ML} = \bar x$, легко можем сделать это. Заодно изучим прикольный приём под названием трубочка. Если мы хотим взять синус от числа, а потом прологорифмировать, мы можем сделать это так:
log(sin(8))
Если мы хотим от этого взять ещё косинус, то слева нужно дописать новую функцию.
cos(log(sin(8)))
Постоянно надо будет плодить новые скобки и ... это не удобно. Легко забыть или запутаться. Особенно когда ты работаешь с таблицами. Для того, чтобы не путаться, внутри пакета dplyr
лежит трубочка: удобная опирация для последовательных дйствий.
8 %>% sin %>% log %>% cos
Это выглядит так как будто мы берём число $5$ и пропускаем его через трубочку по очереди в каждую функцию. Никаких скобок и удобно. Обычно трубочками работают с таблицами. Например, ниже мы отфильтровали в таблице все строки с мужчинами, а затем посчитали среднее по колонке Survived
.
df_train %>% filter(Sex == 'male') %>% summarise(male_p = mean(Survived))
df_train %>% filter(Sex == 'female') %>% summarise(female_p = mean(Survived))
male_p |
---|
0.1906874 |
female_p |
---|
0.7251908 |
Такая оценка будет намного лучше нашей первой. Если вы женщина, вероятность выжить для вас довольно высока. Но стоит ли останавливаться на этом? Нет! Идея требует своего развития. Кажется, что вероятность выжить это функция от всех параметров, которыми обладает человек
$$p = F(age, sex, ...).$$Предположим, что она линейная, то есть каждый фактор, пол, возраст и тп, взвешиваются с какими-то весами, $\beta$, и получается вероятность
$$p = \beta_0 + \beta_1 sex + \beta_2 age + \ldots.$$Например, как мы уже выяснили, женщины выживают чаще мужчин, более старые люди выживают редко. Возможно, что $p = 0.2*sex - 0.01 age_i$. Тогда получится, что если переменная $sex$ принимает значение $1$, вероятность выжить подскакивает на $0.2$. При этом, с каждым прожитым годом, вероятность выжить падает. Осталось только понять откуда взять эти две бэты.
Найти методом максимального правдоподобия! Берём и заменям в нашей изначальной формуле все $p$ на $F(x)$ и занимаемся её оптимизацией по параметрам $\beta$
$$ \ln L = \sum_{i=1}^n [y_i \ln F(x_i) + (1-y_i) \ln(1-F(x_i))]. $$Попробуем для начала сделать это только с полом:
$$ \ln L = \sum_{i=1}^n [y_i \ln (\beta_0 + \beta_1 sex_i) + (1-y_i) \ln(1-(\beta_0 + \beta_1 sex_i))]. $$В каждом слагаемом находится своё уникальное значение $sex_i$. Из-за этого у нас не получится так просто вбить функцию правдоподобия, как мы сделали это выше. Хорошо, что у нас есть матрицы и линал, которые помогут :) Предположим, что у нас есть вот такая вот табличка из четырёх наблюдений:
$$ X = \begin{pmatrix} 1 & sex_1 \\ 1 & sex_2 \\ 1 & sex_3 \\ 1 & sex_4 \end{pmatrix} $$Во втором столбике записано то, какого пола каждый человек. В первом столбике написаны единички. Это небольшой технический приём, чтобы получить потом $\beta_0$. Запишем вектор из $\beta$ в виде столбика:
$$ \beta = \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix} $$Тогда если мы перемножим матрицу $X$ и столбик $\beta$, то мы получим
$$ X \cdot \beta = \begin{pmatrix} \beta_0 + \beta_1 sex_1 \\ \beta_0 + \beta_1 sex_2 \\ \beta_0 + \beta_1 sex_3 \\\beta_0 + \beta_1 sex_4 \end{pmatrix} $$Домножим каждой слагаемое на свой $y_i$, тогда:
$$ y*X \cdot \beta = \begin{pmatrix} y_1 \cdot(\beta_0 + \beta_1 sex_1) \\ y_2 \cdot(\beta_0 + \beta_1 sex_2) \\ y_3 \cdot(\beta_0 + \beta_1 sex_3) \\ y_4 \cdot(\beta_0 + \beta_1 sex_4) \end{pmatrix} $$Если мы просуммируем всё это, то получим первый кусочек суммы. Аналогично можно получить второй кусочек. Давайте вобьём это в R.
X <- cbind(1,runif(5)) # cbind объединяет столбики
X # матрица из наблюдений
1 | 0.3590383327 |
1 | 0.7923830762 |
1 | 0.4680810818 |
1 | 0.0004320503 |
1 | 0.0935492811 |
b <- c(2, 4)
b # вектор из параметров
X%*%b # матричное умножение
3.436153 |
5.169532 |
3.872324 |
2.001728 |
2.374197 |
y <- c(0,1,0,1,1)
y*log(X%*%b) # домножаем поэлементно на соответствующие игреки и получаем первое слагаемое
0.0000000 |
1.6427822 |
0.0000000 |
0.6940109 |
0.8646593 |
(1-y)*log(1 - X%*%b) # при поиске второго слагаемого возникли проблемы :(
Warning message in log(1 - X %*% b): “созданы NaN”
NaN |
NaN |
NaN |
NaN |
NaN |
# Откуда проблемы? Отрицательные числа под логарифмами.
1 - X%*%b
-2.436153 |
-4.169532 |
-2.872324 |
-1.001728 |
-1.374197 |
# Итоговое правдоподобие выглядело бы вот так:
lnL <- function(beta,y,X){
ans <- sum(y*log(X%*%b) + (1-y)*log(1 - X%*%b))
return(ans)
}
# Можно даже подбрать ситуации, когда с ним всё окей и нет ошибок:
y <- c(0,1)
X <- cbind(1,c(1/16,1/16))
b <- c(0.5,0.5)
lnL(b,y,X)
Когда мы предположили, что $p = \beta_0 + \beta_1 sex_i$, мы оказались неправы. Наш главный косяк был в том, что вероятность должна принимать значения от нуля до единицы. При предположенной нами зависимости, возникают ситуации, когда она пробивает этот диапазон. Из-за этого под логарифм попадают отрицательные числа и вся концепция рушится. Нужно исправить это. Надо перевести нашу линейную комбинацию в шкалу от нуля до единицы. Для этого нужна вот такая функция:
Она будет отображать действительные числа на отрезок от нуля до единицы. Обычно для этого используют функцию под названием сигмоида:
$$ F(z) = \frac{e^z}{1 + e^z} = \frac{1}{1 + e^{-z}}. $$Почему именно её? Потому что она няшная! Внимательный студент узнал в ней функцию распределения для логистической случайной величины, которая фигугрировала в первой домашке. При её использовании модель получается интерпретируемой. Но об этом вы узнаете позже. Сама по себе модель с такой функцией называется логистической регрессией. Вероятность того, что $y_i = 1$ выпишится как:
$$ P(y_i = 1) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 sex_i)}} $$А логарифм правдоподобия примет вид:
$$ \ln L = \sum_{i=1}^n y_i \ln \left( \frac{1}{1 + e^{-(\beta_0 + \beta_1 sex_i)}}\right) + (1-y_i) \ln\left(1- \frac{1}{1 + e^{-(\beta_0 + \beta_1 sex_i)}}\right). $$Осталось только закодить.
# Это сигмоид!
sigmoid <- function(z){
return(exp(z)/(1 + exp(z)))
}
sigmoid(5)
# это логарифм правдоподобия!
lnL <- function(b, y, X){
answ = sum(y*log(sigmoid(-1*X%*%b)) + (1-y)*log(1 - sigmoid(-1*X%*%b)))
return(answ)
}
y <- c(0,1)
X <- cbind(1,c(1/16,1/16))
b <- c(0.5,0.5)
lnL(b,y,X)
# Пол нужно сделать числовым вида 1 - мужчины, 0 - женщины
X = cbind(1, as.numeric(df_train$Sex) - 1)
lnL(b,y_sur,X)
res <- maxLik(lnL, start=c(0.5,0.5), y = y_sur, X = X)
summary(res)
-------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 4 iterations Return code 1: gradient close to zero Log-Likelihood: -373.7871 2 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] -0.9704 0.1384 -7.013 2.33e-12 *** [2,] 2.4159 0.1831 13.194 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 --------------------------------------------
Всё работает. Теперь вероятность выживания это функция от пола. Расчитать вероятность выживания для мужчин можно вот так:
$$ P(y_i = 1 \mid sex_i = male) = \frac{1}{1 + e^{-(-0.97 + 2.42 \cdot 1)}} = 0.19 $$А вероятность выживания для женщин так:
$$ P(y_i = 1 \mid sex_i = female) = \frac{1}{1 + e^{-(-0.97 + 2.42 \cdot 0)}} = 0.73 $$Получается удобно. Давайте добавим в модель возраст. Получим новую модель :)
# Пол нужно сделать числовым вида 1 - мужчины, 0 - женщины
X = cbind(1, as.numeric(df_train$Sex) - 1, df_train$Age)
res <- maxLik(lnL, start=c(0.5,0.5,2), y = y_sur, X = X)
summary(res)
-------------------------------------------- Maximum Likelihood estimation Newton-Raphson maximisation, 6 iterations Return code 2: successive function values within tolerance limit Log-Likelihood: -373.653 3 free parameters Estimates: Estimate Std. error t value Pr(> t) [1,] -1.068298 0.235441 -4.537 5.69e-06 *** [2,] 2.408176 0.183565 13.119 < 2e-16 *** [3,] 0.003508 0.006786 0.517 0.605 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 --------------------------------------------
Вероятность выживания $30$-летнего мужчины тогда будет расчитываться по формуле:
$$ P(y_i = 1 \mid sex_i = male, age_i = 30) = \frac{1}{1 + e^{-(-1.068 + 2.408 \cdot 1 +0.004 \cdot 30 )}} = 0.026 $$Похоже, что Джек (если его так звали) был обречён.
Всё то же самое можно сделать более удобно внутри команды glm
, которая расшифровывается как general linear model. Внутри этой команды лежит много няшных моделей. В том числе и логистическая регрессия.
# сделаем целевую переменную факторной (мы до этого сделали её числовой)
df_train <- mutate(df_train, Survived = as.factor(Survived))
# обучаем модель! Первый аргумент - табличка, второй - уравнение
# мы хотим прогнозировать выживаемость по полу, возрасту, классу и цене за билет
# третий аргумент - семейство моделей, в нашем случае это логистическая регрессия
m_logit <- glm(data=df_train, Survived~Sex+Age+Pclass+Fare,
family=binomial(link="logit"))
summary(m_logit)
Call: glm(formula = Survived ~ Sex + Age + Pclass + Fare, family = binomial(link = "logit"), data = df_train) Deviance Residuals: Min 1Q Median 3Q Max -2.5667 -0.7183 -0.4367 0.6560 2.3834 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) 3.129534 0.444319 7.043 1.88e-12 *** Sexmale -2.484069 0.203714 -12.194 < 2e-16 *** Age -0.029922 0.007894 -3.790 0.00015 *** Pclass2 -0.905203 0.310072 -2.919 0.00351 ** Pclass3 -2.118817 0.306747 -6.907 4.94e-12 *** Fare 0.001230 0.002255 0.546 0.58528 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 951.76 on 712 degrees of freedom Residual deviance: 664.78 on 707 degrees of freedom AIC: 676.78 Number of Fisher Scoring iterations: 4
Получаем протокол с оценками. Мы уже получали много протоколов. Давайте наконец научимся читать его. Первая колонка - оценка коэффициентов. Intercept это константа $\beta_0$. Второй столбец - оценка стандартных ошибок для каждого коэфициента. С их помощью можно построить доверительные интервалы для коэффициентов. Например, коэффициент перед полом с вероятностью $95\%$ будет лежать в интервале:
-2.484069 - qnorm(1 - 0.05/2)*0.203714
-2.484069 + qnorm(1 - 0.05/2)*0.203714
Этот доверительный интервал будет асимптотическим. Он строится на основе того, что любая оценка максимального правдоподобия имеет асимптотически нормальное распределение. Кстати говоря, можно посмотреть на ковариационную матрицу, которую мы в прошлый раз получали из гессиана:
vcov(m_logit)
(Intercept) | Sexmale | Age | Pclass2 | Pclass3 | Fare | |
---|---|---|---|---|---|---|
(Intercept) | 0.1974191158 | -3.775000e-02 | -2.565048e-03 | -0.0883898835 | -0.1120282475 | -4.785463e-04 |
Sexmale | -0.0377500011 | 4.149945e-02 | 7.192217e-05 | 0.0090984504 | 0.0167800795 | 4.654275e-05 |
Age | -0.0025650478 | 7.192217e-05 | 6.231689e-05 | 0.0006425882 | 0.0010204559 | 1.596519e-06 |
Pclass2 | -0.0883898835 | 9.098450e-03 | 6.425882e-04 | 0.0961444056 | 0.0644492777 | 2.983607e-04 |
Pclass3 | -0.1120282475 | 1.678008e-02 | 1.020456e-03 | 0.0644492777 | 0.0940938999 | 3.424597e-04 |
Fare | -0.0004785463 | 4.654275e-05 | 1.596519e-06 | 0.0002983607 | 0.0003424597 | 5.083539e-06 |
По диагонали стоят дисперсии. Взяв из них квадратный корень мы получим вторую колонку со стандартными ошибками из нашего протокола.
sqrt(4.149945e-02)
df_test[1,]
Survived | Pclass | Sex | Age | SibSp | Parch | Fare | fit | se.fit | residual.scale |
---|---|---|---|---|---|---|---|---|---|
0 | 3 | male | 22 | 1 | 0 | 7.25 | -2.122707 | 0.1784885 | 1 |
На пересечениях стоят ковариации между оценками. В оставшихся двух колонках проверяется гипотеза о том, что каждый из коэффициентов равен нулю. Но об этом позже. Построим прогнозы для тестовой выборки.
pr_logit <- predict(m_logit, df_test[,-1], se=TRUE) # Строим прогнозы по нашей модели на тесте
df_test <- cbind(df_test, prob_surv = sigmoid(pr_logit$fit)) # объеядиним прогнозы с тестовой табилчкой
head(df_test)
Survived | Pclass | Sex | Age | SibSp | Parch | Fare | fit | se.fit | residual.scale | sigmoid(pr_logit$fit) | pr_logit$se.fit | prob_surv | prob_surv | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 0 | 3 | male | 22 | 1 | 0 | 7.2500 | -2.1227075 | 0.1784885 | 1 | 0.10690929 | 0.1784885 | 0.10690929 | 0.10690929 |
15 | 0 | 3 | female | 14 | 0 | 0 | 7.8542 | 0.6014780 | 0.1869467 | 1 | 0.64599437 | 0.1869467 | 0.64599437 | 0.64599437 |
31 | 0 | 1 | male | 40 | 0 | 0 | 27.7208 | -0.5172934 | 0.2171901 | 1 | 0.37348534 | 0.2171901 | 0.37348534 | 0.37348534 |
32 | 1 | 1 | female | 28 | 1 | 0 | 146.5208 | 2.4720006 | 0.2909189 | 1 | 0.92215550 | 0.2909189 | 0.92215550 | 0.92215550 |
33 | 1 | 3 | female | 28 | 0 | 0 | 7.7500 | 0.1824471 | 0.1694246 | 1 | 0.54548566 | 0.1694246 | 0.54548566 | 0.54548566 |
43 | 0 | 3 | male | 28 | 0 | 0 | 7.8958 | -2.3014426 | 0.1816917 | 1 | 0.09100355 | 0.1816917 | 0.09100355 | 0.09100355 |
В колонке prob_surv
находится прогноз вероятности выжить. Давайте будем считать, что все, у кого она больше $0.5$ выживут.
y_real = as.numeric(as.character(df_test$Survived))
sum(y_real == 1*(df_test$prob_surv > 0.5))/nrow(df_test)
Точность прогноза $83\%$. Неплохо. Помните мы в самом начале построили картинку для частоты выживших и умерших и поняли, что выборка несбалансирована. Давайте посмотим какой будет тончость прогноза на тестовой выборке, если мы скажем, что все живы или все мертвы.
sum(y_real == 1)/nrow(df_test)
sum(y_real == 0)/nrow(df_test)
Если сказать, что все умрут, то точность модели составит $63\%$. Имеет смысл верить модели только если её точность намного больше этого числа. Такая метрика, количество верных прогнозов, называется Accuracy. Она интуитивно понятная, очевидная и почти неиспользуемая на практике.
Её главная проблема в том, что она бесполезна в задачах с неравными классами. Например, пусть у нас есть $50$ больных и $950$ здоровых. Мы хотим научиться их различать. Пусть наш алгоритм предсказывает, что все здоровы. В этом случае доля правлиьных ответов составит $95\%$, но алгорим окажется абсолютно бесполезным. Чтобы избежать таких эксцессов, а также учитывать, что разные типы ошибок могут иметь разную цену, строят другие две метрики. Про них можно найти более подробную информацию в почиташках.
Последний шаг для этого скрипта. Зная качество прогнозов на тестовой выборке, давайте построим прогноз для себя.
# Вбейте данные о себе:
you_data <- data.frame(Sex='male', Age=20, Pclass='3', Fare=70)
# Вероятность того, что вы окажетесь в живых
cat('Вероятность выживавания:', sigmoid(predict(m_logit, you_data, se=TRUE)$fit))
Вероятность выживавания: 0.1207168
Логисчтическая регрессия очень часто используется на практике. Многие модели, которые работают у крупных фирм в продакшене это именно она. Правда немного модернизированная.