Алла Тамбовцева, НИУ ВШЭ
Продолжим работать с таблицей, содержащие данные по преподавателям и их заработной плате:
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
Импортируем функцию ols
:
from statsmodels.formula.api import ols
Построим модель, где зависимой переменной является заработная плата, а независимыми – число лет опыта работа и число лет после получения степени. Учтём обе независимые переменные в модели (константа всегда добавляется по умолчанию):
model1 = ols('salary ~ service + phd', df).fit()
Посмотрим на результаты:
print(model1.summary())
OLS Regression Results ============================================================================== Dep. Variable: salary R-squared: 0.188 Model: OLS Adj. R-squared: 0.184 Method: Least Squares F-statistic: 45.71 Date: Mon, 31 Dec 2018 Prob (F-statistic): 1.40e-18 Time: 21:44:51 Log-Likelihood: -4617.9 No. Observations: 397 AIC: 9242. Df Residuals: 394 BIC: 9254. Df Model: 2 Covariance Type: nonrobust ============================================================================== coef std err t P>|t| [0.025 0.975] ------------------------------------------------------------------------------ Intercept 8.991e+04 2843.560 31.620 0.000 8.43e+04 9.55e+04 service -629.1014 254.469 -2.472 0.014 -1129.389 -128.814 phd 1562.8889 256.820 6.086 0.000 1057.981 2067.797 ============================================================================== Omnibus: 14.927 Durbin-Watson: 1.867 Prob(Omnibus): 0.001 Jarque-Bera (JB): 15.947 Skew: 0.429 Prob(JB): 0.000344 Kurtosis: 3.478 Cond. No. 69.6 ============================================================================== Warnings: [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
По значениям в столбце P>|t|
видно, что оба коэффициента при предикторах являются статистически значимыми на уровне значимости 5%, гипотеза о равенстве коэффициентов нулю отвергается. То же показывает и p-value для F-statistic
. Запишем уравнение модели на основе полученных результатов (коэффициенты округлены):
В отличие от модели парной регрессии, здесь коэффициент при service
уже отрицательный. Получается, что число лет опыта работы негативно сказывается на заработной плате! С одной стороны, это кажется контринтуитивным, но с другой стороны, вполне может оказаться, что другой фактор, включённый в модель, phd
является более значимым и в сочетании с service
«оттягивает» значимость на себя. Давайте внимательно посмотрим на уравнение модели.
В проигрыше те преподаватели, которые имеют большой опыт работы, но при этом недавно получили степень. Ситуация достаточно логичная; можно предположить, что если заработная плата зависит от научных достижений, то сотрудники, работающие давно и при этом не стремящиеся как можно скорее получить научную степень, защитить диссертацию и писать статьи, получают меньшую заработную плату. Здесь мы рискуем задуматься и отдалиться от регрессионной модели, поэтому давайте вернёмся к таблице.
Качество модели, если мы оцениваем его по значению коэффициента детерминации $R^2$, не кажется высоким, особенно, учитывая, что это довольно стандартная модель, поэтому $R^2$ в данном случае можно свободно интерпретировать как долю дисперсии заработной платы, объяснённую моделью (предикторами «опыт работы» и «число лет после получения степени»). Самое время проверить другие показатели качества модели. Из уже имеющихся в таблице показателей, в частности, из p-value для значения F-статистики, можем заключить, что модель нужна: гипотеза о равенстве всех коэффициентов регресии нулю отвергается.
Какие аспекты в приведённой выше таблицы не учтены? Во-первых, проверка модели на наличие потенциальных проблем, возникающих во множественной регрессии.
Первая проблема – мультиколлинеарность, сильная линейная связь между предикторами. К чему эта проблема может привести? К неверным выводам о значимости коэффициентов регрессии, которые следуют из слишком «раздутых» стандартных ошибок коэффициентов, имеющих место при мультиколлинеарности. Как распознать мультиколлинеарность в модели? Рассмотрим некоторые признаки.
Пожалуй, самый однозначный результат можно получить, воспользовавшись последним способом. Однако посмотрим и на первые два. Второй способ сразу исключим, поскольку в нашем случае высокого $R^2$ и незначимых предикторов не наблюдается. Для проверки первым способом посчитаем коэффициент корреляции Пирсона между двумя предикторами в модели:
Для вычисления VIF я позволила себе воспользоваться готовой функцией со stackoverflow, чтобы вывести значения VIF для всех переменных сразу. Встроенная в statsmodels
функция variance_inflation_factor()
расчитывает значения VIF для конкретных переменных и при этом исключает из модели константу, что отражается на результатах (они отличаются от привычных значений VIF, которые обычно используются и интерпретируются в эконометрике).
# взято с https://stackoverflow.com/questions/42658379/variance-inflation-factor-in-python
from statsmodels.regression.linear_model import OLS
from statsmodels.tools.tools import add_constant
def variance_inflation_factors(exog_df):
'''
Parameters
----------
exog_df : dataframe, (nobs, k_vars)
design matrix with all explanatory variables, as for example used in
regression.
Returns
-------
vif : Series
variance inflation factors
'''
exog_df = add_constant(exog_df)
vifs = pd.Series(
[1 / (1. - OLS(exog_df[col].values,
exog_df.loc[:, exog_df.columns != col].values).fit().rsquared)
for col in exog_df],
index=exog_df.columns,
name='VIF'
)
return vifs
variance_inflation_factors(df[['service', 'phd']])
const 4.289177 service 5.795810 phd 5.795810 Name: VIF, dtype: float64
По значениям VIF видно, что они превышают принятое пороговое значение $5$, но несильно. Можно считать, что данная модель не страдает от мультиколлинеарности. Тем не менее, давайте обсудим, как с этой проблемой бороться. Итак, методы борьбы с мультиколлинеарностью:
Вторая возможная проблема – гетероскедастичность – непостоянство дисперсии остатков. Одно из допущений линейной модели (одно из условий Гаусса-Маркова) – постоянство дисперсии остатков/ошибок модели, то есть $D(\varepsilon_i) = const$ для всех $i$. Чтобы распознать гетероскедастичность в модели, можно воспользоваться специальным тестом. Однако обычно рекомендуется думать об источниках гетероскедастичности заранее, если их по содержательным соображениям быть не должно, то не стоит слепо доверять результатам тестов на гетероскедастичность, так как они могут работать нестабильно, зависеть от числа наблюдений в выборке и структуры данных.
Рассмотрим тест Бройша-Пагана, один из популярных тестов на гетероскедастичность. Нулевая гипотеза, проверяемая с помощью этого теста – дисперсия остатков модели постоянна:
$H_0: D(\varepsilon_i) = D(\varepsilon_j) = \sigma^2$ для всех $i \ne j$ (гомоскедастичность).
Применим тест к нашей модели.
from statsmodels.stats.diagnostic import het_breuschpagan
Функция het_breuschpagan
из statsmodels
принимает на вход остатки модели и набор независимых переменных.
het_breuschpagan(model1.resid, df[['phd', 'service']])
(149.44370609921543, 2.293790755993904e-34, 119.22594044982796, 3.0866291683237355e-41)
Этот тест выводит значения F-статистики и p-value для каждой переменной. Значения p-value в нашем случае близки к нулю, поэтому чисто формально гипотезу о гомоскедастичности следует отвергнуть. Как справляться с гетероскедастичностью в модели? Обычно выбирают другой тип стандартных ошибок коэффициентов, устойчивый к гетероскедастичности (heteroskedasticity consistent или точнее Huber–White standard errors или White standard errors). Такие ошибки ещё часто называют робастными, однако робастность – это более общее свойство, относящееся к различным типам скорректированных ошибок (устойчивые к гетероскедастичности, к автокорреляции, к возмущениям в модели, к наличию различий по группам наблюдений).
Есть несколько подтипов ошибок, устойчивых к гетероскедастичности, на практике часто используют 2 типа: hc0 (для больших выборок) и hc3 (для маленьких выборок).
Учесть тип ошибок модели можно при оценке модели, указав их в качестве аргумента методв .fit()
.
model_hc = ols('salary ~ service + phd', df).fit(cov_type='hc3')
print(model_hc.summary())
OLS Regression Results ============================================================================== Dep. Variable: salary R-squared: 0.188 Model: OLS Adj. R-squared: 0.184 Method: Least Squares F-statistic: 38.02 Date: Mon, 31 Dec 2018 Prob (F-statistic): 7.96e-16 Time: 21:44:51 Log-Likelihood: -4617.9 No. Observations: 397 AIC: 9242. Df Residuals: 394 BIC: 9254. Df Model: 2 Covariance Type: hc3 ============================================================================== coef std err z P>|z| [0.025 0.975] ------------------------------------------------------------------------------ Intercept 8.991e+04 2440.680 36.839 0.000 8.51e+04 9.47e+04 service -629.1014 309.075 -2.035 0.042 -1234.877 -23.326 phd 1562.8889 284.489 5.494 0.000 1005.302 2120.476 ============================================================================== Omnibus: 14.927 Durbin-Watson: 1.867 Prob(Omnibus): 0.001 Jarque-Bera (JB): 15.947 Skew: 0.429 Prob(JB): 0.000344 Kurtosis: 3.478 Cond. No. 69.6 ============================================================================== Warnings: [1] Standard Errors are heteroscedasticity robust (hc3)
Теперь мы видим, что тип ошибок выбран другой, и это указано в параметрах модели (Covariance Type
). Значения коэффициентов не отличаются от значений в предыдущей модели, но вот стандартные ошибки изменились, и, соответственно, значимость тоже (и распределение статистики критерия для проверки гипотезы о значимости коэффициент уже не Стьюдента, а стандартное нормальное). К выбору типа ошибок стоит относиться внимательно, поскольку изменение типа стандартных ошибок ведёт к изменению p-value и выводов относительно нулевой гипотезы.
Ещё одна проблема, которая может отразиться на качестве модели и на выводах относительно значимости коэффициентов, – это наличие влиятельных наблюдений. Про нетипичные наблюдения (выбросы, outliers), вы уже, скорее всего знаете, это просто значения, сильно отклоняющиеся от среднего и медианы. Формально границы типичных значений в выборке определяются как $[Q1-1.5\Delta; Q3+1.5\Delta]$, где $Q_1$ и $Q_3$ – нижний и верхний квартили, а $\Delta$ – межквартильный размах ($Q3-Q1$). Влиятельные наблюдения – наблюдения в данных, которые влияют на расположение регрессионной прямой или гиперплоскости. Если рассмотрим случай парной регрессии, то при удалении влиятельного наблюдения наклон прямой может радикально измениться, коэффициент может даже сменить знак. Это плохо, потому что устойчивая модель, которая даёт стабильные результаты, не должна зависеть от таких локальных изменений. Влиятельные наблюдения часто удаляют из данных, а затем оценивают модель заново.
Один из распространённых методов выявления влиятельных наблюдений – оценить меру Кука (расстояние Кука, Cook's distance) для каждого наблюдения. Пороговое значение меры Кука, превышение которого свидетельствует о влиятельности наблюдения, равно $\frac{4}{n-p-1}$, где $n$ – число наблюдений, $p$ – число предикторов. В нашем случае $\frac{4}{397-2-1} \approx 0.01$. Посмотрим на значения меры Кука для наблюдений в нашей модели:
from statsmodels.stats import outliers_influence
influence = outliers_influence.OLSInfluence(model_hc)
cookd = influence.cooks_distance[0]
cookd
array([1.64522825e-03, 4.52821510e-03, 7.67725627e-04, 2.02234416e-03, 1.17868559e-03, 7.52446855e-06, 4.55331291e-03, 1.61554561e-03, 1.36027966e-04, 9.89635971e-04, 4.85722212e-04, 1.11238278e-03, 7.99906728e-04, 9.11209131e-04, 3.59013882e-05, 3.11851393e-04, 8.55457751e-05, 1.78278233e-03, 3.79453042e-04, 2.63923502e-04, 1.73787282e-03, 1.41803049e-03, 1.99192957e-03, 6.61574087e-06, 1.62937397e-03, 5.10666541e-04, 8.50739797e-05, 5.91093845e-04, 3.69253049e-03, 3.95579366e-04, 1.18950521e-03, 1.09938318e-03, 2.77360944e-04, 7.54966568e-04, 7.54966568e-04, 1.42389639e-03, 3.30179480e-03, 3.93764330e-04, 3.54699179e-04, 1.96625028e-05, 5.70147560e-03, 6.50747773e-04, 6.02095482e-03, 4.78629948e-02, 3.46789441e-04, 4.97417002e-05, 6.24085034e-03, 1.65294493e-03, 1.52975151e-03, 1.86592263e-03, 2.67178142e-04, 1.10646886e-04, 2.52500332e-03, 1.17213009e-05, 8.82618650e-06, 1.57735851e-03, 3.64837068e-05, 1.99322625e-04, 1.60669311e-07, 1.70391632e-03, 1.95279432e-04, 1.27356336e-03, 3.59272256e-04, 3.01816089e-05, 2.43535644e-03, 6.50052773e-06, 5.54844412e-04, 1.55318298e-03, 6.66233309e-05, 3.19973904e-03, 7.27267789e-04, 1.43536119e-03, 1.24277687e-03, 2.96640122e-03, 4.55024234e-05, 1.94997263e-03, 6.53957635e-03, 7.61829318e-03, 2.32902172e-04, 5.66670625e-04, 8.95948865e-04, 1.47821135e-03, 1.52539535e-03, 2.52517225e-04, 7.62838263e-04, 1.49547513e-03, 3.18776225e-03, 4.97460861e-05, 7.39623668e-03, 1.35317447e-04, 6.43584745e-05, 3.48537970e-05, 4.41938299e-05, 7.01160580e-03, 3.00856456e-04, 5.13606947e-04, 2.38147807e-04, 1.07471907e-03, 4.02619607e-03, 6.35856533e-03, 3.69239728e-04, 8.06482098e-05, 5.56902353e-03, 2.87028750e-04, 1.38303342e-03, 9.71757924e-05, 8.09837841e-04, 7.31927542e-04, 1.31152120e-03, 1.06887594e-05, 2.06069265e-06, 9.65597116e-04, 1.73286948e-03, 1.62767332e-03, 6.03337842e-05, 4.59814952e-05, 2.03300000e-03, 3.91884871e-04, 1.82919470e-03, 1.62811637e-03, 3.39609486e-04, 6.03356120e-05, 3.89112516e-04, 3.60668912e-03, 4.20627489e-04, 3.59640479e-02, 2.53768812e-03, 1.69963155e-03, 1.35965992e-04, 1.68093778e-03, 7.24419625e-04, 5.10829493e-02, 1.20817219e-03, 1.73286948e-03, 4.66016930e-03, 9.64416770e-04, 7.34167528e-06, 5.56596004e-06, 1.62640652e-03, 2.10730796e-03, 7.77952335e-05, 9.97160144e-04, 6.28467271e-05, 3.14676617e-05, 1.35043983e-05, 2.09353076e-05, 1.04566601e-05, 4.00361221e-03, 6.83133503e-04, 2.34340497e-06, 1.00625982e-03, 1.04566601e-05, 1.14848214e-04, 5.21769033e-06, 6.01816163e-05, 1.75316829e-04, 1.46594806e-03, 5.35194448e-05, 3.89944628e-08, 2.88623945e-03, 2.23176207e-05, 4.61410708e-03, 2.44882309e-03, 3.14676617e-05, 3.18329590e-05, 5.47482649e-04, 8.12838098e-03, 1.43381468e-03, 1.73624020e-05, 5.44168836e-03, 4.07818279e-05, 3.79888295e-03, 7.19394188e-03, 1.28371327e-03, 1.84348036e-04, 5.87878180e-05, 7.11415486e-05, 2.26246697e-05, 2.16656926e-03, 2.08576391e-06, 4.58884528e-03, 2.89909799e-03, 2.07568524e-05, 1.39620501e-03, 2.17483433e-04, 2.81519532e-04, 6.25818791e-08, 3.21329952e-05, 2.28028014e-04, 1.64796100e-03, 1.45327557e-02, 1.18669066e-03, 3.17200765e-04, 9.02215414e-04, 1.99253806e-02, 4.81728599e-04, 3.45808860e-06, 1.04566601e-05, 1.31559571e-02, 3.01248719e-03, 3.45808860e-06, 2.74490955e-04, 3.75886153e-03, 4.43500310e-03, 1.13760708e-02, 5.91602321e-03, 4.09517851e-03, 3.51586140e-04, 1.95311402e-04, 5.84339010e-03, 3.90563134e-05, 1.01321760e-03, 7.90324274e-04, 7.96006332e-05, 5.28812903e-04, 1.99053252e-03, 2.33956613e-03, 3.99988757e-04, 9.96658779e-06, 4.90074719e-04, 7.18166603e-03, 2.50166453e-03, 1.00936017e-06, 3.64251763e-06, 4.95127684e-03, 3.73770416e-04, 3.39897818e-03, 2.14189537e-03, 5.60660392e-04, 1.96062950e-04, 1.37098209e-03, 2.79305251e-03, 1.26306529e-03, 2.06827227e-03, 2.49346812e-03, 4.57356784e-03, 1.17276505e-06, 3.44029871e-03, 1.78352907e-02, 1.58079426e-03, 2.31521743e-03, 2.76190972e-05, 2.09055107e-03, 2.04948070e-05, 1.74325502e-03, 5.25373257e-04, 1.60075012e-03, 1.31586303e-04, 1.24314833e-04, 6.92754561e-02, 1.35605414e-03, 6.73480628e-04, 1.28344704e-05, 1.06125240e-03, 1.90080488e-03, 7.91322505e-04, 1.76538229e-03, 3.83369795e-03, 2.05863433e-03, 1.82778623e-03, 7.30216457e-03, 3.83166722e-03, 1.07042885e-06, 1.42555184e-04, 2.29561015e-03, 1.35861355e-04, 1.10485146e-03, 2.12397542e-06, 1.25810223e-02, 8.54000773e-06, 7.53201542e-04, 4.34355679e-02, 1.75064595e-03, 1.70838810e-03, 1.27310188e-03, 2.88251947e-04, 8.81761078e-03, 3.01353250e-03, 9.71123537e-05, 7.10258447e-03, 3.70458187e-04, 2.41563971e-03, 5.69213871e-02, 2.84916601e-03, 2.62522560e-04, 2.33613389e-02, 1.59669036e-06, 2.60334669e-04, 1.92019338e-03, 1.62776738e-03, 2.73176657e-02, 3.96516971e-03, 4.19528958e-02, 9.92549030e-06, 2.16003796e-04, 3.69581858e-03, 7.64308238e-04, 2.21475045e-03, 2.76722807e-02, 1.99312414e-02, 5.37006857e-03, 6.77336938e-05, 1.05014501e-02, 4.85299774e-05, 5.81487501e-04, 2.73012756e-04, 1.56322998e-03, 5.75486550e-06, 1.86592044e-03, 3.61228895e-03, 2.44722401e-03, 6.35675738e-06, 1.94089837e-03, 1.74969603e-03, 4.30516841e-03, 4.77752672e-04, 1.92444294e-03, 2.74800183e-02, 1.83921258e-03, 1.56833146e-03, 1.25127203e-04, 2.04585448e-05, 9.99483886e-04, 3.38951372e-03, 6.27382528e-03, 6.03952900e-04, 1.05956957e-04, 2.95342638e-03, 2.53344609e-11, 1.10221446e-03, 8.01388208e-02, 2.87241684e-06, 7.05738342e-05, 1.14511857e-03, 5.22816581e-04, 4.12874082e-03, 3.89776621e-03, 3.69064835e-03, 3.38626102e-05, 1.29253224e-05, 1.59559667e-05, 6.97504392e-04, 4.95479064e-04, 5.23987423e-03, 4.88274615e-03, 1.29561802e-06, 1.50249487e-04, 1.04842517e-05, 7.27158678e-04, 2.42657839e-03, 1.49380168e-02, 4.31645120e-03, 2.20484885e-03, 9.72199502e-04, 9.97020468e-04, 1.03545941e-03, 1.74110675e-02, 1.49253166e-03, 9.85968219e-04, 1.60855509e-03, 4.78678804e-04, 7.07592023e-05, 1.15047058e-03, 1.06773216e-03, 3.26577401e-02, 4.35348166e-03, 1.06903652e-04, 4.02354018e-05, 8.69657896e-05, 3.89583899e-04, 1.28951586e-03, 9.54006613e-06, 5.65788862e-05, 4.21207802e-04, 4.58781395e-04, 3.25131882e-04, 1.47414775e-03, 1.22585292e-03, 2.71554336e-03, 8.53084396e-06, 1.57898851e-03, 3.98427265e-03, 3.90359456e-04, 4.11248965e-03, 5.18978462e-05, 7.20146203e-05, 8.11778805e-04, 1.21437740e-03, 1.79203016e-04, 1.69958038e-02, 8.45476630e-03, 1.95500632e-03, 8.26512145e-04, 1.95747811e-03, 1.22635169e-02, 1.28236875e-03, 9.06330816e-04])
Посмотрим, если ли влиятельные наблюдения:
cookd[cookd > 0.01]
array([0.04786299, 0.03596405, 0.05108295, 0.01453276, 0.01992538, 0.01315596, 0.01137607, 0.01783529, 0.06927546, 0.01258102, 0.04343557, 0.05692139, 0.02336134, 0.02731767, 0.0419529 , 0.02767228, 0.01993124, 0.01050145, 0.02748002, 0.08013882, 0.01493802, 0.01741107, 0.03265774, 0.0169958 , 0.01226352])
Такие есть и, в принципе, можно определить индексы таких наблюдений и удалить их. Однако давайте посмотрим на график, который обычно используется для визуализации нетипичных и влиятельных наблюдений, и примем окончательное решение. По оси абсцисс в таком графике указано значение влиятельности (leverage), а по оси ординат – стьюдентизированные остатки модели, которые показывают нетипичность наблюдений. Тип графика – пузырьковая диаграмма (bubble plot), диаграмма рассеяния с точками различного размера. Величина точки изменяется пропорционально мере Кука, чем больше мера Кука, тем больше точка.
from statsmodels.api import graphics
from matplotlib import pyplot as plt
fig, ax = plt.subplots(figsize=(12, 8))
fig = graphics.influence_plot(model1, ax=ax, criterion="cooks")
График не очень красивый для демонстрации, но из него можно вынести много полезной информации. Обычно из массива данных для построения регрессии удаляют не все значения с высокой влиятельностью, а те, которые являются и нетипичными, и влиятельными одновременно. На графике выше это те точки, которые лежат в правом верхнем углу. Решение об удалении влиятельных наблюдений обычно субъективно, зависит от исследователя. В данном случае очень нетипичными и влиятельными наблюдениями являются наблюдения с индексами 271, 290, 292 и 330. В дальнейшем их желательно удалить и переоценить модель.