Отрицательная биномиальная регрессия: пошаговое руководство

Отрицательная биномиальная регрессия: пошаговое руководство

Плюс учебник Python по отрицательной биномиальной регрессии

В этой статье мы рассмотрим следующие темы:

1. Мы познакомимся с моделью отрицательной биномиальной регрессии (NB). Модель NB может быть невероятно полезной для прогнозирования данных на основе подсчета.
2. Мы рассмотрим пошаговое руководство по созданию, обучению и тестированию модели отрицательной биномиальной регрессии в Python с использованием класса GLM из statsmodels. .

Мотивация использования модели отрицательной биномиальной регрессии

В моей предыдущей статье мы познакомились с моделью регрессии Пуассона и увидели, как применить ее к данным на основе подсчета, таким как набор данных подсчета велосипедистов на Бруклинском мосту:

Мы также увидели, что модель регрессии Пуассона оказалась неадекватной для моделирования нашего набора данных о велосипедистах.

Хотя модель регрессии Пуассона делала предсказания, которые визуально удовлетворяли ...:

… Его результаты были статистически неудовлетворительными:

Низкая производительность модели объясняется тем, что данные не подчиняются критерию дисперсия = среднее, требуемому моделью регрессии Пуассона.

Этому довольно строгому критерию часто не удовлетворяют реальные данные. Часто дисперсия больше среднего, это свойство называется избыточной дисперсией, а иногда дисперсия меньше среднего, что называется недостаточной дисперсией. В таких случаях необходимо использовать регрессионную модель, которая не будет делать предположение равной дисперсии ienot Предположим, что дисперсия = среднее.

Модель регрессии Отрицательное биномиальное (NB) - это одна из таких моделей, которая не делает предположения о данных дисперсия = среднее.

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

Макет статьи

Статья оформлена следующим образом:

1. Мы познакомимся с набором реальных данных о счетчиках, который мы будем использовать в оставшейся части этой статьи.
2. Мы определим нашу цель регрессии на этом наборе данных.
3. Мы сформулируем стратегию регрессии, используя модель NB в качестве нашей регрессионной модели.
4. Мы настроим модель NB, обучим ее на наборе данных и сделаем некоторые прогнозы на тестовом наборе данных. Все это мы сделаем с помощью библиотеки statsmodels Python.
5. Наконец, мы проверим, действительно ли производительность модели NB превосходит производительность модели Пуассона.

Набор данных из реального мира

В следующей таблице указано количество велосипедистов, проезжающих по различным мостам Нью-Йорка. Учет проводился ежедневно с 1 апреля по 31 октября 2017 года.

Мы сосредоточим наш анализ на количестве велосипедистов, ежедневно пересекающих Бруклинский мост. Вот временной график подсчета велосипедистов на Бруклинском мосту:

Наша цель регресса

Наша цель регрессии - предсказать количество велосипедистов, пересекающих Бруклинский мост в любой день.

Наша стратегия регрессии

Учитывая значения набора регрессионных переменных для данного дня, мы будем использовать модель NB для прогнозирования количества велосипедистов на Бруклинском мосту в этот день.

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

y = вектор количества велосипедистов, наблюдаемых в дни с 1 по n.
Таким образом, y = [y_1, y_2, y_3,…, y_n].
y_i - количество велосипедистов в день i.

X = матрица предикторов, также называемых регрессорами, также известными как независимые переменные, также известные как переменные регрессии. Размер матрицы X равен (nxm), поскольку в данных n независимых наблюдений (строк) set, и каждая строка содержит значения m независимых переменных.

λ = вектор частоты событий. Вектор λ является основной характеристикой наборов данных на основе подсчета. λ - вектор размера (n x 1). Он содержит n коэффициентов [λ_0, λ_1, λ_2,…, λ_n] , соответствующих n наблюдаемых подсчетов в векторе отсчетов y. Предполагается, что коэффициент λ_i для наблюдения 'i' определяет фактическое наблюдаемое количество y_i в векторе отсчетов y. Столбец λ отсутствует во входных данных. Вместо этого вектор λ представляет собой выведенную переменную, которая вычисляется регрессионной моделью на этапе обучения.

Для данных подсчета велосипедистов каждое из значений λ_i определяется как количество велосипедистов, пересекающих мост за «единицу» времени в день i. Единица времени может быть 1 секунда, 1 час, 1 день, 1 неделя - независимо от того, за какой единичный временной интервал мы хотим измерить скорость. Предполагается, что этот показатель λ_i определяет наблюдаемое количество велосипедистов y_i в день i.

Следующий рисунок иллюстрирует эти определения на подмножестве нашего набора данных подсчета велосипедистов:

Алгоритм обучения модели отрицательной биномиальной регрессии подгонит наблюдаемые значения y к матрице регрессии X .

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

Вспомните, что модель отрицательной биномиальной регрессии не делает допущение дисперсия = среднее, как модель регрессии Пуассона.

Вместо этого модель NB требует от нас определения нового параметра α, который используется для выражения дисперсии в терминах среднего значения следующим образом :

На практике это уравнение принимает одну из двух часто встречающихся форм:

Когда p = 1:

В литературе по регрессии случай p = 1 упоминается как модель NB1. См. Cameron, A.C., and P.K. Триведи (1986), «Эконометрические модели, основанные на данных подсчета: сравнения и применения некоторых оценок

Когда p = 2:

Случай p = 2 называется моделью NB2.

Мы будем использовать модель NB2.

Библиотека Python statsmodels также поддерживает модель NB2 как часть предлагаемого ею класса обобщенной линейной модели.

Фактически, в пакете statsmodels.genmod.families.family есть целый класс, посвященный модели NB2:

class statsmodels.genmod.families.family.NegativeBinomial(link=None, alpha=1.0)

Обратите внимание, что значение по умолчанию alpha = 1, которое принимает этот класс, не всегда является правильным значением для всех наборов данных. Итак, как мы можем определить правильное значение α для нашего набора данных подсчета велосипедистов?

Нахождение правильного значения α

И снова господа Кэмерон и Триведи приходят нам на помощь. В своей книге Регрессионный анализ данных подсчета Кэмерон и Триведи предлагают хитрый способ вычисления α, используя технику, которую они называют вспомогательной регрессией МНК без константы. Уравнение регрессии, которое они рекомендуют , выглядит следующим образом:

Вы сразу можете увидеть взаимосвязь уравнения вспомогательной МНК с уравнением прямой регрессии: Y = B_1 * X + B_0.

Если вам интересно, уравнение для оценки α для модели NB1 выглядит следующим образом:

В оставшейся части этой статьи мы будем использовать модель NB2.

Мы сможем найти значение α, как только приспособим вспомогательное уравнение регрессии с помощью метода регрессии обыкновенных наименьших квадратов к нашему набору данных подсчетов . Мы увидим, как это сделать. сделай это как можно скорее.

Но как найти λ_i, который содержится в вспомогательном уравнении регрессии OLS ?

Чтобы найти λ_i , мы подгоняем модель регрессии Пуассона к нашему набору данных! Фактически, это дает нам полный вектор скорости λ = [λ_1, λ_2, λ_3,…, λ_n], соответствующий всем n наблюдениям. в наборе данных.

Теперь у нас есть все ингредиенты для стратегии регрессии NB2. Подведем итог.

Резюме стратегии регрессии NB2

• ШАГ 1. Подберите модель регрессии Пуассона к набору данных. Это даст нам вектор подходящих коэффициентов λ.
• ШАГ 2. Подберите регрессионную модель aux OLS к набору данных. Это даст нам значение α.
• ШАГ 3: Используйте α из ШАГА 2, чтобы подогнать модель регрессии NB2 к набору данных.
• ШАГ 4. Используйте подобранную модель NB2, чтобы делать прогнозы относительно ожидаемого числа в наборе тестовых данных.
• ШАГ 5. Проверьте соответствие модели NB2.

Теперь, когда наша стратегия регрессии намечена, давайте реализуем ее с помощью Python, Pandas и статистических моделей.

Как сделать отрицательную биномиальную регрессию в Python

Начнем с импорта всех необходимых пакетов.

import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

Затем создайте pandas DataFrame для набора данных счетчиков.

df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=[0], index_col=[0])

Мы добавим несколько производных регрессионных переменных в матрицу X.

ds = df.index.to_series()
df['MONTH'] = ds.dt.month
df['DAY_OF_WEEK'] = ds.dt.dayofweek
df['DAY'] = ds.dt.day

Мы не будем использовать переменную Date в качестве регрессора, поскольку она содержит абсолютное значение даты, но нам не нужно делать ничего особенного для удаления Date, поскольку она уже используется как индекс фрейма данных pandas. Поэтому он не будет доступен нам в матрице X.

Создадим наборы данных для обучения и тестирования.

mask = np.random.rand(len(df)) < 0.8
df_train = df[mask]
df_test = df[~mask]
print('Training data set length='+str(len(df_train)))
print('Testing data set length='+str(len(df_test)))

ШАГ 1. Теперь мы настроим и настроим модель регрессии Пуассона на обучающем наборе данных.

Задайте выражение регрессии в нотации patsy. Мы говорим Пэтси, что BB_COUNT - наша зависимая переменная, и она зависит от переменных регрессии: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T и PRECIP.

expr = """BB_COUNT ~ DAY  + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""

Настройте матрицы X и y для наборов данных для обучения и тестирования. Пэтси делает это очень просто.

y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')

Используя класс statsmodels GLM, обучите модель регрессии Пуассона на наборе обучающих данных.

poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()

На этом обучение модели регрессии Пуассона завершено. Чтобы увидеть результат обучения, вы можете распечатать его резюме.

print(poisson_training_results.summary())

Это распечатает следующее:

Наш реальный интерес заключается в векторе подобранных коэффициентов λ, полученном в процессе обучения. Этот вектор скорости содержится в параметре poisson_training_results.mu.

print(poisson_training_results.mu)
print(len(poisson_training_results.mu))

Следующий вывод показывает несколько первых и последних нескольких значений подобранного вектора λ:

[1920.39434028 2544.81549207 2651.79330653  743.45309242 1072.77132837 1892.9428398  2320.70412868 3269.73598361 3313.21764921 2915.25363322 2614.39482509 2594.44594144 2415.29471195 3181.91998369 2154.15471026 2268.25625592 1793.29625903 2535.31903414 2566.70835529 970.82159668
 2510.70775659 3016.19901465 2260.265944 3365.04650316 2695.6143122...
...2340.12964253 2568.40001641 2232.26752534 2604.97128321  145.92037793 2060.10442187 2518.70470296]

На этом завершается ШАГ1 подгонки модели регрессии Пуассона.

ШАГ 2. Теперь мы подгоним вспомогательную регрессионную модель OLS к набору данных и воспользуемся подобранной моделью, чтобы получить значение α.

Импортируйте пакет api.

import statsmodels.formula.api as smf

Добавьте вектор λ в качестве нового столбца под названием «BB_LAMBDA» во фрейм данных обучающего набора данных. Вспомните, что размеры λ ’ равны (n x 1). В нашем примере это будет (161 x 1). Также помните, что вектор λ доступен в poisson_training_results.mu:

df_train['BB_LAMBDA'] = poisson_training_results.mu

Затем давайте добавим производный столбец под названием «AUX_OLS_DEP» во фрейм данных pandas. В этом новом столбце будут храниться значения зависимой переменной регрессии OLS. Это левая часть уравнения регрессии OLS ниже:

df_train['AUX_OLS_DEP'] = df_train.apply(lambda x: ((x['BB_COUNT'] - x['BB_LAMBDA'])**2 - x['BB_LAMBDA']) / x['BB_LAMBDA'], axis=1)

В приведенном выше фрагменте кода часть, выделенная жирным шрифтом, является левой частью приведенного выше уравнения вспомогательного OLSR.

Давайте воспользуемся patsy, чтобы сформировать спецификацию модели для OLSR. Мы хотим сказать Пэтси, что AUX_OLS_DEP является зависимой переменной, и это объясняется BB_LAMBDA (который является вектором скорости λ). «-1» в конце выражения - это банальный синтаксис, обозначающий: не использовать перехват регрессии; т.е. просто поместите прямую линию, проходящую через начало координат, как предложили господа Кэмерон и Триведи.

ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""

Теперь мы готовы подобрать модель OLSR.

Сконфигурируйте и установите модель OLSR:

aux_olsr_results = smf.ols(ols_expr, df_train).fit()

Выведите параметры регрессии:

print(aux_olsr_results.params)

Вы увидите, как распечатывается следующий единственный коэффициент, соответствующий единственной переменной регрессии BB_LAMBDA. Этот коэффициент и есть α, который мы искали:

BB_LAMBDA    0.037343

Является ли α статистически значимым?

Теперь нам нужно ответить на очень важный вопрос. Является ли это значение α ( 0,037343 ) статистически значимым? Или его можно считать нулевым для всех практических целей?

Почему так важно это выяснить? Вспомните, что если α равно нулю, то следующее уравнение:

… Сокращается до Дисперсия = среднее. Это функция дисперсии модели регрессии Пуассона.

Если значение α статистически не значимо, то модель отрицательной биномиальной регрессии не может лучше подбирать набор обучающих данных, чем модель регрессии Пуассона.

Объект OLSResults содержит t-оценку коэффициента регрессии α. Давайте распечатаем:

aux_olsr_results.tvalues

Это распечатывает:

BB_LAMBDA    4.814096

Из калькулятора t-значения мы можем видеть, что критическое t-значение при уровне достоверности 99% (правый хвост) и степенях свободы = (161 наблюдение) - (1 параметр дисперсии α) = 160 составляет 2,34988. Это удобно меньше, чем t-статистика α, которая составляла 4,814096 . Мы заключаем, что,

α = 0,037343 статистически значимо.

На этом завершается ШАГ 2: определение α.

ШАГ 3. Мы передаем значение альфа, найденное на ШАГЕ 2, в класс statsmodels.genmod.families.family.NegativeBinomial и обучаем модель NB2 на наборе обучающих данных.

В статистических моделях это одноэтапная операция:

nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params[0])).fit()

Как и раньше, распечатаем сводку по обучению:

print(nb2_training_results.summary())

Что печатает следующее резюме:

ШАГ 4. Давайте сделаем некоторые прогнозы, используя нашу обученную модель NB2.

Прогнозирование снова представляет собой одноэтапную процедуру в статистических моделях:

nb2_predictions = nb2_training_results.get_prediction(X_test)

Распечатываем прогнозы:

predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)

Ниже приведены первые несколько строк вывода:

Давайте также отобразим прогнозируемое количество в сравнении с фактическим количеством для тестовых данных.

predicted_counts=predictions_summary_frame['mean']
actual_counts = y_test['BB_COUNT']
fig = plt.figure()
fig.suptitle('Predicted versus actual bicyclist counts on the Brooklyn bridge')
predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted counts')
actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual counts')
plt.legend(handles=[predicted, actual])
plt.show()

Вот результат:

Не плохо! Модель NB2, похоже, более или менее отслеживает тенденцию в подсчете количества велосипедов. И так же, как и в случае с характеристиками регрессионной модели Пуассона, в некоторых случаях ее прогнозы не соответствуют фактическим значениям.

Вот полный исходный код Python для обучения модели отрицательной биномиальной регрессии и тестирования ее прогнозов:

Последний вопрос, который стоит перед нами:

Статистически наша регрессионная модель NB2 работает лучше, чем регрессионная модель Пуассона?

Давайте разберемся.

ШАГ 5: Измерение степени соответствия модели NB2

С точки зрения согласия, в сводке по обучению модели NB2 есть три интересных момента. На рисунке ниже они заключены в красный прямоугольник. Мы рассмотрим каждый из них, начиная с логарифма вероятности.

Давайте сначала сравним сводку обучения модели NB2 с моделью регрессии Пуассона на том же наборе данных:

Первый статистический показатель, на который следует обратить внимание, - это значение Log-Likelihood. Максимальное логарифмическое правдоподобие было получено с помощью метода оценки максимального правдоподобия (MLE), который выполнялся статистическими моделями во время обучения моделей Пуассона и NB2. Метод MLE используется для фиксации значений всех коэффициентов модели на некоторых оптимальных значениях, которые максимизируют вероятность увидеть вектор отсчетов y в наборе обучающих данных. Чтобы узнать больше о MLE и о том, как он используется в обучении моделей, обратитесь к моей статье Модель регрессии Пуассона.

Тест отношения правдоподобия (LR)

Тест отношение правдоподобия используется для сравнения того, насколько две модели соответствуют данным.

Статистика теста LR просто отрицательна, вдвое больше, чем разница в подобранных логарифмических вероятностях двух моделей.

В нашем случае логарифмическое правдоподобие для NB2 составляет -1383,2, а для модели регрессии Пуассона - -12616. Таким образом, статистика теста LR равна 2 * (12616–1383,2) = 22465,6. Это значение значительно превышает критическое значение χ2 (1) на уровне значимости 1%, которое составляет 5,412.

Согласно тесту LR, обученная модель регрессии NB2 продемонстрировала гораздо лучшее соответствие набору данных велосипедистов по сравнению с моделью регрессии Пуассона.

Теперь давайте сравним степень согласия регрессионной модели NB2 в абсолютном выражении.

Девианс и статистика хи-квадрат Пирсона

Сообщенные значения девианса и хи-квадрат Пирсона для модели NB2 составляют 330,99 и 310 соответственно. Чтобы количественно определить степень согласия на некотором уровне достоверности, скажем, 95% (p = 0,05), мы ищем значение в таблице χ2 для p = 0,05 и степеней свободы. остатков = 165. Мы сравниваем это значение хи-квадрат с наблюдаемой статистикой - в данном случае это отклонение или значение хи-квадрат Пирсона, указанное в GLMResults. Мы находим, что при p = 0,05 и DF Residuals = 165 значение хи-квадрат из стандартной таблицы хи-квадрат составляет 195,973, что меньше заявленных статистических значений 330,99 и 310. Следовательно, согласно этому тесту, регрессия NB2 Модель, несмотря на то, что демонстрирует гораздо лучшее соответствие, чем модель регрессии Пуассона, по-прежнему неоптимальна. Мы могли добиться большего.

Заключение и следующие шаги

Модели Пуассона и модели отрицательной биномиальной регрессии используются для моделирования наборов данных на основе подсчетов. Обе модели дают следующие результаты:

• Объяснимо
• Сопоставимый
• Оправданный
• Пригодный к употреблению

Обе модели опираются на сильную и хорошо изученную статистическую теорию.

Для выполнения регрессии на наборах данных, основанных на подсчете, хорошая стратегия - начать с модели регрессии Пуассона, а затем посмотреть, сможете ли вы получить лучшие результаты, используя модель отрицательной биномиальной регрессии.

Если ни Пуассон, ни NB2 не подходят для вашего набора данных, рассмотрите возможность использования более продвинутых методов, таких как:

1. Сложные варианты модели регрессии Пуассона, такие как Модель с нулевым раздутием.
2. Модель препятствий
3. Регрессионная модель на основе случайного леса
4. Модель регрессии на основе нейронной сети Long-Short Term Memory (LSTM).

Связанное чтение

Спасибо за внимание! Если вам понравилась эта статья, подпишитесь на меня в Sachin Date, чтобы получать советы, инструкции и советы по программированию по темам, посвященным регрессии, анализу временных рядов и прогнозированию.