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

# # Когнитивные технологии
# 
# *Алла Тамбовцева, НИУ ВШЭ*
# 
# ## Фиктивные (дамми) переменные и эффекты взаимодействия в линейных моделях

# Продолжим работать с базой по заработной плате: загрузим необходимые библиотеки и сам файл, переименуем столбцы с точкой в названии:

# In[1]:


from statsmodels.formula.api import ols 
import pandas as pd


# In[2]:


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

# опять переименуем столбцы с точкой
cols = list(df.columns)
cols[3], cols[4] = 'phd', 'service'
df.columns = cols

# посмотрим
df.head()


# Попробуем учесть в модели пол преподавателя. В нашей таблице переменная `sex` является текстовой (два значения `female` и `male`), поэтому не очень ясно, как включать её в модель. Обычно в таких случаях переменную превращают в бинарную. Но в данном случае этого делать не нужно, Python выполнит преобразования автоматически: упорядочит два значения по алфавиту, первому присвоит значение $0$, второму – значение $1$.
# 
# В предыдущих моделях заработная плата у нас измерялась как есть, в условных единицах, но для качества модели и для интерпретации будет лучше, если разброс значений зависимой переменной будет меньше. Давайте добавим переменную `salary_th` – заработную плату в тысячах, а затем будем её использовать в качестве зависимой переменной.

# In[3]:


df['salary_th'] = df['salary'] / 1000


# Построим модель, где предикторами являются переменные `service`, `phd`, и `sex`:

# In[4]:


model1 = ols('salary_th ~ service + phd + sex', df).fit()
print(model1.summary())


# В модели добавился предиктор `sex[T.Male]`, коэффициент при котором показывает разницу в заработной плате мужчин и женщин. Так как в качестве базовой категории (значения $0$, категории, с которой мы сравниваем остальные) используется значение `female`, коэффициент при `sex` показывает, на сколько больше/меньше заработная плата у мужчин по сравнению с женщинами. 
# 
# Запишем уравнение с учётом коэффициентов модели:

# $$
# salary = 82.9 + 8.46 \cdot sex - 0.65 \cdot service + 1.56 \cdot phd
# $$

# Известно, что для преподавателей женского пола $sex=0$, для преподавателей мужского пола $sex=1$. Учтём этот факт и запишем уравнения отдельно для каждого пола:

# $$
# salary_{female} = 82.9 - 0.65 \cdot service + 1.56 \cdot phd
# $$
# 
# $$
# salary_{male} = 82.9 + 8.46 - 0.65 \cdot service + 1.56 \cdot phd
# $$
# 
# $$
# salary_{male} = 91.36 - 0.65 \cdot service + 1.56 \cdot phd 
# $$

# Видно, что в среднем, при прочих равных, заработная плата у мужчин на 8.46 тысяч выше, чем заработная плата у женщин. Это и есть коэффициент при `sex[T.Male]` в таблице. Теперь пойдём дальше и попробуем учесть не только влияние пола преподавателя, но и тот факт, что влияние опыта работы у мужчин и женщин неодинаково. Добавим в модель *эффект взаимодействия* – предиктор, представляющий собой произведение двух переменных, в данном случае произведение `sex` и `service`. В общем виде уравнение модели выглядит так:
# 
# $$
# salary = \beta_0 + \beta_1 \cdot sex + \beta_2 \cdot service + \beta_3 \cdot sex \cdot service + \beta_4 \cdot phd.
# $$
# 
# *Важно:* чтобы избежать смещения оценок коэффициентов и проблемы пропущенных переменных, переменные которые  входят в эффект взаимодействия, должны быть включены в модель по отдельности. Так, здесь мы не можем убрать слагаемые $\beta_1 \cdot sex$ и $\beta_2 \cdot service$.
# 
# В `statsmodels` переменная взаимодействия записывается через `:`:

# In[5]:


model2 = ols('salary_th ~ service + phd + sex +  sex:service', df).fit()
print(model2.summary())


# Уравнение модели:
# 
# $$
# salary = 73.6 + 18.5 \cdot sex + 0.17 \cdot service - 0.85 \cdot sex \cdot service + 1.54 \cdot phd
# $$
# 
# Запишем уравнения модели отдельно для женщин и мужчин.
# 
# `sex = 0`
# 
# $$
# salary_{female} = 73.6 + 0.17 \cdot service + 1.54 \cdot phd
# $$
# 
# `sex = 1`
# 
# $$
# salary_{male} = 73.6 + 18.5 + 0.17 \cdot service - 0.85 \cdot service + 1.54 \cdot phd
# $$
# 
# $$
# salary_{male} = 92.1 - 0.68 \cdot service + 1.54 \cdot phd
# $$

# Что мы здесь видим? Во-первых, заработная плата мужчин выше, чем заработная плата женщин. Во-вторых, опыт работы у женщин оказывает положительное влияние на заработную плату, а у мужчин – отрицательное. 
# 
# Можем записать уравнение, описывающее влияние переменной `service` в зависимости от пола, *предельный эффект* опыта работа на заработную плату. Предельный эффект какой-либо переменной определяется  как частная производная уравнения регрессии по этой переменной:
# 
# $$
# \frac{\partial salary}{\partial service} = 0.17 \cdot service - 0.85 \cdot sex 
# $$
# 
# И здесь мы опять видим, что влияние опыта работы у женщин более значимо отражается на заработной плате. 

# Посмотрим на значимость коэффициентов. Интересно, что при добавлении эффекта взаимодействия, соответствующий ему предиктор «оттянул» на себя значимость переменной `service`. Теперь сам по себе опыт работы не оказывает статистически значимого влияния на заработную плату, а вот с учётом пола преподавателя – оказывает (значим на 10% уровня значимости). Кроме того, значимое влияние на заработную плату оказывает пол преподавателя и число лет после защиты диссертации. 
# 
# Теперь давайте учтём в модели статус преподавателя, его должность. Посмотрим на уникальные значения переменной `rank`:

# In[6]:


df['rank'].unique()


# В отличие от переменной «пол» здесь уже не два значения, а три. Как быть? Включить в модель не саму переменную `rank`, а набор фиктивных переменных (дамми-переменных), которые являются бинарными.
# 
# Так, в нашем случае вместо `rank` будут созданы три дамми-переменных: `rank[AssocProf]`, `rank[T.AsstProf]` и `rank[T.Prof]`. Проиллюстрируем их смысл на небольшом фрагменте таблицы:
# 
#                 rank    rank[AssocProf]    rank[T.AsstProf]    rank[T.Prof]
#                 Prof                  0                   0               1  
#                 AssocProf             1                   0               0
#                 AsstProf              0                   1               0
# 
# 
# На самом деле, нет необходимости создавать дамми-переменные самим, Python опять это сделает самостоятельно, но важно понимать, что именно происходит. Посмотрим на модель:

# In[7]:


model3 = ols('salary_th ~ service + phd + sex + rank', df).fit()
print(model3.summary())


# Видно, что в модель были включены не три дамми-переменные для должности, а только две. Так будет всегда: чтобы избежать строгой мультиколлинеарности (два предиктора абсолютно скоррелированы, коэффициент корреляции равен 1), в модель будет включаться на одну дамми-переменную меньше. В модель не включена переменная `rank[AssocProf]`, следовательно, в качестве базового уровня выбрана должность доцента (`AssocProf`). С ней и будем сравнивать заработную плату других категорий преподавателей. Проинтерпретируем полученные результаты.
# 
# 1. Заработная плата профессора статистически значимо отличается от заработной платы доцента: при прочих равных зарплата профессора в среднем на 33 тысячи выше зарплаты доцента. 
# 
# 2. Заработная плата обычного преподавателя значимо отличается от заработной платы доцента: при прочих равных зарплата преподавателя в среднем на 14 тысяч ниже зарплаты доцента.
# 
# 3. Заработная плата преподавателей-мужчин в среднем на 5.5 тысяч выше, но если принимать во внимание должность преподавателя, эта разница не является статистически значимой.
# 
# 4. Если принимать во внимание должность преподавателя, то число лет после получения степени и опыт работы не оказывают значимого влияния на заработную плату. Должность «оттягивает» на себя значимость других предикторов и оказывает решающее влияние на размер заработной платы.
# 
# 5. Значение константы можно проинтерпретировать так: средняя заработная плата преподавателей равна 88 тысяч (без учёта пола преподавателя, числа лет опыта работы и других признаков, считаем все предикторы равными нулю). 

# Качество модели не очень высокое, но гораздо лучше, чем в предыдущих моделях (в двух предыдущих ноутбуках), наша последняя модель описывает 40% вариации заработной платы. Самостоятельно можете проверить модель на наличие мультиколлинеарности, гетероскедастичности и влиятельные наблюдения.