#!/usr/bin/env python
# coding: utf-8

# # Когнитивные технологии
# 
# *Алла Тамбовцева*
# 
# ## Парная регрессия, проверка гипотез о значимости коэффициентов

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

# In[67]:


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


# Таблица `df` содержит показатели по преподавателям в США за 2008-2009 годы, включая их пол, статус, преподаваемую дисциплину, число лет после получения степени, опыт работы и заработную плату. Подробное описание см. [здесь](https://vincentarelbundock.github.io/Rdatasets/doc/carData/Salaries.html).

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

# In[68]:


df.info()


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

# In[40]:


df.head()


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

# In[70]:


cols = list(df.columns)
cols[3], cols[4] = 'phd', 'service'
df.columns = cols


# In[71]:


df.head()


# Теперь приступим к созданию регрессионной модели, которая будет моделировать взаимосвязь между опытом работы в годах и заработной платой преподавателя. В статистике часто формулируют подобные отношения в терминах влияния: «как опыт работы влияет на заработную плату преподавателя». Это терминологически верно, однако стоит иметь в виду, что под влиянием имеется наличие значимой статистической направленной связи, а не причинно-следственную связь. Причинно-следственную связь позволяет выявить только специальный экспериментальный или псевдо-экспериментальный дизайн и соответствующие более сложные методы.
# 
# В чём принципиальная разница между моделированием в машинном обучении и в классической статистике/эконометрике? В задачах. Задача машинного обучения часто сводится к созданию модели или алгоритма, который позволяет предсказывать некоторый результат наилучшим образом. Другими словами, в машинном обучении важна предсказательная сила модели. В эконометрике предсказательная сила модели и её качество тоже, безусловно, важны, однако большее внимание уделяется содержательным аспектам и эффектам каждого признака на зависимую переменную. Для сравнения приведём задачи, которые могут интересовать специалиста по машинному обучению и исследователя-статистика.
# 
# Есть массив данных, включающий пол, возраст, статус, преподаваемую дисциплину, число лет после получения степени, опыт работы и заработную плату преподавателей.
# 
# *Машинное обучение.*  Как на основе имеющихся данных построить модель, которая будет наилучшим образом предсказывать зарплату преподавателя? Нужно отобрать те переменные, добавление которых уменьшает ошибку прогноза, проверить, не страдает ли модель от переобучения, провести кросс-валидацию ...
# 
# *Статистика и исследования.* Как сравнить влияние опыта работы и числа лет после получения степени на заработную плату? Останется ли влияние опыта работы на заработную плату значимым, если в модели учесть пол преподавателя? Возраст преподавателя связан с зарплатой линейно или 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)


# In[75]:


print(model.summary())


# Выдача очень подробная, давайте рассмотрим всё по порядку. Начнём со второй таблицы. В ней указаны оценки коэффициентов регрессии, их стандартные ошибки, значения статистики критерия Стьюдента, используемого для проверки гипотезы о равенстве коэффициента нулю, значимость коэффициентов. 
# 
# Столбец `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()


# Или так:

# In[81]:


model.summary2()


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

# In[82]:


model.conf_int()  # 95% доверительные интервалы для коэффициентов


# In[83]:


model.conf_int(0.90)  # 90% доверительные интервалы для коэффициентов


# In[85]:


model.pvalues  # p-value для коэффициентов


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