¿Cuándo Poisson y binomial negativa regresiones ajuste de los mismos coeficientes?

Question

Me he dado cuenta de que en R, Poisson y binomial negativa (NB) regresiones siempre parecen ajustarse a los mismos coeficientes para categórica, pero no continuo, los predictores.

Por ejemplo, he aquí una regresión con un predictores categóricos:

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))

Aquí está un ejemplo con un predictor continuo, donde la distribución de Poisson y NB adaptarse a los distintos coeficientes:

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))

(Por supuesto estos no son los datos de recuento, y los modelos no son significativas...)

Entonces me recodificar el predictor en un factor, y los dos modelos más aptos los mismos coeficientes:

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))

Sin embargo, José Hilbe de la Regresión Binomial Negativa da un ejemplo (6.3, pg 118-119) donde un predictores categóricos, sexo, está en condiciones ligeramente diferentes coeficientes de Poisson ($b=0.883$) y NB ($b=0.881$). Él dice: "La tasa de incidencia de las relaciones entre la distribución de Poisson y NB modelos son bastante similares. Esto no es sorprendente dada la proximidad de $\alpha$ [correspondiente a $1/\theta$ R] a cero".

Entiendo que esto, pero en los ejemplos anteriores, summary(rs2) nos dice que $\theta$ se estima en 9.16 y 7.93 respectivamente.

Entonces, ¿por qué son los coeficientes exactamente el mismo? Y por qué sólo para los predictores categóricos?

---

Edición #1

Aquí está un ejemplo con dos no ortogonal predictores. De hecho, los coeficientes no son el mismo:

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))

Y, como otro predictor causas de los modelos para adaptarse a los distintos coeficientes incluso cuando el nuevo predictor es continua. Así que, es algo para hacer con la ortogonalidad de las variables ficticias que he creado en mi ejemplo original?

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))

Answer 1

Usted ha descubierto un íntimo, pero genérico, propiedad de GLMs ajuste por máxima verosimilitud. El resultado se retira una vez que se considera el caso más simple de todos: el Ajuste de un solo parámetro para una sola observación!

Una frase de respuesta: Si todos nos preocupamos por que es conveniente separar los medios para distintos subgrupos de la muestra, a continuación, GLMs siempre será el rendimiento de $\hat\mu_j = \bar y_j$ para cada subconjunto $j$, por lo que el error real de la estructura y configuración de parámetros de la densidad de ambos se vuelven irrelevantes para el (punto) estimación!

Un poco más: Ajuste ortogonal categórica factores de máxima verosimilitud es equivalente a la colocación de distintas vías para distintos subgrupos de la muestra, esto explica el por qué de Poisson y binomial negativa GLMs el rendimiento de la misma estimaciones de los parámetros. De hecho, lo mismo sucede si utilizamos Poisson, negbin, Gaussiana inversa Gaussiana o Gamma de regresión (ver más abajo). En la distribución de Poisson y negbin caso, el valor predeterminado función de enlace es el $\log$ enlace, sino que es un arenque rojo; mientras que los rendimientos de la misma primas estimaciones de los parámetros, veremos más adelante que esta propiedad no tiene realmente nada que ver con la función de enlace.

Cuando estamos interesados en un parametrización con más estructura, o que depende de la continua predictores, entonces el error que se asume la estructura se convierte en relevante debido a la media de la varianza de la relación de la distribución, ya que se refiere a los parámetros y la función no lineal utilizado para el modelado de los medios condicionales.

GLMs y exponencial de la dispersión de las familias: curso intensivo

Un exponencial de la dispersión de la familia en forma natural es uno de esos en que el registro de densidad es de la forma $$ \log f(y;\,\theta,\nu) = \frac{\theta y - b(\theta)}{\nu} + a(y,\nu) \>. $$

Aquí $\theta$ es el natural parámetro y $\nu$ es el parámetro de dispersión. Si $\nu$ eran conocidos, se trata de una norma de un parámetro exponencial de la familia. Todos los GLMs considera por debajo de asumir un modelo de error de esta familia.

Se toma una muestra de una sola observación de esta familia. Si encajamos $\theta$ por máxima verosimilitud, tenemos que $ $ y = b'(\hat\theta)$, independientemente del valor de $\nu$. Esto fácilmente se extiende para el caso de un alcoholímetro de la muestra ya que el registro de las probabilidades agregar, produciendo $\bar y = b'(\hat\theta)$.

Pero, también sabemos que, debido a la agradable regularidad del registro de densidad como función de la $\theta$, que $$ \frac{\partial}{\partial \theta} \mathbb E \log f(Y;\theta,\nu) = \mathbb E \frac{\partial}{\partial \theta} \log f(Y;\theta,\nu) = 0 \>. $$ Así que, de hecho $b'(\theta) = \mathbb E Y = \mu$.

Desde el máximo de estimaciones de probabilidad son invariantes bajo transformaciones, esto significa que $ \bar y = \hat\mu $ para esta familia de densidades.

Ahora, en un GLM, el modelo de $\mu_i$ como $\mu_i = g^{-1}(\mathbf x_i^T \beta)$ donde $g$ es la función de enlace. Pero si $\mathbf x_i$ es un vector de ceros excepto un solo 1 en la posición $j$, entonces $\mu_i = g(\beta_j)$. La probabilidad de que el GLM, a continuación, factorizes de acuerdo a la $\beta_j$'s y procedemos como en el anterior. Este es precisamente el caso de factores ortogonales.

¿Qué tan diferente acerca de la continua predictores?

Cuando los predictores son continuos o son categóricos, pero no puede ser reducida a una forma ortogonal, entonces la probabilidad de no más factores en términos individuales con zona media dependiendo de parámetros por separado. En este punto, el error de la estructura y función de enlace hacer entrar en juego.

Si uno de los pedales a través de la (tedioso) el álgebra, la probabilidad de ecuaciones son $$ \sum_{i=1}^n \frac{(y_i - \mu_i)x_{ij}}{\sigma_i^2}\frac{\partial \mu_i}{\parcial \lambda_i} = 0\>, $$ para todos $j = 1,\ldots,p$ donde $\lambda_i = \mathbf x_i^T \beta$. Aquí, el $\beta$ y $\nu$ accede a los parámetros de forma implícita a través del vínculo de la relación $\mu_i = g(\lambda_i) = g(\mathbf x_i^T \beta)$ y variación de $\sigma_i^2$.

De esta manera, la función de enlace y asumió el modelo de error de ser pertinente para la estimación.

Ejemplo: El modelo de error de (casi) no importa

En el ejemplo siguiente, se genera una binomial negativa datos aleatorios dependiendo de tres factores categóricos. Cada observación proviene de una sola categoría y el mismo parámetro de dispersión ($k = 6$).

Nosotros, a continuación, ajuste a estos datos usando cinco GLMs, cada uno con un $\log$ enlace: (un) binomial negativa, (b) Poisson, (c) Gaussiano, (d) Inversa Gaussiana y (e) Gamma GLMs. Todos estos son ejemplos de la exponencial de la dispersión de las familias.

A partir de la tabla, podemos ver que las estimaciones de los parámetros son idénticos, aunque algunas de estas GLMs son para datos discretos y otros son de continua,y algunas son para no negativo de datos, mientras que otros no lo son.

      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

La advertencia en el título proviene del hecho de que el procedimiento de ajuste se producirá si las observaciones no caen dentro del dominio de la particular densidad. Por ejemplo, si tuviéramos $0$ recuentos generado aleatoriamente en los datos de arriba, a continuación, la Gamma GLM no convergen desde Gamma GLMs requieren estrictamente positivo de datos.

Ejemplo: La función de enlace (casi) no importa

Utilizando los mismos datos, volvemos a repetir el procedimiento de ajuste de los datos con una distribución de Poisson GLM con tres diferentes funciones de enlace: (una) $\log$ link; (b) la identidad de enlace y (c) raíz cuadrada de enlace. La siguiente tabla muestra los coeficientes estimados después de la conversión de la espalda para el registro de la parametrización. (Así, en la segunda columna showns $\log(\hat \beta)$ y la tercera muestra $\log(\hat \beta^2)$ utilizando el raw $\hat\beta$ de cada uno de los ataques). De nuevo, las estimaciones son idénticos.

> 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

La advertencia de que en el título se refiere simplemente al hecho de que las primas de las estimaciones varían con la función de enlace, pero que implicaba la media de las estimaciones de los parámetros no.

R código de

# 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))

Answer 2

Con el fin de ver lo que está pasando aquí, es útil primero para hacer la regresión sin la intersección, ya que una intercepción en una regresión categórica con un solo predictor es de sentido:

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

Desde Poisson y binomial negativa regresiones especificar el registro de la media del parámetro, luego de regresión categórica, exponentiating los coeficientes le dará la real significa parámetro para cada categoría:

>  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

Estos parámetros corresponden a los medios sobre los diferentes valores de la categoría:

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

Así que lo que está sucediendo es que la media del parámetro $\lambda$ en cada caso que maximiza la probabilidad es igual a la media de la muestra para cada categoría.

Para la distribución de Poisson es claro por qué ocurre esto. Sólo hay un parámetro para el ajuste, y maximizar así el global de la probabilidad de que un modelo con un único predictores categóricos es equivalente, de forma independiente, encontrar un $\lambda$, lo que maximiza la probabilidad de las observaciones en cada categoría en particular. El estimador de máxima verosimilitud para la distribución de Poisson es simplemente la media de la muestra, que es la razón de los coeficientes de regresión son precisamente los (registros) de los de la muestra significa para cada categoría.

Para binomial negativa no es tan simple, porque hay dos parámetros de ajuste: $\lambda$ y la forma de parámetros de $\theta$. Por otra parte, la regresión se ajusta a una sola de $\theta$ que cubre todo el conjunto de datos, por lo que en esta situación una regresión categórica no es simplemente equivalente a la colocación de un modelo para cada categoría. Sin embargo, mediante el examen de la función de probabilidad, podemos ver que para cualquier theta, la probabilidad función es la de volver a ser maximizada por la configuración de $\lambda$ a la media de la muestra:

\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) y= \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) y= \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} por lo que el máximo se alcanza cuando $\lambda=\bar{x}$.

La razón por la que usted no consigue los mismos coeficientes para datos continuos se debe a que en una continua regresión, $\text{log}(\lambda)$ ya no va a ser una constante a trozos función de las variables predictoras, pero uno lineal. Maximizar la probabilidad de la función en este caso no va a reducir, de forma independiente, el ajuste de un valor de $\lambda$ para distintos subconjuntos de los datos, sino que más bien va a ser un problema no trivial que se resuelve numéricamente, y es probable que produzca resultados diferentes para las diferentes funciones de probabilidad.

Del mismo modo, si usted tiene múltiples los predictores categóricos, a pesar del hecho de que el modelo ajustado en última instancia especificar $\lambda$ como una función constante a trozos, en general, no habrá suficientes grados de libertad para permitir que $\lambda$ a se determina de forma independiente para cada constante segmento. Por ejemplo, supongamos que usted tiene $2$ predictores con $5$ cada una de las categorías. En este caso, el modelo tiene $10$ grados de libertad, considerando que no hay $5*5=25$ únicos diferentes combinaciones de las categorías, cada una de las cuales tiene sus propias equipado valor de $\lambda$. Así que, asumiendo que las intersecciones de estas categorías son no vacíos (o al menos que $11$ de ellos son no vacías), la probabilidad de problema de maximización de nuevo se vuelve trivial y por lo general se producen diferentes resultados de Poisson frente binomial negativa o de cualquier otra distribución.