Иллюстрированное руководство по модели регрессии Пуассона с нулевым раздутием

Плюс учебник Python по обучению модели ZIP на наборах данных с лишними нулями.

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

Наборы данных подсчета - это наборы, в которых зависимой переменной является событие, например:

  • Количество транспортных средств, пересекающих перекресток в час.
  • Количество посещений скорой помощи каждый месяц
  • Количество поданных претензий по страхованию автотранспортных средств в год
  • Количество дефектов, обнаруженных в печатной плате серийного производства.

Многие явления в реальном мире приводят к почти всегда нулевым счетам. Например:

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

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

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

Например, если вы предполагаете, что явление подчиняется следующему процессу Пуассона (5), вы ожидаете увидеть нулевые значения не более 0,67% времени:

Если вы наблюдаете нулевые отсчеты гораздо чаще, значит, набор данных содержит избыток нулей.

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

Так что же делать моделисту, когда он сталкивается с такими данными с лишними нулями?

Модель регрессии Пуассона с нулевым раздутием

К счастью, есть способ изменить стандартную модель подсчета, такую ​​как Пуассон или Отрицательный бином, чтобы учесть наличие лишних нулей. На самом деле, есть как минимум два способа сделать это. Один метод известен как модель препятствий, а второй метод известен как модель с нулевым накачиванием.

В этой статье мы более подробно рассмотрим модель регрессии с нулевым раздутием. В частности, мы сосредоточимся на модели регрессии Пуассона с нулевым раздутием, часто называемой моделью ZIP.

Структура модели ZIP

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

Представьте себе набор данных, содержащий n выборок и p переменных регрессии для каждой выборки. Следовательно, переменные регрессии X могут быть представлены матрицей размера (nxp) и каждой строкой x_i в матрице X - это вектор размером (1 xp), соответствующий значению зависимой переменной y_i :

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

  1. Функция Пуассона P M ass F (PMF), обозначаемая как P (y_i = k), используемая для расчета вероятность наблюдения k событий в любом единичном интервале при средней частоте событий λ событий / единицу времени.
  2. Функция связи, которая используется для выражения средней скорости λ как функции переменных регрессии X.

Это показано на рисунке ниже:

Обычно мы предполагаем, что существует некий базовый процесс, производящий наблюдаемые подсчеты согласно Пуассону PMF: P (y_i = k).

Интуиция, лежащая в основе модели Пуассона с нулевым надуванием, заключается в том, что существует второй базовый процесс, который определяет, является ли счетчик нулевым или ненулевым. Как только счетчик определен как ненулевой, обычный процесс Пуассона берет на себя определение его фактического ненулевого значения на основе PMF процесса Пуассона.

Таким образом, регрессионная модель ZIP состоит из трех частей:

  1. PMF P (y_i = 0), который используется для вычисления вероятности наблюдения нулевого счета.
  2. Второй PMF P (y_i = k), который используется для вычисления вероятности наблюдения k событий, учитывая, что k ›0.
  3. Функция связи, которая используется для выражения средней скорости λ как функции переменных регрессии X.

Это показано на следующем рисунке:

Как и раньше, y_i - это случайная величина, обозначающая наблюдаемое количество, соответствующее строке регрессионных переменных x_i = [x_i1, x_i2, x_i3,…, x_ip].

ϕ_i - это мера доли лишних нулей, соответствующих i-й строке ( y_i , x_i) в набор данных; .

Знакомство ϕ_i

Простой способ понять ϕ_i заключается в следующем:

Представьте, что вы делаете 1000 наблюдений y_i, каждое с одинаковым сочетанием значений регрессионных переменных x_i = [x_i1, x_i2, x_i3,…, x_ip]. Поскольку y_i является случайной величиной, которая следует распределению Пуассона, вы можете увидеть разные значения y_i в каждом из 1000 наблюдений.

Предположим, что из 1000 наблюдаемых значений y_i вы наблюдаете 874 нулевых значения. Вы определяете, что из этих 874 нулевых значений обычное распределение Пуассона, которое вы приняли для y_i, сможет объяснить только до 7 нулевых значений. Таким образом, оставшиеся 867 нулевых значений являются избыточными нулевыми наблюдениями. Итак, для строки ith в вашем наборе данных ϕ_i = 867/1000 = 0,867.

Когда в наборе данных нет лишних нулей в зависимой переменной, значение ϕ оказывается равным нулю, и PMF модели ZIP уменьшается до PMF стандартная модель Пуассона (вы можете легко проверить это, установив для ϕ значение 0 в PMF модели ZIP).

Как оценить ϕ?

Итак, как мы можем оценить значение ϕ_i?

Простой и грубый способ оценить ϕ_i - установить для каждого ϕ _i следующее соотношение:

Возможно, более реалистичным способом вычисления ϕ_i является его оценка как функция регрессионных переменных X. Обычно это делается путем преобразования переменной y в двоичную случайную величину 0/1 y ' ( y_prime), который принимает значение 0, если базовый y равен 0, и 1 во всех остальных случаях. Затем мы применяем модель логистической регрессии к преобразованному y ’. Затем мы обучаем модель логистической регрессии на наборе данных [X, y '], и она дает вектор подобранных вероятностей µ_fitted. = [µ_1, µ_2, µ_3,…, µ_n], (потому что это то, что делает модель логистической регрессии) .

Как только мы получили вектор µ_fitted, мы просто устанавливаем его на вектор ϕ.

Таким образом, [ϕ_1 = µ_1, ϕ_2 = µ_2, ϕ_3 = µ_3,…, ϕ_n = µ_n].

Вышеупомянутый процесс оценки ϕ проиллюстрирован ниже :

После оценки вектора ϕ мы вставляем его в функции вероятности модели ZIP и используем то, что известно как M aximum Техника стимуляции L вероятности E (MLE) для обучения модели ZIP на наборе данных с избыточными счетчиками.

Please see my article on Poisson Regression Model for an explanation of how MLE works.

На следующем рисунке показана последовательность обучения модели ZIP:

К счастью, существует множество статистических пакетов, которые автоматизируют всю эту процедуру оценки ϕ и использования предполагаемого ϕ для обучения модель ZIP с использованием техники MLE для вашего набора данных.

В оставшейся части этой статьи мы будем использовать библиотеку Python statsmodels для создания и обучения модели ZIP в одной строке кода.

Как обучить модель ZIP с помощью Python

В нашем руководстве Python по модели ZIP мы будем использовать набор данных о походах, совершенных 250 группами людей:

Набор данных доступен здесь. Вот несколько важных особенностей этого набора данных:

  • Кемперы могли или не могли ловить рыбу во время поездки.
  • Если бы группа рыбачила, они бы не поймали ни одной рыбы или больше.
  • Мы хотим оценить не только, сколько рыбы было поймано (если рыбалка была сделана кемпинговой группой), но и вероятность того, что кемпинговая группа вообще поймала какую-либо рыбу.

Таким образом, задействованы два различных процесса генерации данных:

  1. Процесс, определяющий, успешно ли занималась рыбная ловля кемпинговая группа: модель ZIP будет внутренне использовать модель логистической регрессии, которая была объяснена ранее, для моделирования этого бинарного процесса.
  2. Второй процесс, который определяет, сколько рыбы было поймано кемпинговой группой, учитывая, что группа поймала хотя бы одну рыбу: модель ZIP будет использовать обычную модель Пуассона для моделирования этого второго процесса.

Переменные в наборе данных

Набор данных о походах содержит следующие переменные:

FISH_COUNT: количество пойманной рыбы. Это будет наша зависимая переменная y.

LIVE_BAIT: двоичная переменная, указывающая, использовалась ли живица.

КАМПЕР: использовала ли рыболовная группа автофургон.

ЛИЦА: общее количество человек в рыболовной группе. Обратите внимание, что в некоторых группах никто из них, возможно, не ловил рыбу.

ДЕТИ: количество детей в кемпинге.

Вот частотное распределение зависимой переменной FISH_COUNT:

Как мы видим, в этом наборе данных могут быть лишние нули. Мы обучим модель ZIP на этом наборе данных, чтобы проверить эту теорию и, надеюсь, добиться лучшего соответствия, чем обычная модель Пуассона.

Цель регрессии

Наши цели регрессии для этого набора данных следующие:

Предскажите количество пойманной рыбы (FISH_COUNT) кемпинговой группой на основе значений переменных LIVE_BAIT, CAMPER, PERSONS и CHILDREN.

Стратегия регрессии

Наша стратегия регрессии будет следующей:

  1. FISH_COUNT будет зависимой переменной y, а [LIVE_BAIT, CAMPER, PERSONS and CHILDREN] - независимыми переменными X .
  2. Мы будем использовать библиотеку statsmodels Python для обучения модели регрессии ZIP на данных (y, X). установленный.
  3. Мы сделаем некоторые прогнозы, используя модель ZIP на тестовом наборе данных, который модель не видела во время обучения.

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

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

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

df = pd.read_csv('fish.csv', header=0)

Давайте напечатаем несколько верхних строк набора данных:

print(df.head(10))

Давайте также распечатаем частотное распределение значений FISH_COUNT:

df.groupby('FISH_COUNT').count()

Создайте наборы данных для обучения и тестирования. Обратите внимание, что на данный момент мы не выполняем стратифицированное случайное разбиение:

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)))
>> Training data set length=196
>> Testing data set length=54

Задайте выражение регрессии в нотации Пэтси. Мы говорим Пэтси, что FISH_COUNT является нашей зависимой переменной y и зависит от переменных регрессии LIVE_BAIT, CAMPER, PERSONS и CHILDREN:

expr = 'FISH_COUNT ~ LIVE_BAIT  + CAMPER + CHILDREN + PERSONS'

Давайте воспользуемся Пэтси, чтобы выделить матрицы 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 ZeroInflatedPoisson, давайте построим и обучим модель регрессии ZIP на наборе обучающих данных.

Но прежде чем мы это сделаем, позвольте мне объяснить, как использовать два параметра, которые принимает конструктор класса:

  • инфляция: класс модели ZeroInflatedPoisson будет внутренне использовать модель LogisticRegression для оценки параметра ϕ. Поэтому мы устанавливаем параметр модели инфляция в 'logit'. Мы также можем поэкспериментировать с установкой для него других функций биномиальных ссылок, таких как «пробит».
  • exog_infl: Мы также хотим попросить модель ZIP оценить ϕ как функцию того же набора регрессионных переменных, что и родительская модель, а именно : LIVE_BAIT, CAMPER, PERSONS и CHILDREN. Следовательно, мы устанавливаем для параметра exog_infl значение X_train. Если вы хотите использовать только подмножество X_train, вы можете это сделать или можете установить для exog_infl совершенно другой набор регрессионных переменных.

В приведенной ниже строке строится и обучается модель ZIP на нашем наборе обучающих данных в одной строке кода.

zip_training_results = sm.ZeroInflatedPoisson(endog=y_train, exog=X_train, exog_infl=X_train, inflation='logit').fit()

Распечатайте сводку обучения:

print(zip_training_results.summary())

Вот краткое изложение обучения (я выделил важные элементы в выходных данных):

Интерпретация результатов тренировки

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

Обратите внимание, что модель логистической регрессии не нашла переменных Intercept, LIVE_BAIT и CAMPER, полезных для оценки ϕ. Их коэффициенты регрессии оказались НЕ статистически значимыми при уровне достоверности 95%, на что указывают соответствующие значения p:
inflate_Intercept = 0,425,
inflate_LIVE_BAIT = 0,680 и
inflate_CAMPER = 0,240,
, которые все больше 0,05 (т. е. порог ошибки 5% ).

Наблюдение 1

Единственными двумя переменными, которые модель логистической регрессии определила как полезные для оценки вероятности того, была ли поймана какая-либо рыба, были ДЕТИ и ЧЕЛОВЕКИ.

Наблюдение 2

Коэффициент регрессии PERSONS отрицательный (inflate_PERSONS -1.2193), что означает, что по мере увеличения количества людей в группе кемпинга вероятность того, что эта группа не поймает рыбу, уменьшается. Это соответствует нашей интуиции.

В красном поле содержится информация о переменных, которые родительская модель Пуассона использовала для оценки FISH_COUNT при условии, что FISH_COUNT> 0.

Наблюдение 3

Мы можем видеть, что коэффициенты для всех 5 регрессионных переменных статистически значимы на уровне достоверности 99%, о чем свидетельствует их значение p, которое меньше 0,01. Фактически, значение p меньше 0,001 для всех 5 переменных, следовательно, оно отображается как 0,000.

Наблюдение 4

Коэффициент для ДЕТЕЙ отрицательный (ДЕТИ -1,0810), что означает, что по мере увеличения количества детей в группе кемпинга количество рыбы, пойманной этой группой, уменьшается!

Наблюдение 5

Максимальное логарифмическое правдоподобие этой модели составляет -566,43. Это значение полезно для сравнения степени согласия модели с другими моделями (см. Забавное упражнение в конце статьи).

Наблюдение 6

Наконец, обратите внимание, что алгоритм обучения модели ZIP не смог сойтись на наборе обучающих данных, о чем свидетельствует следующее:

Если бы они сошлись, возможно, получилось бы лучше.

Прогноз

Мы получим прогнозы модели ZIP на тестовом наборе данных и вычислим среднеквадратическую ошибку w.r.t. фактические значения:

zip_predictions = zip_training_results.predict(X_test,exog_infl=X_test)
predicted_counts=np.round(zip_predictions)
actual_counts = y_test[dep_var]

print('ZIP RMSE='+str(np.sqrt(np.sum(np.power(np.subtract(predicted_counts,actual_counts),2)))))
>> ZIP RMSE=55.65069631190611

Давайте изобразим прогнозируемое количество рыбы в сравнении с фактическим:

fig = plt.figure()
fig.suptitle('Predicted versus actual counts using the ZIP model')
predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted')
actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual')
plt.legend(handles=[predicted, actual])
plt.show()

Мы видим следующий сюжет:

На этом мы завершаем рассмотрение модели регрессии Пуассона с нулевым раздутием.

Веселое упражнение

  1. Прочтите следующую статью: Иллюстрированное руководство по модели регрессии Пуассона
  2. Используя Python и statsmodels, обучите стандартную модель Пуассона на наборе данных о походах и сравните производительность модели с производительностью модели ZIP. Вы можете быть удивлены тем, что найдете.
  3. Сравните степень согласия двух моделей на наборе обучающих данных, используя значение максимального логарифмического правдоподобия, сообщаемое statsmodels в сводке по обучению обеих моделей. Чем больше Максимальный LL, тем лучше степень соответствия.
  4. Сравните оценки RMSE двух моделей на наборе тестовых данных.

Удачного моделирования!

Предлагаемые темы для дальнейшего чтения

Спасибо за внимание! Я пишу на темы науки о данных, уделяя особое внимание регрессии и анализу временных рядов.

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