В этом параграфе мы продолжим использовать данные о персонажах из вселенной «Игры Престолов». Давайте рассчитаем, какие факторы увеличивают шанс выжить для персонажа. Загрузим набор данных и попробуем кратко охарактеризовать его.
import pandas as pd
data = pd.read_csv('character_predictions_pose.csv')
Всего в выборке 1946 наблюдений. Этого вполне достаточно для построения логистической регрессии, если не учитывать отсутствующие значения в тех переменных, которые мы будем использовать.
Отберем подходящие для предсказания переменные. Нам нужны гипотезы, из которых мы будем исходить при отборе:
- шансы женщин на выживание выше;
- у замужних и женатых персонажей шансы выжить больше;
- чем больше персонаж убил людей, тем меньше его шанс выжить;
- титул дворянина увеличивает шансы выжить;
- чем популярнее персонаж, тем больше его шанс выжить.
Гипотезы определяют выбранные переменные. В нашем случае будет одна зависимая переменная isalive и пять независимых переменных ('male', 'ismarried', 'numdeadrelations', 'isnoble', 'popularity'), каждая из которых соответствует одной гипотезе. Они помогут нам предсказать значения шансов на выживание.
got = data[['isalive', 'male', 'ismarried', 'numdeadrelations', 'isnoble', 'popularity']]
print(got)
Для логистической регрессии зависимая переменная должна быть бинарной. Поскольку мы хотим предсказать шанс смерти, то наша зависимая переменная — isalive, которая сообщает о том, выжил персонаж (True) или нет (False). Преобразуем все наши переменные в категориальные.
got['isalive'] = got['isalive'].astype('int')
got['male'] = got['male'].astype('int')
got['isnoble'] = got['isnoble'].astype('int')
got['ismarried'] = got['ismarried'].astype('int')
print(got)
Сбалансированность данных
Если мы говорим о содержании переменных, то нас интересуют две проблемы:
- наличие значений в выбранных для моделирования переменных;
- отсутствие ситуаций, когда у нас одна из групп представлена всего несколькими значениями.
Про первую мы говорили ранее. В таком случае мы будем должны полностью исключить наблюдения с пропущенными значениями в выбранных нами для анализа переменных или заполнить пропуски, используя методы импутации данных. Про вторую проблему стоит поговорить подробнее.
🔍 Несбалансированность данных — свойство распределения категориальных данных, где один класс представлен существенно больше, чем все остальные в выборке.
В таких случаях нам бы пришлось бы исключить некоторые наблюдения из всей выборки, чтобы сделать сравнение более корректным. Это достаточно сложная техника работы с данными, мы не будем останавливаться на ней в этом учебнике. Если мы уверены в том, что собранные данные отражают генеральную совокупность, то может не делать ничего дополнительно. Наши данные содержат всех героев без исключения, поэтому мы можем ничего не делать.
Визуализируем категориальные переменные с помощью столбчатой диаграммы. С помощью метода set добавим название графика и подпишем оси.
import matplotlib.pyplot as plt
got['isalive'].value_counts().plot(kind='bar').set(title='Пропорция выживших персонажей', xlabel='Персонаж выжил?', ylabel='Количество')
Пропорция между выжившими и умершими: примерно один к трём. Такое различие в количестве наблюдений между категориями называется несбалансированностью данных. Потенциально это может означать, что в данных может быть какой-то неучтенный перекос между группами.
Теперь давайте посмотрим на распределения зависимых переменных.
got['male'].value_counts().plot(kind='bar').set(title='Пропорция мужчин и женщин', xlabel='Персонаж - мужчина?', ylabel='Количество')
Соотношение мужчин и женщин в выборке 3 к 2. Такое распределение относительно пропорционально, поэтому мы сможем достоверно оценить связь между гендером и шансом выжить.
Взглянем на свободных и состоящих в браке персонажей.
got['ismarried'].value_counts().plot(kind='bar').set(title='Пропорция замужних персонажей', xlabel='Супруг(а) у персонажа', ylabel='Количество')
Персонажи в браке встречаются реже в серии. Это может быть не таким надежным предиктором в предсказании шансов, но мы можем попробовать его использовать.
Визуализируем соотношение титулованых и нетитулованых персонажей вселенной.
got['isnoble'].value_counts().plot(kind='bar').set(title='Пропорция титулованных персонажей', xlabel='Дворянский титул персонажа', ylabel='Количество')
Относительно равное распределение по классам переменной.
Если у нас есть пара непрерывных переменных, то можно отразить их зависимость на графике отношений. Если их значения расположены пропорционально друг другу, то это может означать, что они коррелированы. Так можно предварительно проверить подозрения на мультиколлинеарность переменных, так как нам нужно, чтобы переменные были независимы друг к другу.
plt.scatter(got['popularity'], got['numdeadrelations'])
plt.title('Отношение популярности и количества совершенных убийств')
plt.xlabel('Рейтинг популярности')
plt.ylabel('Количество убийств совершенных персонажем')
plt.show()
На графике видно, что их распределение друг к другу практически случайно. Трудно выявить линейную зависимость. Дополнительно можно провести корреляционный тест.
График может наглядно показать сильную корреляцию или её полное отсутствие для переменной, но для неочевидных случаев лучше тест.
import numpy as np
np.corrcoef(got['popularity'], got['numdeadrelations'])
array([[1, 0.61498378],
[0.61498378, 1])
Результаты теста показывают, что корреляция средняя. Навскидку можно предположить, что дело в том, какое количество экранного времени есть у персонажа: тогда он может оказаться популярнее и совершить больше убийств. Но тут как таковая связь не прослеживается.
🔍 Мультиколлинеарность — наличие линейной зависимости между объясняющими переменными
Теоретически, мы должны проверить каждую пару переменных между собой. Это может быть сложно сделать для каждой пары, поэтому мы покажем ниже простой способ.
Моделирование
Данные подготовлены. Подгрузим необходимые пакеты и команды, которые помогут построить модель и провести над ней диагностику.
import statsmodels.api as sm
from statsmodels.genmod.generalized_linear_model import GLM
from statsmodels.genmod import families
import statsmodels.stats.tests.test_influence
Мы строим столько моделей, сколько предикторов у нас есть. В нашем случае это пять моделей. Начинаем мы с модели, в которой один предиктор, дальше мы последовательно добавляем оставшиеся. Мы ориентируемся на показатели правдоподобия (Log-Likelihood), чтобы понять, какая модель лучше подходит для объяснения данных.
Заметим, что если при включении дополнительного предиктора, другая независимая переменная превратится из незначимой в значимую, то это будет поводом задуматься о нарушении допущения о некоррелированности предикторов.
Если две переменных сильно связаны между собой, то их коэффициенты будут неопределенными, то есть произвольное взаимное изменение коэффициентов перед ними будет приводить к той же модели.
Давайте попробуем объяснить на пальцах что будет происходить: коэффициенты у этих переменных будут меняться, а сама модель в общем — нет.
Для создания модели используем команду GLM из sm. Запишем результаты команды в объект model_1. Первый параметр – зависимая переменная выживаемости got['isalive']. Во втором параметре указывается лист со всеми независимыми переменными. Начнем с переменной гендера. Последним параметром, укажем family=families.Binomial(), чтобы указать, что мы используем биномиальную логистическую регрессию.
Наконец, мы инициализируем с помощью метода fit() создание модели. Наконец, выведем результаты.
model_1 = sm.GLM(
got["isalive"],
got["male"],
family=families.Binomial(),
).fit()
print(model_1.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1945
Model Family: Binomial Df Model: 0
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1254.3
Date: Wed, 12 Apr 2023 Deviance: 2508.6
Time: 16:27:16 Pearson chi2: 1.95e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1676
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
male 0.8257 0.063 13.190 0.000 0.703 0.948
==============================================================================
Взглянем на нижнюю строчку. Мы помним, что знак коэффициента показывает направление связи: шансы выжить у мужчин значительно выше. Коэффициент логарифмирован, поэтому мы не можем сказать о том, насколько именно изменяются шансы.
Показатель, по которому определяется порог значимости - p-value. Если он больше 0.05, то мы отвергаем гипотезу о том, что переменная незначимо влияет на зависимую переменную. В нашем случае это не так — результат значим.
Пока мы видим только то, что вероятность выжить у мужчин выше, и это значимый результат. А вот насколько именно выше — нет, но мы поговорим про это ниже.
Теперь построим другую модель.
model_2 = sm.GLM(
got["isalive"],
got[["male", "ismarried"]],
family=families.Binomial(),
).fit()
print(model_2.summary())
model_3 = sm.GLM(
got["isalive"],
got[["male", "ismarried", "numdeadrelations"]],
family=families.Binomial(),
).fit()
print(model_3.summary())
model_4 = sm.GLM(
got["isalive"],
got[["male", "ismarried", "numdeadrelations", "isnoble"]],
family=families.Binomial(),
).fit()
print(model_4.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1944
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1249.0
Date: Thu, 13 Apr 2023 Deviance: 2498.0
Time: 14:31:43 Pearson chi2: 1.96e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1613
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
male 0.7780 0.064 12.126 0.000 0.652 0.904
ismarried 0.4327 0.136 3.192 0.001 0.167 0.698
==============================================================================
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1943
Model Family: Binomial Df Model: 2
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1228.1
Date: Thu, 13 Apr 2023 Deviance: 2456.2
Time: 14:31:43 Pearson chi2: 1.97e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1366
Covariance Type: nonrobust
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
male 0.8463 0.066 12.836 0.000 0.717 0.976
ismarried 0.6514 0.145 4.500 0.000 0.368 0.935
numdeadrelations -0.2386 0.041 -5.882 0.000 -0.318 -0.159
====================================================================================
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1942
Model Family: Binomial Df Model: 3
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1205.4
Date: Thu, 13 Apr 2023 Deviance: 2410.7
Time: 14:31:43 Pearson chi2: 1.99e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1103
Covariance Type: nonrobust
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
male 0.5675 0.077 7.352 0.000 0.416 0.719
ismarried 0.4034 0.150 2.691 0.007 0.110 0.697
numdeadrelations -0.2557 0.041 -6.216 0.000 -0.336 -0.175
isnoble 0.6343 0.096 6.641 0.000 0.447 0.821
====================================================================================
Все вновь добавленные переменные оказываются значимы. Их интерпретация аналогична той, которую мы проделали выше с первой переменной. Чтобы выжить лучше быть в браке, дворянином и менее агрессивным.
Давайте добавим последний предиктор.
model_5 = sm.GLM(
got["isalive"],
got[["male", "ismarried", "numdeadrelations", "isnoble", "popularity"]],
family=families.Binomial(),
).fit()
print(model_5.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1942
Model Family: Binomial Df Model: 3
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1205.4
Date: Thu, 13 Apr 2023 Deviance: 2410.7
Time: 14:31:46 Pearson chi2: 1.99e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1103
Covariance Type: nonrobust
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
male 0.5675 0.077 7.352 0.000 0.416 0.719
ismarried 0.4034 0.150 2.691 0.007 0.110 0.697
numdeadrelations -0.2557 0.041 -6.216 0.000 -0.336 -0.175
isnoble 0.6343 0.096 6.641 0.000 0.447 0.821
====================================================================================
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: isalive No. Observations: 1946
Model: GLM Df Residuals: 1941
Model Family: Binomial Df Model: 4
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -1205.3
Date: Thu, 13 Apr 2023 Deviance: 2410.6
Time: 14:31:46 Pearson chi2: 1.99e+03
No. Iterations: 4 Pseudo R-squ. (CS): -0.1103
Covariance Type: nonrobust
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
male 0.5746 0.081 7.109 0.000 0.416 0.733
ismarried 0.4129 0.153 2.693 0.007 0.112 0.713
numdeadrelations -0.2471 0.050 -4.920 0.000 -0.346 -0.149
isnoble 0.6379 0.096 6.622 0.000 0.449 0.827
popularity -0.1206 0.409 -0.295 0.768 -0.923 0.682
====================================================================================
Показатель правдоподобия и коэффициенты изменились не сильно. Более того, последняя переменная не проходит порог значимости в p-value < 0.05. Оно равно 0.768 и это говорит о том, что результат незначим. Это значит, что мы можем исключить её и использовать как финальный результат предпоследнюю модель (model_4).
Выберем лучшую модель, опираясь на показатель правдоподобия Log-Likelihood. Чем ближе эта метрика к нулю тем лучше модель описывает данные. Мы видим что у модели с четыремя предикторами — наилучший показатель правдоподобия.
Далее, посмотрим на отношения шансов без логарифма, чтобы узнать, насколько именно меняется шанс в зависимости от смены значений в предикторах. Для этого импортируем пакет numpy, извлечем из модели её коэффициенты и экспонируем.
model_odds = pd.DataFrame(np.exp(model_5.params), columns= ['Odds Ratio'])
print(model_odds)
|
--- |
Odds Ratio |
|
male |
1.776362 |
|
ismarried |
1.511153 |
|
numdeadrelations |
0.781058 |
|
isnoble |
1.892507 |
|
popularity |
0.886407 |
Давайте их интерпретируем. Если Odds Ratio = 1, то отношения шансов равные для всех классов. Если Odds Ratio > 1, то отношения шансов при повышении значения переменной вырастает. Если Odds Ratio < 1*`, то отношения шансов при повышении значения переменной понижается.
К примеру, если персонаж мужчина — его шанс выжить повышается на 77,6%. Все, что идет после 1 в Odds Ratio можно интерпретировать как проценты.
С каждым убийством других персонажей, шанс героя выжить понижается почти на 12%. Быть более миролюбивым в «Игре Престолов» гораздо выгоднее.
Диагностика и надежность предсказания модели
Проверка допущения о мультиколлинеарности полностью идентична той, которую мы делали в параграфе 10.3
from statsmodels.stats.outliers_influence import variance_inflation_factor
X = got[['male', 'ismarried', 'numdeadrelations', 'isnoble', 'popularity']]
vif_data = pd.DataFrame()
vif_data["feature"] = X.columns
vif_data["VIF"] = [variance_inflation_factor(X.values, i)
for i in range(len(X.columns))]
print(vif_data)
feature VIF
0 male 1.732705
1 ismarried 1.270200
2 numdeadrelations 1.689705
3 isnoble 1.743861
4 popularity 2.244383
Мы видим, что мультиколлинеарности тут нет, так как ни для одной переменной значение не оказывается больше консенсусного значения (четырех).
Дополнительно проверим модель на отклоняющиеся значения. Нарисуем график дистанции Кука. Если дистанция Кука превышает 0.5 - значит какое-то из значений искажает наши предсказания, выбиваясь из общего ряда. Иногда нам важно оставлять такие значения (например, если у нас мало наблюдений), но чаще всего их просто убирают.
🔍 Дистанция Кука - это мера, определяющая отклоняющую силу каждого предсказанного значения.
С помощью get_influence запишем в объект infl коэффициенты влияния предсказанных значений из модели model_5. После этого, в fig запишем эти коэффициенты и с помощью plot_index отразим их на графике. Внутри метода укажем параметры y_var='cooks'.
infl = model5.get_influence()
fig = infl.plot_index(y_var="cooks")
fig.tight_layout()
По шкале Y отражена дистанция Кука для каждого предсказанного значения в модели и индекс в датафрейме. Значение считается влиятельным, если его значение больше 0.5. Самое влиятельное наблюдение отклоняется всего лишь на 0.0, поэтому в нашем случае делать ничего не надо.
Вот и всё! Мы составили модель, которая предсказывает вероятность выжить для персонажа из «Игры престолов». Если персонаж женатый мужчина-дворянин, то его шансы выжить существенно выше, особенно если он убил не слишком много людей.
В следующем параграфе мы немного выдохнем и узнаем, как работать с текстовыми, а не табличными данными — в том числе в Python.