Продолжим исследовать данные из параграфа 10.2 о популярности песен на YouTube.
В прошлом параграфе мы говорили о том, что мы выдвигаем гипотезы до того как приступаем к моделированию. Давайте попробуем исследовать гипотезу: «Правда ли, что более танцевальные песни обычно короче?». Мы выдвигаем такую гипотезу, потому что предполагаем, что танцевать шесть минут под одну песню может быть довольно скучно.
Для проверки мы должны установить с помощью регрессии наличие или отсутствие связи между длительностью песни и тем, что она характеризуется как танцевальная. При этом нам важно рассмотреть и другие факторы, которые могут быть связаны с зависимой переменной (т.е длиной песни).
Для этого предложим ещё несколько гипотез:
короткиепеснигромче, потому что артисты стремятся сделать их максимально привлекающими внимание на короткий промежуток времени;- если в песне есть
слова, то она будеткороче, потому что в ней будет меньше инструментальных проигрышей; энергичныепесникороче, потому что в таких песнях будет меньше медленных частей.
Если суммировать, то мы предполагаем, что с одной стороны у нас будут короткие, громкие и энергичные песни, а с другой — длинные, тихие и меланхоличные. Но прежде чем проверить наши гипотезы — давайте подготовим данные.
Подготовка данных
Для начала составим впечатление о данных. Раньше мы выводили несколько строк со всеми переменными, но в больших датафреймах это не так удобно — потому здесь мы выведем названия колонок:
import pandas as pd
df = pd.read_csv('Spotify_Youtube.csv')
df.columns
Index(['Artist', 'Url_spotify', 'Track', 'Album', 'Album_type',
'Uri', 'Danceability', 'Energy', 'Key', 'Loudness', 'Speechiness',
'Acousticness', 'Instrumentalness', 'Liveness', 'Valence', 'Tempo',
'Duration_ms', 'Url_youtube', 'Title', 'Channel', 'Views', 'Likes',
'Comments', 'Description', 'Licensed', 'official_video', 'Stream'],
dtype='object')
Всего в наборе данных 26 переменных. Полное описание набора данных доступно тут. Мы не будем использовать все — возьмём только те, что нужны для анализа:
- Duration.
- Danceability.
- Loudness.
- Speechiness.
- Energy.
Сделаем два преобразования: одно для удобства, другое — содержательное.
Длительность у нас указана в миллисекундах. Поэтому давайте для удобства превратим миллисекунды в секунды.
df['Duration'] = df['Duration_ms'] / 1000
Содержательное изменение следующее: так как в данных есть и очень короткие, и очень длинные песни, то нам стоит исключить все треки, которые длятся меньше 30 секунд и больше 300 секунд. Нам это нужно, чтобы избежать выбросов — то есть влияния нетипичных песен на результаты регрессии.
df = df[(df['Duration'] < 300) & (df['Duration'] > 30)]
Далее выпишем переменные (колонки), в которых собраны относящиеся к нашим гипотезам данные (порядок совпадает с гипотезами). Затем с помощью метода dropna() сразу удалим пропуски, которые мы не сможем учесть в модели. И, наконец, выведем описательные статистики для наших переменных. Последнее нужно и для удобства восприятия данных и для проверки допущений для проведения линейной регрессии.
df = df[['Duration', 'Danceability', 'Loudness', 'Speechiness', 'Energy']].dropna()
print(df.describe())
Duration Danceability Loudness Speechiness Energy
count 18646.000000 18646.000000 18646.000000 18646.000000 18646.000000
mean 206.110522 0.626312 -7.536997 0.097840 0.638738
std 45.703952 0.161589 4.577807 0.113562 0.213155
min 30.985000 0.000000 -46.251000 0.000000 0.000055
25% 176.144000 0.527000 -8.649750 0.035900 0.511000
50% 207.264500 0.643000 -6.422000 0.051300 0.670000
75% 237.559000 0.744000 -4.860000 0.107000 0.801000
max 299.960000 0.975000 0.920000 0.964000 1.000000
Все наши независимые переменные непрерывные и расположены в диапазоне от нуля до единицы — кроме громкости, она измеряется по своей шкале в децибелах. Все независимые переменные варьируются в каком-то диапазоне, поэтому мы можем предположить наличие какой-то связи между ними.
Теперь мы можем перейти непосредственно к моделированию.
Моделирование
Для расчётов нам нужно загрузить пакет statsmodels и объект OLS, благодаря которому мы и рассчитаем линейную регрессию. Второе необязательно, но позволяет не писать sm.OLS.
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
Для моделирования мы сохраняем данные в две переменные Y и X. Первая содержит зависимую переменную, вторая — независимые. Мы делаем это для того, чтобы код, который описывает модель, можно было записать без указания на какие-то конкретные колонки в датафрейме. Так нам будет проще переиспользовать его, когда мы будем запускать моделирование несколько раз подряд.
С помощью метода sm.add_constant мы добавляем b_0 (константу). Это техническое действие, но важное, потому что без него результаты моделирования будут некорректными.
X = df[['Danceability', 'Loudness', 'Speechiness', 'Energy']]
Y = df['Duration']
X = sm.add_constant(X, prepend=False)
Далее мы можем создать объект модели, а потом с помощью метода fit вычислить наилучшую регрессионную линию. Наконец, мы сохраняем результаты моделирования и выводим их с помощью метода summary.
Важно:
fitотносится к объекту модели, аsummaryк результатам.
model = OLS(Y, X)
res = model.fit()
print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Duration R-squared: 0.046
Model: OLS Adj. R-squared: 0.046
Method: Least Squares F-statistic: 225.9
Date: Wed, 03 May 2023 Prob (F-statistic): 1.01e-189
Time: 11:44:03 Log-Likelihood: -97284.
No. Observations: 18646 AIC: 1.946e+05
Df Residuals: 18641 BIC: 1.946e+05
Df Model: 4
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Danceability -27.5230 2.208 -12.463 0.000 -31.851 -23.194
Loudness 2.0397 0.111 18.413 0.000 1.823 2.257
Speechiness -39.7622 2.954 -13.461 0.000 -45.552 -33.972
Energy -3.9776 2.288 -1.739 0.082 -8.462 0.507
const 245.1527 2.712 90.402 0.000 239.837 250.468
==============================================================================
Omnibus: 69.121 Durbin-Watson: 1.304
Prob(Omnibus): 0.000 Jarque-Bera (JB): 69.670
Skew: -0.147 Prob(JB): 7.44e-16
Kurtosis: 2.938 Cond. No. 97.9
==============================================================================
Выше мы видим множество показателей. Большинство из них нужны только для глубокой интерпретации результатов, которую мы делать не будем, потому что это все же вводный учебник, а не продвинутый.
Давайте обратим внимание на те показатели, о которых мы говорили в прошлом параграфе. Нас интересуют:
R-squared— это коэффициент детерминации;- колонка
coefс коэффициентамиb; - колонка
P>|t|сp_value— его мы используем для оценки значимости коэффициента.
Давайте начнём с простого шага в интерпретации: мы видим, что p_value у переменной Energy незначим (он больше 0.05). Значит, мы исключим её из модели.
X = df[['Danceability', 'Loudness', 'Speechiness']]
Y = df['Duration']
X = sm.add_constant(X, prepend=False)
model = OLS(Y, X)
res = model.fit()
print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Duration R-squared: 0.046
Model: OLS Adj. R-squared: 0.046
Method: Least Squares F-statistic: 300.2
Date: Wed, 03 May 2023 Prob (F-statistic): 2.43e-190
Time: 11:49:32 Log-Likelihood: -97286.
No. Observations: 18646 AIC: 1.946e+05
Df Residuals: 18642 BIC: 1.946e+05
Df Model: 3
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Danceability -27.2695 2.204 -12.375 0.000 -31.589 -22.950
Loudness 1.8999 0.076 24.945 0.000 1.751 2.049
Speechiness -40.1560 2.945 -13.634 0.000 -45.929 -34.383
const 241.4378 1.670 144.591 0.000 238.165 244.711
==============================================================================
Omnibus: 71.496 Durbin-Watson: 1.305
Prob(Omnibus): 0.000 Jarque-Bera (JB): 72.155
Skew: -0.150 Prob(JB): 2.15e-16
Kurtosis: 2.945 Cond. No. 83.2
==============================================================================
Как видите, после удаления Energy коэффициент детерминации совсем не изменился. Это показывает, что переменная Energy совсем никак не объясняла изменчивость нашей зависимой переменной. Простыми словами — «энергичность» песни никак не влияет на её длительность.
Затем нам нужно построить модель с единственной независимой переменной, которая была в нашей самой первоначальной гипотезе. То есть мы проверим связь между длительностью и танцевальностью песни. Нам это нужно, чтобы понять, насколько добавление дополнительных переменных оказалось полезным:
X = df[['Danceability']]
Y = df['Duration']
X = sm.add_constant(X, prepend=False)
model = OLS(Y, X)
res = model.fit()
print(res.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Duration R-squared: 0.003
Model: OLS Adj. R-squared: 0.003
Method: Least Squares F-statistic: 49.39
Date: Wed, 31 May 2023 Prob (F-statistic): 2.17e-12
Time: 12:01:05 Log-Likelihood: -97701.
No. Observations: 18646 AIC: 1.954e+05
Df Residuals: 18644 BIC: 1.954e+05
Df Model: 1
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Danceability -14.5388 2.069 -7.028 0.000 -18.594 -10.484
const 215.2163 1.338 160.840 0.000 212.594 217.839
==============================================================================
Omnibus: 369.810 Durbin-Watson: 1.245
Prob(Omnibus): 0.000 Jarque-Bera (JB): 394.235
Skew: -0.341 Prob(JB): 2.47e-86
Kurtosis: 3.206 Cond. No. 8.66
==============================================================================
Тут мы снова смотрим на коэффициент детерминации. Он стал значительно меньше (0.003 вместо 0.046) — значит, от дополнительных переменных (Loudness и Speechiness) всё же была польза.
Поэтому давайте рассмотрим коэффициенты наилучшей (второй) модели. Продублируем их ниже для наглядности.
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Danceability -27.2695 2.204 -12.375 0.000 -31.589 -22.950
Loudness 1.8999 0.076 24.945 0.000 1.751 2.049
Speechiness -40.1560 2.945 -13.634 0.000 -45.929 -34.383
const 241.4378 1.670 144.591 0.000 238.165 244.711
==============================================================================
Напомним, что значение коэффициента соответствует изменению зависимой переменной при изменении независимой на единицу при прочих (других независимых переменных) равных. Знак показывает направление изменения. Получается, что если значение танцевальности песни увеличится на единицу, то длительность уменьшится на 27 секунд. Это значит, что мы можем подтвердить нашу основную гипотезу: песни, которые подходят для танцев короче.
Гипотезу о связи энергичности и длины мы отвергли, потому что коэффициент был незначим. У нас были ещё две дополнительные гипотезы:
короткиепеснигромче;- если в песне есть
слова, то онакороче.
Мы были правы относительно обеих: знак положительный, значит уменьшение громкости на один децибел связано с уменьшением длины песни на две секунды. Отсутствие слов связано с увеличением длины песни на сорок секунд.
Проверка допущений
Не менее важный этап в моделировании — диагностика. Здесь мы смотрим, насколько наша модель надежная в предсказаниях, и соблюдает ли она все допущения метода.
Проверим допущение о мультиколлинеарности. Переменные не должны быть коррелированы между друг другом. Мы можем использовать VIF (Variance Inflation Factor). Если он будет больше четырех (ещё одно консенсусное значение, такое же как p-value или минимальное число наблюдений в группе), то это значит, что эта переменная коррелирована с какой-то из присутствующих и её нужно убрать из модели.
Что нужно сделать:
- импортируем метод
variance_inflation_factor; - сохраним в переменную
Xлист с независимыми переменными из модели; - создадим пустой объект
vif_data, а потом в колонку"feature"включим переменные изX; - для вычисления нужного коэффициента создадим колонку
"VIF"в объектеvif_data; - рассчитаем для каждого значения в
Xпоказатель, включив параметр длиныi, который рассчитывается с помощью последующей простого цикла.
from statsmodels.stats.outliers_influence import variance_inflation_factor
X = df[['Danceability', 'Loudness', 'Speechiness']]
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 Danceability 1.186238
1 Loudness 1.137306
2 Speechiness 1.046675
3 const 26.087156
Мультиколлинеарности не наблюдается, мы можем оставить все переменные.
Поздравляем, вы научились вычислять линейную регрессию! Вы сделали очень большой шаг — от дескриптивного анализа и простых методов к сложным моделям, которые учитывают несколько различных факторов.
На этом хорошие новости заканчиваются — следующая глава будет ещё чуть сложнее. В ней мы попробуем рассмотреть ещё один вариант регрессионого анализа. Но сперва поговорим о том, как оценить качество данных.