Когнитивные технологии

Алла Тамбовцева

Парная регрессия, проверка гипотез о значимости коэффициентов

Импортируем библиотеку pandas и загрузим таблицу с данными по ссылке:

In [67]:
import pandas as pd
df = pd.read_csv("https://vincentarelbundock.github.io/Rdatasets/csv/carData/Salaries.csv")

Таблица df содержит показатели по преподавателям в США за 2008-2009 годы, включая их пол, статус, преподаваемую дисциплину, число лет после получения степени, опыт работы и заработную плату. Подробное описание см. здесь.

Посмотрим на общую информацию о датафрейме:

In [68]:
df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 397 entries, 0 to 396
Data columns (total 7 columns):
Unnamed: 0       397 non-null int64
rank             397 non-null object
discipline       397 non-null object
yrs.since.phd    397 non-null int64
yrs.service      397 non-null int64
sex              397 non-null object
salary           397 non-null int64
dtypes: int64(4), object(3)
memory usage: 21.8+ KB

Посмотрим на первые несколько строк таблицы:

In [40]:
df.head()
Out[40]:
Unnamed: 0 rank discipline yrs.since.phd yrs.service sex salary
0 1 Prof B 19 18 Male 139750
1 2 Prof B 20 16 Male 173200
2 3 AsstProf B 4 3 Male 79750
3 4 Prof B 45 39 Male 115000
4 5 Prof B 40 41 Male 141500

Переименуем столбцы, содержащие точку в названии (это не мешает работе с датафреймом, но будет мешать при включении переменных в модель, которая задается в виде текстовой формулы):

In [70]:
cols = list(df.columns)
cols[3], cols[4] = 'phd', 'service'
df.columns = cols
In [71]:
df.head()
Out[71]:
Unnamed: 0 rank discipline phd service sex salary
0 1 Prof B 19 18 Male 139750
1 2 Prof B 20 16 Male 173200
2 3 AsstProf B 4 3 Male 79750
3 4 Prof B 45 39 Male 115000
4 5 Prof B 40 41 Male 141500

Теперь приступим к созданию регрессионной модели, которая будет моделировать взаимосвязь между опытом работы в годах и заработной платой преподавателя. В статистике часто формулируют подобные отношения в терминах влияния: «как опыт работы влияет на заработную плату преподавателя». Это терминологически верно, однако стоит иметь в виду, что под влиянием имеется наличие значимой статистической направленной связи, а не причинно-следственную связь. Причинно-следственную связь позволяет выявить только специальный экспериментальный или псевдо-экспериментальный дизайн и соответствующие более сложные методы.

В чём принципиальная разница между моделированием в машинном обучении и в классической статистике/эконометрике? В задачах. Задача машинного обучения часто сводится к созданию модели или алгоритма, который позволяет предсказывать некоторый результат наилучшим образом. Другими словами, в машинном обучении важна предсказательная сила модели. В эконометрике предсказательная сила модели и её качество тоже, безусловно, важны, однако большее внимание уделяется содержательным аспектам и эффектам каждого признака на зависимую переменную. Для сравнения приведём задачи, которые могут интересовать специалиста по машинному обучению и исследователя-статистика.

Есть массив данных, включающий пол, возраст, статус, преподаваемую дисциплину, число лет после получения степени, опыт работы и заработную плату преподавателей.

Машинное обучение. Как на основе имеющихся данных построить модель, которая будет наилучшим образом предсказывать зарплату преподавателя? Нужно отобрать те переменные, добавление которых уменьшает ошибку прогноза, проверить, не страдает ли модель от переобучения, провести кросс-валидацию ...

Статистика и исследования. Как сравнить влияние опыта работы и числа лет после получения степени на заработную плату? Останется ли влияние опыта работы на заработную плату значимым, если в модели учесть пол преподавателя? Возраст преподавателя связан с зарплатой линейно или U-образно? Чем старше преподаватель, тем выше зарплата, или возраст влияет положительно только на определенном возрастном промежутке (скажем, от 30 до 40 лет)?

Начнём разбираться с подобными вопросами. Как вы уже могли заметить, выдачи Python по итогам статистических тестов очень лаконичные. Иногда их хватает, но в случае регрессионного анализа хотелось бы получать подробные результаты сразу, чтобы не вычислять дополнительные параметры отдельно. Поэтому мы не будем пользоваться функциями для линейной регрессии из стандартных модулей scipy, а импортируем функцию для обычной линейной регрессии из statsmodels:

In [54]:
from statsmodels.formula.api import ols

Для начала построим простую парную линейную регрессию, оценённую методом наименьших квадратов (ordinary least squares). Работаем с базой df, зависимая переменная salary (зарплата), независимая – service (опыт работы в годах).

In [77]:
model = ols("salary ~ service", df).fit()

Посмотрим на результаты:

In [74]:
print(model)
<statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x1a22b623c8>
In [75]:
print(model.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                 salary   R-squared:                       0.112
Model:                            OLS   Adj. R-squared:                  0.110
Method:                 Least Squares   F-statistic:                     49.85
Date:                Sun, 30 Dec 2018   Prob (F-statistic):           7.53e-12
Time:                        03:06:44   Log-Likelihood:                -4635.7
No. Observations:                 397   AIC:                             9275.
Df Residuals:                     395   BIC:                             9283.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept   9.997e+04   2416.605     41.370      0.000    9.52e+04    1.05e+05
service      779.5691    110.417      7.060      0.000     562.491     996.647
==============================================================================
Omnibus:                       25.187   Durbin-Watson:                   1.843
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               28.656
Skew:                           0.593   Prob(JB):                     5.99e-07
Kurtosis:                       3.570   Cond. No.                         36.9
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Выдача очень подробная, давайте рассмотрим всё по порядку. Начнём со второй таблицы. В ней указаны оценки коэффициентов регрессии, их стандартные ошибки, значения статистики критерия Стьюдента, используемого для проверки гипотезы о равенстве коэффициента нулю, значимость коэффициентов.

Столбец coef содержит сами значения коэффициентах при предикторах (независимых переменных, признаках). По полученным результатам можем составить уравнение модели (коэффициенты округлены):

$$ salary = 99970 + 779.6 \times service $$

Видно, что опыт работы оказывает положительный эффект на заработную плату (коэффициент положителен) и при увеличении опыта работа на $1$ год заработная плата, в среднем, увеличивается на $779.6$.

В таблице с результатами мы видим не истинные значения коэффициентов (как было бы в модели для всех преподавателей в США вообще), а их оценки, полученные на основе имеющейся выборки в 397 человек. Каждый коэффициент «снабжён» стандартной ошибкой – стандартным отклонением коэффициента, посчитанным по нашей выборке. Чтобы понять, сильно ли истинное значение коэффициента отклоняется от нуля, мы можем поделить полученный коэффициент на стандартную ошибку и оценить, насколько нетипичным является полученное значение в предположении, что предиктор не влияет на зависимую переменную. И тут мы вплотную подходим к проверке гипотезы о значимости коэффициента регрессии.

$H_0: \beta = 0$ (истинный коэффициент при предикторе незначим, равен $0$)

$H_1: \beta \ne 0$

Для проверки нулевой гипотезы используется критерий Стьюдента для проверки значимости коэффициента регрессии. Статистика критерия имеет распределение Стьюдента с $df=n-2$ степенями свободы, где $n$ – число наблюдений в выборке. Значение статистики критерия и есть t из таблицы выше. Если значение t слишком большое или слишком маленькое, то есть находится далеко в «хвостах» графика плотности распределения, то оно считается нетипичным для случая, если нулевая гипотеза верна. И тогда нулевая гипотеза отвергается, коэффициент признаётся значимым. Как понять, что значение статистики является нетипичным? Посмотреть на p-value, которое в таблице обозначено P>|t|. Если p-value маленькое, близко к нулю, это означает, что значений справа или слева (зависит от знака) от полученного t очень мало, то есть оно находится далеко в «хвостах» распределения. Если p-value большое, вероятность P>|t| большая, и это сигнализирует о том, что значение t достаточно сильно удалено от «хвостов», то есть находится в области типичных значений.

В данном случае p-value для коэффициента при переменной service равно $0$, что меньше конвенциональных $0.05$, поэтому можно сделать вывод о том, что на имеющихся данных на 5%-ном уровне значимости есть основания отвергнуть нулевую гипотезу о равенстве коэффициента нулю. Опыт работы оказывает значимое влияние на заработную плату.

Кроме самих коэффициентов и p-value в таблице приведены 95%-ные доверительные интервалы для коэффициентов регрессии. Рассмотрим доверительный интервал для коэффициента при service: [562.491, 996.647]. Проинтерпретируем. Если мы будем проводить аналогичные исследования на разных выборках объёма 397 человек и строить такую же регрессионную модель много раз, в 95% случаев коэффициент при переменной service будет принадлежать интервалу [562.491, 996.647]. Так как доверительный интервал не накрывает значение $0$, это тоже служит подтверждениме того, что коэффициент значимо отличен от нуля.

Перейдём к первой таблице. В левой части содержится сводная информация о модели, о числе наблюдений, о типе стандартных ошибок (по умолчанию обычные, не робастные, то есть не устойчивые к различным проблемам моделей). В правой части содержатся показатели качества модели.

R-squared ($R^2$) – коэффициент детерминации, показывает, какую долю дисперсии зависимой переменной объясняет построенная модель. В случае парной регрессии $R^2$ совпадает с квадратом коэффициента корреляции между предиктором и зависимой переменной. Чётких границ для интерпретации значения $R^2$ не существует, в разных науках по-разному, но $R^2$ больше 0.7 считается высоким. Здесь он явно низкий, следовательно, модель не очень хорошо объясняет изменчивость зависимой переменной.

Значения F-statistic и Prob (F-statistic) относятся к проверке гипотезы о целесообразности модели. В случае парной регрессии это не очень актуально, поскольку достаточно проверить значимость ровно одного коэффициента (значимость константы обычно не играет никакой роли), а в случае множественной регресии это полезно.

$H_0: \beta_1 = \beta_2 = \dots = \beta_k = 0$

$H_1: \exists j: \beta_j \ne 0 $

Нулевая гипотеза по сути гласит, что модель не нужна: все коэффициенты равны $0$. Естественно, при построении модели мы заинтересованы в том, чтобы такая нулевая гипотеза отвергалась. Что здесь и происходит.

Значения Log-Likelihood, AIC, BIC относятся к правдоподобию модели, это натуральный логарифм правдоподобия модели, информационный критерий Акаике и информационный критерий Байеса. Обычно данные показатели не очень информативны сами по себе, то есть для одной модели, обычно их используют для сравнения двух и более моделей. Например, есть две модели

$$ y = \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 $$$$ y = \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 $$

и нас интересует, какая модель лучше. Можем выбрать ту, у который информационные критерии AIC и BIC меньше: так как при расчёте критериев используется логарифм правдободобия, взятый со знаком минус, очевидно, что меньшее значение IC соответствует большему правдоподобию модели.

Про третью таблицу поговорим в следующей лекции, а сейчас осталось пояснить только одно: два связанных момента Covariance Type: nonrobust и Warnings:[1] Standard Errors assume that the covariance matrix of the errors is correctly specified. Эти аспекты относятся к спецификации стандартных ошибок коэффициентов. Дело в том, что по умолчанию используются обычные стандартные ошибки коэффициентов, которые не являются устойчивыми к потенциальным проблемам модели. В парной регрессии таких проблем не возникает, а вот модель множественной регрессии может страдать от мультиколлинеарности (сильной линейной связи между предикторами), от гетероскедастичности (определенных паттернов в распределении ошибок модели) и других сложностей, связанных со структурой модели. Про эти вещи мы поговорим чуть позже, а пока посмотрим, как можно выгрузить результаты регрессии.

Например, так:

In [78]:
model.summary()
Out[78]:
OLS Regression Results
Dep. Variable: salary R-squared: 0.112
Model: OLS Adj. R-squared: 0.110
Method: Least Squares F-statistic: 49.85
Date: Sun, 30 Dec 2018 Prob (F-statistic): 7.53e-12
Time: 15:54:02 Log-Likelihood: -4635.7
No. Observations: 397 AIC: 9275.
Df Residuals: 395 BIC: 9283.
Df Model: 1
Covariance Type: nonrobust
coef std err t P>|t| [0.025 0.975]
Intercept 9.997e+04 2416.605 41.370 0.000 9.52e+04 1.05e+05
service 779.5691 110.417 7.060 0.000 562.491 996.647
Omnibus: 25.187 Durbin-Watson: 1.843
Prob(Omnibus): 0.000 Jarque-Bera (JB): 28.656
Skew: 0.593 Prob(JB): 5.99e-07
Kurtosis: 3.570 Cond. No. 36.9


Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Или так:

In [81]:
model.summary2()
Out[81]:
Model: OLS Adj. R-squared: 0.110
Dependent Variable: salary AIC: 9275.3764
Date: 2018-12-30 15:55 BIC: 9283.3443
No. Observations: 397 Log-Likelihood: -4635.7
Df Model: 1 F-statistic: 49.85
Df Residuals: 395 Prob (F-statistic): 7.53e-12
R-squared: 0.112 Scale: 8.1669e+08
Coef. Std.Err. t P>|t| [0.025 0.975]
Intercept 99974.6529 2416.6053 41.3699 0.0000 95223.6362 104725.6695
service 779.5691 110.4169 7.0602 0.0000 562.4908 996.6475
Omnibus: 25.187 Durbin-Watson: 1.843
Prob(Omnibus): 0.000 Jarque-Bera (JB): 28.656
Skew: 0.593 Prob(JB): 0.000
Kurtosis: 3.570 Condition No.: 37

При желании можно запросить отдельные показатели модели:

In [82]:
model.conf_int()  # 95% доверительные интервалы для коэффициентов
Out[82]:
0 1
Intercept 95223.636243 104725.669490
service 562.490765 996.647462
In [83]:
model.conf_int(0.90)  # 90% доверительные интервалы для коэффициентов
Out[83]:
0 1
Intercept 99670.783695 100278.522038
service 765.685050 793.453178
In [85]:
model.pvalues  # p-value для коэффициентов
Out[85]:
Intercept    1.186078e-145
service       7.528739e-12
dtype: float64

И многое другое. Все доступные атрибуты и методы можно посмотреть, набрав точку и нажав Tab. Многие из них мы ещё разберём далее, но на более интересных моделях.