Когда пуассоновская и отрицательная биномиальные регрессии соответствуют одинаковым коэффициентам?

Я заметил, что в R регрессии Пуассона и отрицательная биномиальная (NB) всегда соответствуют одинаковым коэффициентам для категориальных, но не непрерывных предикторов.

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

data(warpbreaks)
library(MASS)

rs1 = glm(breaks ~ tension, data=warpbreaks, family="poisson")
rs2 = glm.nb(breaks ~ tension, data=warpbreaks)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Вот пример с непрерывным предиктором, в котором Пуассон и NB соответствуют разным коэффициентам:

data(cars)
rs1 = glm(dist ~ speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

(Конечно, это не данные подсчета, и модели не имеют смысла ...)

Затем я перекодирую предиктор в фактор, и две модели снова соответствуют одинаковым коэффициентам:

library(Hmisc)
speedCat = cut2(cars$speed, g=5) 
#you can change g to get a different number of bins

rs1 = glm(cars$dist ~ speedCat, family="poisson")
rs2 = glm.nb(cars$dist ~ speedCat)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

Однако в « Отрицательной биномиальной регрессии» Джозефа Хильбе приведен пример (6.3, стр. 118–119), в котором категорический предиктор «пол» соответствует слегка отличающимся коэффициентам по Пуассону ( ) и NB ( ). Он говорит: «Коэффициенты заболеваемости между моделями Пуассона и NB довольно похожи. Это неудивительно, учитывая близость [соответствующего в R] к нулю ».b = 0,881 α 1 /$b=0.883$$b=0.881$$\alpha$$1/\theta$

Я понимаю это, но в приведенных выше примерах summary(rs2)говорится, что оценивается в 9,16 и 7,93 соответственно.$\theta$

Так почему же коэффициенты одинаковы? И почему только для категориальных предикторов?

---

Редактировать # 1

Вот пример с двумя неортогональными предикторами. Действительно, коэффициенты уже не одинаковы:

data(cars)

#make random categorical predictor
set.seed(1); randomCats1 = sample( c("A","B","C"), length(cars$dist), replace=T)
set.seed(2); randomCats2 = sample( c("C","D","E"), length(cars$dist), replace=T)

rs1 = glm(dist ~ randomCats1 + randomCats2, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + randomCats2, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

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

rs1 = glm(dist ~ randomCats1 + speed, data=cars, family="poisson")
rs2 = glm.nb(dist ~ randomCats1 + speed, data=cars)

#compare coefficients
cbind("Poisson"=coef(rs1), "NB"=coef(rs2))

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

Ответ из одного предложения : если все, о чем мы заботимся, это подгоняет отдельные средства для непересекающихся подмножеств нашей выборки, то GLM всегда будут давать для каждого подмножества , поэтому фактическая структура ошибок и параметризация плотности становятся не имеет отношения к (точечной) оценке!$\hat\mu_j = \bar y_j$$j$

Немного больше : подгонка ортогональных категориальных факторов по максимальной вероятности эквивалентна подгонке отдельных средств к непересекающимся подмножествам нашей выборки, поэтому это объясняет, почему пуассоновские и отрицательные биномиальные GLM дают одинаковые оценки параметров. Действительно, то же самое происходит независимо от того, используем ли мы пуассоновскую, негбиновую, гауссовскую, обратную гауссовскую или гамма-регрессию (см. Ниже). В случае Пуассона и Негбина функция ссылки по умолчанию - ссылка , но это красная сельдь; в то время как это дает те же самые необработанные оценки параметров, мы увидим ниже, что это свойство на самом деле не имеет ничего общего с функцией связи.$\log$

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

GLM и экспоненциальные дисперсионные семейства: ускоренный курс

Экспоненциальной дисперсия семьи в естественной форме является одним из таких , что плотность лог имеет вид

$$\log f(y;\,\theta,\nu) = \frac{\theta y - b(\theta)}{\nu} + a(y,\nu) \>.$$

Здесь - естественный параметр, а - параметр дисперсии . Если бы были известны, это было бы просто одноэлементное экспоненциальное семейство. Все GLM, рассмотренные ниже, предполагают модель ошибок из этого семейства.ν $\theta$$\nu$$\nu$

Рассмотрим пример одного наблюдения из этой семьи. Если мы подберем по максимальной вероятности, мы получим, что , независимо от значения . Это легко распространяется на случай iid-примера, поскольку логарифмические правдоподобия складываются, приводя к .у = Ь ' ( θ ) N , ˉ у = Ь ' ( θ$\theta$$y = b'(\hat\theta)$$\nu$$\bar y = b'(\hat\theta)$

Но мы также знаем, что из-за хорошей регулярности плотности бревна как функции от это Итак, на самом деле .$\theta$b ' ( θ ) = E Y =

$$\frac{\partial}{\partial \theta} \mathbb E \log f(Y;\theta,\nu) = \mathbb E \frac{\partial}{\partial \theta} \log f(Y;\theta,\nu) = 0 \>.$$ $b'(\theta) = \mathbb E Y = \mu$

Поскольку оценки максимального правдоподобия инвариантны относительно преобразований, это означает, что для этого семейства плотностей.$\bar y = \hat\mu$

Теперь в GLM мы моделируем как где - функция связи. Но если является вектором всех нулей, кроме одного 1 в позиции , то . Затем вероятность GLM разлагается в соответствии с , и мы действуем, как указано выше. Это как раз случай ортогональных факторов.μ i = g - 1 ( x T i β ) g x i j μ i = g ( β j ) β $\mu_i$$\mu_i = g^{-1}(\mathbf x_i^T \beta)$$g$$\mathbf x_i$$j$$\mu_i = g(\beta_j)$$\beta_j$

Чем же отличаются непрерывные предикторы?

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

Если пройти через (утомительную) алгебру, уравнения правдоподобия становятся для всех где . Здесь параметры и вводятся неявно через отношение ссылки и дисперсию .j = 1 , … , p λ i = x T i β β ν μ i = g ( λ i ) = g ( x T i β ) σ 2

$$\sum_{i=1}^n \frac{(y_i - \mu_i)x_{ij}}{\sigma_i^2}\frac{\partial \mu_i}{\partial \lambda_i} = 0\>,$$ $j = 1,\ldots,p$$\lambda_i = \mathbf x_i^T \beta$$\beta$$\nu$$\mu_i = g(\lambda_i) = g(\mathbf x_i^T \beta)$$\sigma_i^2$

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

Пример: модель ошибки (почти) не имеет значения

В приведенном ниже примере мы генерируем отрицательные биномиальные случайные данные в зависимости от трех категориальных факторов. Каждое наблюдение происходит из одной категории, и используется один и тот же параметр дисперсии ( ).$k = 6$

Затем мы подгоняем к этим данным пять различных GLM, каждый со ссылкой : ( a ) отрицательный бином, ( b ) пуассоновский, ( c ) гауссовский, ( d ) обратный гауссовский и ( e ) гамма-GLM. Все это примеры экспоненциальных дисперсионных семейств.$\log$

Из таблицы видно, что оценки параметров идентичны , даже если некоторые из этих GLM предназначены для дискретных данных, а другие - для непрерывных, а некоторые - для неотрицательных данных, а другие - нет.

      negbin  poisson gaussian invgauss    gamma
XX1 4.234107 4.234107 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033 4.841033 4.841033

Предостережение в заголовке исходит из того факта, что процедура подбора потерпит неудачу, если наблюдения не попадают в область конкретной плотности. Например, если бы у нас было отсчетов, случайно сгенерированных в данных выше, то Gamma GLM не сможет сойтись, так как Gamma GLM требует строго положительных данных.$0$

Пример: функция ссылки (почти) не имеет значения

Используя те же данные, мы повторяем процедуру подбора данных с помощью Пуассона GLM с тремя различными функциями связи: ( a ) link, ( b ) идентификационная ссылка и ( c ) квадратно-корневая ссылка. В таблице ниже приведены оценки коэффициентов после преобразования обратно в логарифмическую параметризацию. (Таким образом, во втором столбце показано а в третьем столбце показано с использованием необработанного из каждого совпадения). Опять же, оценки идентичны.журнал ( β ) журнал ( β 2 ) $\log$$\log(\hat \beta)$$\log(\hat \beta^2)$$\hat\beta$

> coefs.po
         log       id     sqrt
XX1 4.234107 4.234107 4.234107
XX2 4.790820 4.790820 4.790820
XX3 4.841033 4.841033 4.841033

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

Код R

# Warning! This code is a bit simplified for compactness.
library(MASS)
n <- 5
m <- 3
set.seed(17)
b <- exp(5+rnorm(m))
k <- 6

# Random negbin data; orthogonal factors
y <- rnbinom(m*n, size=k, mu=rep(b,each=n))
X <- factor(paste("X",rep(1:m,each=n),sep=""))

# Fit a bunch of GLMs with a log link
con <- glm.control(maxit=100)
mnb <- glm(y~X+0, family=negative.binomial(theta=2))
mpo <- glm(y~X+0, family="poisson")
mga <- glm(y~X+0, family=gaussian(link=log), start=rep(1,m), control=con)
miv <- glm(y~X+0, family=inverse.gaussian(link=log), start=rep(2,m), control=con)
mgm <- glm(y~X+0, family=Gamma(link=log), start=rep(1,m), control=con)    
coefs <- cbind(negbin=mnb$coef, poisson=mpo$coef, gaussian=mga$coef
                   invgauss=miv$coef, gamma=mgm$coef)

# Fit a bunch of Poisson GLMs with different links.
mpo.log  <- glm(y~X+0, family=poisson(link="log"))
mpo.id   <- glm(y~X+0, family=poisson(link="identity"))
mpo.sqrt <- glm(y~X+0, family=poisson(link="sqrt"))   
coefs.po <- cbind(log=mpo$coef, id=log(mpo.id$coef), sqrt=log(mpo.sqrt$coef^2))

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

> rs1 = glm(breaks ~ tension-1, data=warpbreaks, family="poisson")
> rs2 = glm.nb(breaks ~ tension-1, data=warpbreaks)

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

>  exp(cbind("Poisson"=coef(rs1), "NB"=coef(rs2)))
          Poisson       NB
tensionL 36.38889 36.38889
tensionM 26.38889 26.38889
tensionH 21.66667 21.66667

Эти параметры соответствуют фактическим средним значениям для разных значений категории:

> with(warpbreaks,tapply(breaks,tension,mean))
       L        M        H 
36.38889 26.38889 21.66667 

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

Для распределения Пуассона понятно, почему это происходит. Подходит только один параметр, и, следовательно, максимизация общей вероятности модели с одним категориальным предиктором эквивалентна независимому поиску который максимизирует вероятность для наблюдений в каждой конкретной категории. Оценка максимального правдоподобия для распределения Пуассона - это просто среднее значение выборки, поэтому коэффициенты регрессии являются именно (логарифмами) выборочного среднего значения для каждой категории.$\lambda$

Для отрицательного бинома это не так просто, потому что есть два параметра: и параметр формы . Более того, регрессия соответствует одной , охватывающей весь набор данных, поэтому в этой ситуации категориальная регрессия не просто эквивалентна подгонке совершенно отдельной модели для каждой категории. Однако, изучая функцию правдоподобия, мы можем увидеть, что для любой данной тэты функция правдоподобия снова максимизируется путем установки в качестве значения образца:θ θ $\lambda$$\theta$$\theta$$\lambda$

λ=ˉ

$$\begin{align} L(X,\lambda,\theta) &= \prod \left(\frac{\theta}{\lambda+\theta}\right)^\theta \frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\left(\frac{\lambda}{\lambda+\theta}\right)^{x_i}\\ \log L(X,\lambda,\theta) &= \sum\theta\left(\text{log}\theta-\text{log}(\lambda+\theta)\right) +x_i\left(\text{log}\lambda-\text{log}(\lambda+\theta)\right) +\log\left(\frac{\Gamma(\theta + x_i)}{x_i!\Gamma(\theta)}\right)\\ \frac{d}{d\lambda}\log L(X,\lambda,\theta) &= \sum \frac{x_i}{\lambda}-\frac{\theta+x_i}{\lambda+\theta} =n\left(\frac{\bar{x}}{\lambda}-\frac{\bar{x}+\theta}{\lambda+\theta}\right), \end{align}$$ поэтому максимум достигается, когда .$\lambda=\bar{x}$

Причина, по которой вы не получаете одинаковые коэффициенты для непрерывных данных, заключается в том, что в непрерывной регрессии больше не будет кусочно-постоянной функцией переменных предиктора, а будет линейной. Максимизация функции правдоподобия в этом случае не приведет к независимой подгонке значения для непересекающихся подмножеств данных, но скорее будет нетривиальной задачей, которая решается численно и, вероятно, даст разные результаты для разных функций правдоподобия.$\text{log}(\lambda)$$\lambda$

Точно так же, если у вас есть несколько категориальных предикторов, несмотря на тот факт, что подобранная модель в конечном итоге будет указывать как кусочно-постоянную функцию, в общем случае не будет достаточно степеней свободы, чтобы позволить быть определенным независимо для каждого постоянного сегмента. Например, предположим, что у вас есть предиктора с категориями в каждом. В этом случае ваша модель имеет степеней свободы, в то время как существует уникальных различных комбинаций категорий, каждая из которых будет иметь свое собственное подогнанное значение . Таким образом, предполагая, что пересечения этих категорий непусты (или, по крайней мере, чтоλ 2 5 10 5 ∗ 5 = 25 λ $\lambda$$\lambda$$2$$5$$10$$5*5=25$$\lambda$$11$ из них непустые), проблема максимизации правдоподобия снова становится нетривиальной и, как правило, приводит к различным результатам для Пуассона по сравнению с отрицательным биномиальным или любым другим распределением.