Пользователь
0,0
рейтинг
17 декабря 2013 в 08:36

Разработка → Пример решения задачи множественной регрессии с помощью Python

Введение


Добрый день, уважаемые читатели.
В прошлых статьях, на практических примерах, мной были показаны способы решения задач классификации (задача кредитного скоринга) и основ анализа текстовой информации (задача о паспортах). Сегодня же мне бы хотелось коснуться другого класса задач, а именно восстановления регрессии. Задачи данного класса, как правило, используются при прогнозировании.
Для примера решения задачи прогнозирования, я взял набор данных Energy efficiency из крупнейшего репозитория UCI. В качестве инструментов по традиции будем использовать Python c аналитическими пакетами pandas и scikit-learn.

Описание набора данных и постановка задачи


Дан набор данных, котором описаны следующие атрибуты помещения:
Поле Описание Тип
X1 Относительная компактность FLOAT
X2 Площадь FLOAT
X3 Площадь стены FLOAT
X4 Площадь потолка FLOAT
X5 Общая высота FLOAT
X6 Ориентация INT
X7 Площадь остекления FLOAT
X8 Распределенная площадь остекления INT
y1 Нагрузка при обогреве FLOAT
y2 Нагрузка при охлаждении FLOAT

В нем — характеристики помещения на основании которых будет проводиться анализ, а — значения нагрузки, которые надо спрогнозировать.

Предварительный анализ данных


Для начала загрузим наши данные и посмотрим на них:

from pandas import read_csv, DataFrame
from sklearn.neighbors import KNeighborsRegressor
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
from sklearn.cross_validation import train_test_split

dataset = read_csv('EnergyEfficiency/ENB2012_data.csv',';')
dataset.head()

X1 X2 X3 X4 X5 X6 X7 X8 Y1 Y2
0 0.98 514.5 294.0 110.25 7 2 0 0 15.55 21.33
1 0.98 514.5 294.0 110.25 7 3 0 0 15.55 21.33
2 0.98 514.5 294.0 110.25 7 4 0 0 15.55 21.33
3 0.98 514.5 294.0 110.25 7 5 0 0 15.55 21.33
4 0.90 563.5 318.5 122.50 7 2 0 0 20.84 28.28

Теперь давайте посмотрим не связаны ли между собой какие-либо атрибуты. Сделать это можно рассчитав коэффициенты корреляции для всех столбцов. Как это сделать было описано в предыдущей статье:

dataset.corr()

X1 X2 X3 X4 X5 X6 X7 X8 Y1 Y2
X1 1.000000e+00 -9.919015e-01 -2.037817e-01 -8.688234e-01 8.277473e-01 0.000000 1.283986e-17 1.764620e-17 0.622272 0.634339
X2 -9.919015e-01 1.000000e+00 1.955016e-01 8.807195e-01 -8.581477e-01 0.000000 1.318356e-16 -3.558613e-16 -0.658120 -0.672999
X3 -2.037817e-01 1.955016e-01 1.000000e+00 -2.923165e-01 2.809757e-01 0.000000 -7.969726e-19 0.000000e+00 0.455671 0.427117
X4 -8.688234e-01 8.807195e-01 -2.923165e-01 1.000000e+00 -9.725122e-01 0.000000 -1.381805e-16 -1.079129e-16 -0.861828 -0.862547
X5 8.277473e-01 -8.581477e-01 2.809757e-01 -9.725122e-01 1.000000e+00 0.000000 1.861418e-18 0.000000e+00 0.889431 0.895785
X6 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 1.000000 0.000000e+00 0.000000e+00 -0.002587 0.014290
X7 1.283986e-17 1.318356e-16 -7.969726e-19 -1.381805e-16 1.861418e-18 0.000000 1.000000e+00 2.129642e-01 0.269841 0.207505
X8 1.764620e-17 -3.558613e-16 0.000000e+00 -1.079129e-16 0.000000e+00 0.000000 2.129642e-01 1.000000e+00 0.087368 0.050525
Y1 6.222722e-01 -6.581202e-01 4.556712e-01 -8.618283e-01 8.894307e-01 -0.002587 2.698410e-01 8.736759e-02 1.000000 0.975862
Y2 6.343391e-01 -6.729989e-01 4.271170e-01 -8.625466e-01 8.957852e-01 0.014290 2.075050e-01 5.052512e-02 0.975862 1.000000

Как можно заметить из нашей матрицы, коррелируют между собой следующие столбы (Значение коэффициента корреляции больше 95%):
  • y1 --> y2
  • x1 --> x2
  • x4 --> x5

Теперь давайте выберем, какие столбцы их наших пар мы можем убрать из нашей выборки. Для этого, в каждой паре, выберем столбцы, которые в большей степени оказывают влияние на прогнозные значения Y1 и Y2 и оставим их, а остальные удалим.
Как можно заметить и матрицы с коэффициентами корреляции на y1,y2 больше значения оказывают X2 и X5, нежели X1 и X4, таким образом мы можем последние столбцы мы можем удалить.

dataset = dataset.drop(['X1','X4'], axis=1)
dataset.head()

Помимо этого, можно заметить, что поля Y1 и Y2 очень тесно коррелируют между собой. Но, т. к. нам надо спрогнозировать оба значения мы их оставляем «как есть».

Выбор модели


Отделим от нашей выборки прогнозные значения:

trg = dataset[['Y1','Y2']]
trn = dataset.drop(['Y1','Y2'], axis=1)

После обработки данных можно перейти к построению модели. Для построения модели будем использовать следующие методы:

Теорию о данным методам можно почитать в курсе лекций К.В.Воронцова по машинному обучению.
Оценку будем производить с помощью коэффициента детерминации (R-квадрат). Данный коэффициент определяется следующим образом:



, где image — условная дисперсия зависимой величины у по фактору х.
Коэффициент принимает значение на промежутке и чем он ближе к 1 тем сильнее зависимость.
Ну что же теперь можно перейти непосредственно к построению модели и выбору модели. Давайте поместим все наши модели в один список для удобства дальнейшего анализа:

	models = [LinearRegression(), # метод наименьших квадратов
	          RandomForestRegressor(n_estimators=100, max_features ='sqrt'), # случайный лес
	          KNeighborsRegressor(n_neighbors=6), # метод ближайших соседей
	          SVR(kernel='linear'), # метод опорных векторов с линейным ядром
	          LogisticRegression() # логистическая регрессия
	          ]

Итак модели готовы, теперь мы разобьем наши исходные данные на 2 подвыборки: тестовую и обучающую. Кто читал мои предыдущие статьи знает, что сделать это можно с помощью функции train_test_split() из пакета scikit-learn:

Xtrn, Xtest, Ytrn, Ytest = train_test_split(trn, trg, test_size=0.4)

Теперь, т. к. нам надо спрогнозировать 2 параметра , надо построить регрессию для каждого из них. Кроме этого, для дальнейшего анализа, можно записать полученные результаты во временный DataFrame. Сделать это можно так:

	#создаем временные структуры
	TestModels = DataFrame()
	tmp = {}
	#для каждой модели из списка
	for model in models:
	    #получаем имя модели
	    m = str(model)
	    tmp['Model'] = m[:m.index('(')]    
	    #для каждого столбцам результирующего набора
	    for i in  xrange(Ytrn.shape[1]):
	        #обучаем модель
	        model.fit(Xtrn, Ytrn[:,i]) 
	        #вычисляем коэффициент детерминации
	        tmp['R2_Y%s'%str(i+1)] = r2_score(Ytest[:,0], model.predict(Xtest))
	    #записываем данные и итоговый DataFrame
	    TestModels = TestModels.append([tmp])
	#делаем индекс по названию модели
	TestModels.set_index('Model', inplace=True)

Как можно заметить из кода выше, для расчета коэффициента используется функция r2_score().
Итак, данные для анализа получены. Давайте теперь построим графики и посмотрим какая модель показала лучший результат:

fig, axes = plt.subplots(ncols=2, figsize=(10,4))
TestModels.R2_Y1.plot(ax=axes[0], kind='bar', title='R2_Y1')
TestModels.R2_Y2.plot(ax=axes[1], kind='bar', color='green', title='R2_Y2')

image

Анализ результатов и выводы


Из графиков, приведенных выше, можно сделать вывод, что лучше других с задачей справился метод RandomForest (случайный лес). Его коэффициенты детерминации выше остальных по обоим переменным:
ля дальнейшего анализа давайте заново обучим нашу модель:

model = models[1]
model.fit(Xtrn, Ytrn)

При внимательном рассмотрении, может возникнуть вопрос, почему в предыдущий раз и делили зависимую выборку Ytrn на переменные(по столбцам), а теперь мы этого не делаем.
Дело в том, что некоторые методы, такие как RandomForestRegressor, может работать с несколькими прогнозируемыми переменными, а другие (например SVR) могут работать только с одной переменной. Поэтому на при предыдущем обучении мы использовали разбиение по столбцам, чтобы избежать ошибки в процессе построения некоторых моделей.
Выбрать модель это, конечно же, хорошо, но еще неплохо бы обладать информацией, как каждый фактор влиет на прогнозное значение. Для этого у модели есть свойство feature_importances_.
С помощью него, можно посмотреть вес каждого фактора в итоговой моделей:

model.feature_importances_

array([ 0.40717901, 0.11394948, 0.34984766, 0.00751686, 0.09158358,
0.02992342])


В нашем случае видно, что больше всего на нагрузку при обогреве и охлаждении влияют общая высота и площадь. Их общий вклад в прогнозной модели около 72%.
Также необходимо отметить, что по вышеуказанной схеме можно посмотреть влияние каждого фактора отдельно на обогрев и отдельно на охлаждение, но т. к. эти факторы у нас очень тесно коррелируют между собой (), мы сделали общий вывод по ним обоим который и был написан выше.

Заключение


В статье я постарался показать основные этапы при регрессионном анализе данных с помощью Python и аналитческих пакетов pandas и scikit-learn.
Необходимо отметить, что набор данных специально выбирался таким образом чтобы быть максимально формализованым и первичная обработка входных данных была бы минимальна. На мой взгляд статья будет полезна тем, кто только начинает свой путь в анализе данных, а также тем кто имеет хорошую теоретическую базу, но выбирает инструментарий для работы.
@kuznetsovin
карма
55,2
рейтинг 0,0
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама

Самое читаемое Разработка

Комментарии (20)

  • 0
    Можете посоветовать, что почитать на тему прогнозирования временных рядов?
    • 0
      Я планирую написать следующую статью по основам анализа временных рядов на python+statsmodels. При ее написании буду использовать лекции по Анализу временных рядов ВШЭ, ссылку на них приложу в статье. Планирую выложить статью в конце этой недели.
  • 0
    Это работает только с числами?
    Как поступить если во входных данных есть категорийный столбец(например цвет)?
    • 0
      вы интересуетесь как это вообще сделать, или именно в Питоне? Если в целом, то этому посвящен раздел эконометрики «Limited and Categorical Predictor Variables». Чаще всего, категорийные переменные кодируются как {0;1} дамми-переменные, например переменная female={0 если мужчина; 0; если ребенок; 1 — если женщина}, и так для всех переменных — мужчина, женщина, ребенок.
      Можно использовать Multiple Categorical кодирование (1-красный, 2-синий, и тп). Еще как вариант, можно закодировать относительную частоту вхождения категорий в выборку.
      Только и оценивать модель уже надо будет другими способами.
      • 0
        Интересует общий подход. До {0;1} я дошел, но иногда бывают столбцы с 100 000 вариантов состояний. Нормально ли там будет работать Multiple Categorical?

        И еще есть вопрос: есть набор полей допустим
        возраст|пол{0;1}|рост || итог{0;1}
        25|0|187|0
        26|0|192|1

        Нужно по входным данным определять выходное значение.

        Проблема в том что есть множество строк с одинаковыми значениями полей. ( Ну т.е в магазин 100 раз зашел один человек и 1 раз совершил покупку). Хотелось бы чтобы модель в таком случае выдавала 0.01 а не 0.
        Но множество строк по одному набору тоже есть не всегда(т.е допустим человек зашел в магазин 1 раз и сразу совершил покупку, не факт что в следующий раз он тоже совершит, хотелось бы как-то брать агрегированную информацию с заходов похожих людей).
        Не подскажете есть какие-то статьи на эту тему, или ключевые слова для поиска?

        • 0
          По первому вопросу, честно, не знаю. Надо порыться в литературе, больше 4-5 категорий я не встречал в приложениях. Как вариант — если у вас n категорий, создаете n-1 дамми (0;1) и их используете в регрессии.
          Насчет второго вопроса я не совсем уверен что понял что именно вы спрашиваете.
          Насколько я понимаю, у вас бинарный результат (купил / не купил) и разношерстные вводные данные. В таком случае — Limited dependent variable модели. Классические варианты — логит или пробит.
  • 0
    Можно его закодировать через hotonecoder или label encoder про данные функции я писал в предыдущих своих статьях, например здесь
  • +1
    У меня тут несколько замечаний по статье.
    1.Теперь давайте посмотрим не связаны ли между собой какие-либо атрибуты. Сделать это можно рассчитав коэффициенты корреляции для всех столбцов.
    Смотреть на корреляцию каждый раз можно, но толку от этого практически ноль. Как говорили древние – «Correlation does not imply causation». Помните картинку с сомалийскими пиратами и температурой на планете? Так вот, там практически стопроцентная корелляция. Но если так уж хочется корелляцию посмотреть, то возьмите тогда уж Granger causality test, там хоть используются предыдущие лаги чтобы определить, влияет ли Х на У и наоборот. Да и F-тест поавторитетнее обычного коэффициента корелляции будет.

    2. Как можно заметить и матрицы с коэффициентами корреляции на y1,y2 больше значения оказывают X2 и X5, нежели X1 и X4, таким образом мы можем последние столбцы мы можем удалить.
    Facepalm. Как вы можете так просто удалять данные только потому что у них коэффициент кореляции низкий? Если Х1 влияет на У в два раза слабее чем Х2, это не повод выбрасывать Х1, он все равно значимый. Для определения этого и существует t-test, и основанная на нем p-value для тестирования нулевой гипотезы, что переменная Х1 не является значимой. Только прогнав регрессию, вы сможете выбрасывать переменные. Кстати, насчет выбрасывания данных я еще в прошлый раз хотел вам на это указать, когда вы сказали что «Это связано с тем, что у одно заемщика быть несколько кредитов и по каждому из них в разных бюро моте быть разная информация» и сделали .drop_duplicates для заемщиков, оставив 50000 записей из 280942 и удалив таким образом 80% информации.

    3. Если уж заниматься первичной обработкой данных, вам следовало обратить внимание, что они очень отличаются по величине. Так, у Х1 все значения меньше единицы, а Х2 больше 500. Хорошо бы было привести их к более-менее одинаковой шкале, трансформировав через логарифм. Во-вторых, если вы посмотрите на гистограммы переменных, то они очень далеки от нормальных, что является еще одним поводом сделать логарифмическую трансформацию. И запускать регрессию не на абсолютных значениях, а в логарифмически масштабированной шкале – log(Х1), log (X2), так званая log-log или log-linear регрессия. Правда тогда нужно не забыть трансформировать полученные коэффициенты обратно.

    4. R-квадрат для оценки результатов можно использовать очень условно. Вы не проверили есть ли у вас в данных:
    a. Мультиколлинеарность коэффициентов. В случае наличия, ваши коэффициенты ошибочны. Именно высокий R2 и есть первым признаком мультиколлинеарности. Кстати, когда я вижу R2 больше 90%, это сразу повод усомниться в валидности вашей модели.
    b. Heteroscedasticity (не знаю как это перевести на русский, но в общем дисперсия ошибок регрессии получается непостоянной. Результат — коэффициенты неэффективны и ошибочны. Зачем вам высокий R2, если найденные параметры неправильные?
    c. Автокорреляция или серийная корреляция ошибок регрессии. Результат – почти все тесты коэффициентов регрессии будут ошибочны.
    В любом из этих случаев, результаты всех моделей которые вы описали будут ошибочными, и R-квадрат как показатель оценки бесполезен. Да и построенный прогноз далек от реалий. Для выявления всех проблем существует множество тестов и измененных моделей.

    5. И еще много чего можно обсудить в этой и прошлых статьях. Я не пытаюсь придираться, просто если вы хотите серьезно заниматься статистикой, начните с хороших книг, а не с питона.
    Данная статья может быть интересна только как пример использования питона для статистических задач, но по-моему, это как микроскопом гвозди забивать. Если вам интересна статистика или эконометрика, возьмитесь за Стату, Матлаб, R, EViews.
    • 0
      Кстати, только что быстро прикинул Метод квадратов (OLS) для ваших данных по основным проблемам регрессии из пункта 4.

      Итого:
      Heteroskedasticity Test: Breusch-Pagan-Godfrey: F-test = 6,761, probability =0.0000
      Итого, нулевая гипотеза про отсутвие Heteroscedasticity может быть отброшена.

      Breusch-Godfrey Serial Correlation LM Test: F-test = 2,759, probability 0.0000
      Ошибки регрессии кореллированы. Горе.

      Jarque–Bera normality test:Test statistic: 20.49, Probability = 0.00003
      О нормальности ошибок можно забыть. Модель ошибочна.

      В общем, все плохо, и это только навскидку.
    • 0
      Данная статья и задумывалась, как пример использования python в качестве альтернативы указанным выше пакетам. По поводу теории буду признателен, если порекомендуете книги.
      • +1
        Если есть база в статистике, то можно начать с Гуджарати — хорошая вводная книга из набора «на пальцах» по многим направлениям в эконометрике. Я тут пока искал ссылки, даже случайно нашел pdf-скан на предыдущую версию.
        Также неплохие пособия для начала — Вербек, Wooldridge

        Немного сложнее уровень, на уровне магистерки — Джонстон и ДиНардо. Тоже обо всем по немногу.

        На уровне аспирантуры — Библия Эконометриста. Но ее можно в основном как справочник использовать.

        Ну а дальше уже куда специализация и интерес заведет.
        По panel data — советую Baltagi, Wooldbridge
        Time series — Hamilton, Bisgaard
        Ну и еще с полторадесятка направлений есть (финансовая эконометрика, Censored data, Bayesian econometrics, Truncated and censored models, Discrete regression and Qualitative choice models, Duration analysis и тд и тп), в каждом свои знаковые авторы есть. Большинство книг можно найти в интернете отсканированными. Предыдущие версии обычно не намного хуже новинок, так как сама база как и в математике остается постоянной.
        Если вас интересует именно практическая сторона, то поищите книги «заточенные» под какой-нибудь пакет. Например R, Stata, EViews
        • 0
          А если нет базы в статистике (а теорвер был давно и забылся за ненадобностью), какие бы книги посоветовали для начала?
          • 0
            увы, статистику я очень давно проходил. Не помню ни авторов, ни книг. Да и половина материала была лекторская. Погуглите форумы, там вроде можно найти отзывы как по русскоязычной, так и иностранной литературе.
    • 0
      2. Как можно заметить и матрицы с коэффициентами корреляции на y1,y2 больше значения оказывают X2 и X5, нежели X1 и X4, таким образом мы можем последние столбцы мы можем удалить.
      Этот шаг как раз и подразумевал избавление от мультиколениарности т.к. между X2 и X1, и X5 и X4 соответствующие коэффициенты парной корреляции r > 90%. Разве не так?
      • 0
        Да, все правильно. Если рассматривать корелляцию между Х1 и Х2, то коэффициент корелляции между ними очень высок (-0.99), значит можно подразумевать мультиколлинеарность. Поэтому теоретически можно выбросить один из них – Х1 или Х2. То же относительно Х4 и Х5 – если они коррелируют, то скрипач одна переменная не нужна. Но даже в этом случае их нельзя выбросить просто так. Помните – про пиратов и температуру? Объяснение этому — Spurious relationship. Посмотрите интересные примеры, например, как кореллируют Индекс экспорта в США и продолжительность жизни мужчин в Австралии. Корреляция еще не означает что переменные взаимосвязаны и взаимозаменяемы. Надо поочередно выбросить по одной из переменной и посмотреть как изменятся результаты регрессии. Если коэффициенты значительно изменились, значит выбросили не напрасно. Или же с использовать что-то из серии Farrar-Glauber Test, Variance Inflation Factor test, и иже с ними.
        Автор в целом был на правильном пути. Он собирался выбросить переменные из кореллирующих пар (X1, X2), (X4,X5). Но при этом говоря, что большее влияние оказывают X2 и X5, поэтому их и оставим. Ок, corr(X1,Y1)=0.62, corr(X2,Y1)=-0.66. Разница в влиянии – минимальна, и трудно сказать почему X2 лучше чем X1, особенно при наличии в модели еще шести переменных. Опять же – просто так взять и выбросить – нельзя. Только проверять.
        • 0
          Тогда такой вопрос, если я правильно понимаю, если мы применим метод главных компонент, то в принципе мы с большой долей вероятности избавимся от мультиколлениарности, так?
          • 0
            да, это поможет. в principal component analysis вы вместо измерения влияния отдельных переменных используете их средневзвешенное влияние.
            Но проще выбросить ненужную переменную)
            • 0
              Я правильно понимаю, что чем ниже p-value, тем больше вероятность неправильности нулевой гипотезы? или же наоборот?
              • +1
                да, все так.
                p-value указывает на вероятность нулевой гипотезы. Чем меньше p-value, тем более высокая вероятность что H0 неверна и ее можно отбросить.
                Например, в обычной регрессии смотрим на p-values коэффициентов. H0 в данном случае утверждает что коэффициент при переменной равен нулю, т.е. наша переменная не имеет никакого влияния. Тогда низкая p-value (обычно меньше 10%, 5%, или 1%, в зависимости от выбранного уровня significancy), означает что вы можете отбросить нулевую гипотезу, и переменная нужна и важна)
                Или, например, в тесте Breusch-Godfrey на серийную корелляцию, H0 — это гипотеза об отсутствии корелляции. В таком случае высокая p-val означает что гипотезу можно принять, и у вас есть корелляция.
                • 0
                  om nom nom nom

Только зарегистрированные пользователи могут оставлять комментарии. Войдите, пожалуйста.