Продолжим работать с базой по заработной плате: загрузим необходимые библиотеки и сам файл, переименуем столбцы с точкой в названии:
from statsmodels.formula.api import ols
import pandas as pd
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()
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 |
Попробуем учесть в модели пол преподавателя. В нашей таблице переменная sex
является текстовой (два значения female
и male
), поэтому не очень ясно, как включать её в модель. Обычно в таких случаях переменную превращают в бинарную. Но в данном случае этого делать не нужно, Python выполнит преобразования автоматически: упорядочит два значения по алфавиту, первому присвоит значение $0$, второму – значение $1$.
В предыдущих моделях заработная плата у нас измерялась как есть, в условных единицах, но для качества модели и для интерпретации будет лучше, если разброс значений зависимой переменной будет меньше. Давайте добавим переменную salary_th
– заработную плату в тысячах, а затем будем её использовать в качестве зависимой переменной.
df['salary_th'] = df['salary'] / 1000
Построим модель, где предикторами являются переменные service
, phd
, и sex
:
model1 = ols('salary_th ~ service + phd + sex', df).fit()
print(model1.summary())
OLS Regression Results ============================================================================== Dep. Variable: salary_th R-squared: 0.195 Model: OLS Adj. R-squared: 0.189 Method: Least Squares F-statistic: 31.75 Date: Fri, 04 Jan 2019 Prob (F-statistic): 2.13e-18 Time: 23:42:59 Log-Likelihood: -1873.8 No. Observations: 397 AIC: 3756. Df Residuals: 393 BIC: 3772. Df Model: 3 Covariance Type: nonrobust =============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------- Intercept 82.8759 4.801 17.264 0.000 73.438 92.314 sex[T.Male] 8.4571 4.656 1.816 0.070 -0.697 17.611 service -0.6498 0.254 -2.558 0.011 -1.149 -0.150 phd 1.5528 0.256 6.062 0.000 1.049 2.056 ============================================================================== Omnibus: 14.548 Durbin-Watson: 1.888 Prob(Omnibus): 0.001 Jarque-Bera (JB): 15.448 Skew: 0.425 Prob(JB): 0.000442 Kurtosis: 3.460 Cond. No. 156. ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
В модели добавился предиктор sex[T.Male]
, коэффициент при котором показывает разницу в заработной плате мужчин и женщин. Так как в качестве базовой категории (значения $0$, категории, с которой мы сравниваем остальные) используется значение female
, коэффициент при sex
показывает, на сколько больше/меньше заработная плата у мужчин по сравнению с женщинами.
Запишем уравнение с учётом коэффициентов модели:
Известно, что для преподавателей женского пола $sex=0$, для преподавателей мужского пола $sex=1$. Учтём этот факт и запишем уравнения отдельно для каждого пола:
Видно, что в среднем, при прочих равных, заработная плата у мужчин на 8.46 тысяч выше, чем заработная плата у женщин. Это и есть коэффициент при sex[T.Male]
в таблице. Теперь пойдём дальше и попробуем учесть не только влияние пола преподавателя, но и тот факт, что влияние опыта работы у мужчин и женщин неодинаково. Добавим в модель эффект взаимодействия – предиктор, представляющий собой произведение двух переменных, в данном случае произведение sex
и service
. В общем виде уравнение модели выглядит так:
Важно: чтобы избежать смещения оценок коэффициентов и проблемы пропущенных переменных, переменные которые входят в эффект взаимодействия, должны быть включены в модель по отдельности. Так, здесь мы не можем убрать слагаемые $\beta_1 \cdot sex$ и $\beta_2 \cdot service$.
В statsmodels
переменная взаимодействия записывается через :
:
model2 = ols('salary_th ~ service + phd + sex + sex:service', df).fit()
print(model2.summary())
OLS Regression Results ============================================================================== Dep. Variable: salary_th R-squared: 0.201 Model: OLS Adj. R-squared: 0.193 Method: Least Squares F-statistic: 24.60 Date: Fri, 04 Jan 2019 Prob (F-statistic): 3.47e-18 Time: 23:42:59 Log-Likelihood: -1872.4 No. Observations: 397 AIC: 3755. Df Residuals: 392 BIC: 3775. Df Model: 4 Covariance Type: nonrobust ======================================================================================= coef std err t P>|t| [0.025 0.975] --------------------------------------------------------------------------------------- Intercept 73.5908 7.385 9.965 0.000 59.072 88.110 sex[T.Male] 18.5159 7.659 2.418 0.016 3.458 33.574 service 0.1696 0.557 0.305 0.761 -0.925 1.265 sex[T.Male]:service -0.8473 0.513 -1.652 0.099 -1.856 0.161 phd 1.5412 0.256 6.028 0.000 1.039 2.044 ============================================================================== Omnibus: 15.188 Durbin-Watson: 1.888 Prob(Omnibus): 0.001 Jarque-Bera (JB): 16.321 Skew: 0.430 Prob(JB): 0.000286 Kurtosis: 3.496 Cond. No. 302. ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Уравнение модели:
$$ salary = 73.6 + 18.5 \cdot sex + 0.17 \cdot service - 0.85 \cdot sex \cdot service + 1.54 \cdot phd $$Запишем уравнения модели отдельно для женщин и мужчин.
sex = 0
sex = 1
Что мы здесь видим? Во-первых, заработная плата мужчин выше, чем заработная плата женщин. Во-вторых, опыт работы у женщин оказывает положительное влияние на заработную плату, а у мужчин – отрицательное.
Можем записать уравнение, описывающее влияние переменной service
в зависимости от пола, предельный эффект опыта работа на заработную плату. Предельный эффект какой-либо переменной определяется как частная производная уравнения регрессии по этой переменной:
И здесь мы опять видим, что влияние опыта работы у женщин более значимо отражается на заработной плате.
Посмотрим на значимость коэффициентов. Интересно, что при добавлении эффекта взаимодействия, соответствующий ему предиктор «оттянул» на себя значимость переменной service
. Теперь сам по себе опыт работы не оказывает статистически значимого влияния на заработную плату, а вот с учётом пола преподавателя – оказывает (значим на 10% уровня значимости). Кроме того, значимое влияние на заработную плату оказывает пол преподавателя и число лет после защиты диссертации.
Теперь давайте учтём в модели статус преподавателя, его должность. Посмотрим на уникальные значения переменной rank
:
df['rank'].unique()
array(['Prof', 'AsstProf', 'AssocProf'], dtype=object)
В отличие от переменной «пол» здесь уже не два значения, а три. Как быть? Включить в модель не саму переменную 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 опять это сделает самостоятельно, но важно понимать, что именно происходит. Посмотрим на модель:
model3 = ols('salary_th ~ service + phd + sex + rank', df).fit()
print(model3.summary())
OLS Regression Results ============================================================================== Dep. Variable: salary_th R-squared: 0.402 Model: OLS Adj. R-squared: 0.394 Method: Least Squares F-statistic: 52.51 Date: Fri, 04 Jan 2019 Prob (F-statistic): 1.29e-41 Time: 23:42:59 Log-Likelihood: -1814.9 No. Observations: 397 AIC: 3642. Df Residuals: 391 BIC: 3666. Df Model: 5 Covariance Type: nonrobust ==================================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------------ Intercept 89.5964 4.891 18.318 0.000 79.980 99.212 sex[T.Male] 5.4991 4.035 1.363 0.174 -2.433 13.431 rank[T.AsstProf] -13.8866 4.333 -3.205 0.001 -22.406 -5.367 rank[T.Prof] 33.0534 3.701 8.932 0.000 25.778 40.329 service -0.3738 0.221 -1.693 0.091 -0.808 0.060 phd 0.2658 0.248 1.072 0.284 -0.222 0.753 ============================================================================== Omnibus: 39.686 Durbin-Watson: 1.769 Prob(Omnibus): 0.000 Jarque-Bera (JB): 60.293 Skew: 0.668 Prob(JB): 8.08e-14 Kurtosis: 4.363 Cond. No. 178. ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Видно, что в модель были включены не три дамми-переменные для должности, а только две. Так будет всегда: чтобы избежать строгой мультиколлинеарности (два предиктора абсолютно скоррелированы, коэффициент корреляции равен 1), в модель будет включаться на одну дамми-переменную меньше. В модель не включена переменная rank[AssocProf]
, следовательно, в качестве базового уровня выбрана должность доцента (AssocProf
). С ней и будем сравнивать заработную плату других категорий преподавателей. Проинтерпретируем полученные результаты.
Заработная плата профессора статистически значимо отличается от заработной платы доцента: при прочих равных зарплата профессора в среднем на 33 тысячи выше зарплаты доцента.
Заработная плата обычного преподавателя значимо отличается от заработной платы доцента: при прочих равных зарплата преподавателя в среднем на 14 тысяч ниже зарплаты доцента.
Заработная плата преподавателей-мужчин в среднем на 5.5 тысяч выше, но если принимать во внимание должность преподавателя, эта разница не является статистически значимой.
Если принимать во внимание должность преподавателя, то число лет после получения степени и опыт работы не оказывают значимого влияния на заработную плату. Должность «оттягивает» на себя значимость других предикторов и оказывает решающее влияние на размер заработной платы.
Значение константы можно проинтерпретировать так: средняя заработная плата преподавателей равна 88 тысяч (без учёта пола преподавателя, числа лет опыта работы и других признаков, считаем все предикторы равными нулю).
Качество модели не очень высокое, но гораздо лучше, чем в предыдущих моделях (в двух предыдущих ноутбуках), наша последняя модель описывает 40% вариации заработной платы. Самостоятельно можете проверить модель на наличие мультиколлинеарности, гетероскедастичности и влиятельные наблюдения.