In [1]:
import statsmodels
import scipy as sc
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.stats.api as sms
from statsmodels.graphics.regressionplots import plot_leverage_resid2
import matplotlib.pyplot as plt
In [2]:
%pylab inline
По 1260 опрошенным имеются следующие данные:
Требуется оценить влияние внешней привлекательности на уровень заработка с учётом всех остальных факторов.
Hamermesh D.S., Biddle J.E. (1994) Beauty and the Labor Market, American Economic Review, 84, 1174–1194.
Данные:
In [3]:
raw = pd.read_csv("beauty.csv", sep=";", index_col=False)
raw.head()
Out[3]:
Посмотрим на матрицу диаграмм рассеяния по количественным признакам:
In [4]:
pd.tools.plotting.scatter_matrix(raw[['wage', 'exper', 'educ', 'looks']], alpha=0.2,
figsize=(15, 15), diagonal='hist')
pylab.show()
Оценим сбалансированность выборки по категориальным признакам:
In [5]:
print raw.union.value_counts()
print raw.goodhlth.value_counts()
print raw.black.value_counts()
print raw.female.value_counts()
print raw.married.value_counts()
print raw.service.value_counts()
У каждого признака все значения встречаются достаточно много раз, так что всё в порядке.
In [6]:
data = raw
Посмотрим на распределение целевого признака — уровня заработной платы:
In [7]:
plt.figure(figsize(16,7))
plt.subplot(121)
data['wage'].plot.hist()
plt.xlabel('Wage', fontsize=14)
plt.subplot(122)
np.log(data['wage']).plot.hist()
plt.xlabel('Log wage', fontsize=14)
pylab.show()
Один человек в выборке получает 77.72\$ в час, остальные — меньше 45\$; удалим этого человека, чтобы регрессия на него не перенастроилась.
In [8]:
data = data[data['wage'] < 77]
Посмотрим на распределение оценок привлекательности:
In [9]:
plt.figure(figsize(8,7))
data.groupby('looks')['looks'].agg(lambda x: len(x)).plot(kind='bar', width=0.9)
plt.xticks(rotation=0)
plt.xlabel('Looks', fontsize=14)
pylab.show()
В группах looks=1 и looks=5 слишком мало наблюдений. Превратим признак looks в категориальный и закодируем с помощью фиктивных переменных:
In [10]:
data['belowavg'] = data['looks'].apply(lambda x : 1 if x < 3 else 0)
data['aboveavg'] = data['looks'].apply(lambda x : 1 if x > 3 else 0)
data.drop('looks', axis=1, inplace=True)
Данные теперь:
In [11]:
data.head()
Out[11]:
Построим линейную модель по всем признакам.
In [12]:
m1 = smf.ols('wage ~ exper + union + goodhlth + black + female + married +'\
'service + educ + belowavg + aboveavg',
data=data)
fitted = m1.fit()
print fitted.summary()
Посмотрим на распределение остатков:
In [13]:
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
pylab.show()
Оно скошенное, как и исходный признак. В таких ситуациях часто помогает перейти от регрессии исходного признака к регрессии его логарифма.
In [14]:
m2 = smf.ols('np.log(wage) ~ exper + union + goodhlth + black + female + married +'\
'service + educ + belowavg + aboveavg', data=data)
fitted = m2.fit()
print fitted.summary()
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
pylab.show()
Теперь стало лучше. Посмотрим теперь на зависимость остатков от непрерывных признаков:
In [15]:
plt.figure(figsize(16,7))
plt.subplot(121)
scatter(data['educ'],fitted.resid)
plt.xlabel('Education', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.subplot(122)
scatter(data['exper'],fitted.resid)
plt.xlabel('Experience', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
pylab.show()
На втором графике видна квадратичная зависимость остатков от опыта работы. Попробуем добавить к признакам квадрат опыта работы, чтобы учесть этот эффект.
In [16]:
m3 = smf.ols('np.log(wage) ~ exper + np.power(exper,2) + union + goodhlth + black + female +'\
'married + service + educ + belowavg + aboveavg', data=data)
fitted = m3.fit()
print fitted.summary()
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
plt.figure(figsize(16,5))
plt.subplot(131)
scatter(data['educ'],fitted.resid)
plt.xlabel('Education', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.subplot(132)
scatter(data['exper'],fitted.resid)
plt.xlabel('Experience', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
plt.subplot(133)
scatter(data['exper']**2,fitted.resid)
plt.xlabel('Experience^2', fontsize=14)
plt.ylabel('Residuals', fontsize=14)
pylab.show()
Используем критерий Бройша-Пагана для проверки гомоскедастичности ошибок:
In [17]:
print 'Breusch-Pagan test: p=%f' % sms.het_breushpagan(fitted.resid, fitted.model.exog)[1]
Ошибки гетероскедастичны, значит, значимость признаков может определяться неверно. Сделаем поправку Уайта:
In [18]:
m4 = smf.ols('np.log(wage) ~ exper + np.power(exper,2) + union + goodhlth + black + female +'\
'married + service + educ + belowavg + aboveavg', data=data)
fitted = m4.fit(cov_type='HC1')
print fitted.summary()
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
pylab.show()
В предыдущей модели незначимы: цвет кожи, здоровье, семейное положение. Удалим их. Индикатор привлекательности выше среднего тоже незначим, но удалять его не будем, потому что это одна из переменных, по которым на нужно в конце ответить на вопрос.
In [19]:
m5 = smf.ols('np.log(wage) ~ exper + np.power(exper,2) + union + female + service + educ +'\
'belowavg + aboveavg', data=data)
fitted = m5.fit(cov_type='HC1')
print fitted.summary()
plt.figure(figsize(16,7))
plt.subplot(121)
sc.stats.probplot(fitted.resid, dist="norm", plot=pylab)
plt.subplot(122)
np.log(fitted.resid).plot.hist()
plt.xlabel('Residuals', fontsize=14)
pylab.show()
Посмотрим, не стала ли модель от удаления трёх признаков значимо хуже, с помощью критерия Фишера:
In [20]:
print "F=%f, p=%f, k1=%f" % m4.fit().compare_f_test(m5.fit())
Не стала.
Проверим, нет ли наблюдений, которые слишком сильно влияют на регрессионное уравнение:
In [21]:
plt.figure(figsize(8,7))
plot_leverage_resid2(fitted)
pylab.show()
In [22]:
data.loc[[1122]]
Out[22]:
In [23]:
data.loc[[269]]
Out[23]:
Итоговая модель объясняет 40% вариации логарифма отклика.
In [24]:
plt.figure(figsize(16,7))
plt.subplot(121)
scatter(data['wage'],np.exp(fitted.fittedvalues))
plt.xlabel('Wage', fontsize=14)
plt.ylabel('Exponentiated predictions', fontsize=14)
plt.xlim([0,50])
plt.subplot(122)
scatter(np.log(data['wage']),fitted.fittedvalues)
plt.xlabel('Log wage', fontsize=14)
plt.ylabel('Predictions', fontsize=14)
plt.xlim([0,4])
pylab.show()
При интересующих нас факторах привлекательности стоят коэффициенты -0.1307 (ниже среднего) и -0.0010 (выше среднего).
Поскольку регрессия делалась на логарифм отклика, интерпретировать их можно как прирост в процентах. С учётом дополнительных факторов представители генеральной совокупности, из которой взята выборка, получают в среднем: