¿Cuándo las regresiones binomiales negativas y de Poisson se ajustan a los mismos coeficientes?

He notado que en R, Poisson y las regresiones binomiales negativas (NB) siempre parecen ajustarse a los mismos coeficientes para predictores categóricos, pero no continuos.

Por ejemplo, aquí hay una regresión con un predictor categórico:

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í hay un ejemplo con un predictor continuo, donde Poisson y NB se ajustan a diferentes 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 datos de conteo, y los modelos no son significativos ...)

Luego recodifico el predictor en un factor, y los dos modelos nuevamente se ajustan a 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, la regresión binomial negativa de Joseph Hilbe da un ejemplo (6.3, pg 118-119) donde un predictor categórico, el sexo, se ajusta con coeficientes ligeramente diferentes por el Poisson ( $b=0.883$ ) y NB ( $b=0.881$ ). Él dice: “Las tasas de incidencia entre los modelos de Poisson y NB son bastante similares. Esto no es sorprendente dada la proximidad de [correspondiente a en R] a cero ".$\alpha$$1/\theta$

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

Entonces, ¿por qué los coeficientes son exactamente iguales? ¿Y por qué solo para predictores categóricos?

---

Editar # 1

Aquí hay un ejemplo con dos predictores no ortogonales. De hecho, los coeficientes ya no son los mismos:

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, al incluir otro predictor, los modelos se ajustan a diferentes coeficientes incluso cuando el nuevo predictor es continuo. Entonces, ¿tiene algo que ver con la ortogonalidad de las variables ficticias que creé 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))

Ha descubierto una propiedad íntima, pero genérica, de GLM ajustada por la máxima probabilidad . El resultado desaparece una vez que se considera el caso más simple de todos: ¡Ajustar un solo parámetro a una sola observación!

Una respuesta frase : Si todo lo que importa es apropiado medios separados de subconjuntos disjuntos de nuestra muestra, a continuación, GLM siempre dará μ j = ˉ y j para cada subconjunto$\hat\mu_j = \bar y_j$ , por lo que la estructura de error real y parametrización de la densidad tanto vuelto irrelevante a la (punto) estimación!$j$

Un poco más : ajustar factores categóricos ortogonales por máxima probabilidad es equivalente a ajustar medios separados para separar subconjuntos de nuestra muestra, por lo que esto explica por qué Poisson y los GLM binomiales negativos producen las mismas estimaciones de parámetros. De hecho, lo mismo sucede si utilizamos la regresión de Poisson, negbin, gaussiana, gaussiana inversa o gamma (ver más abajo). En el caso de Poisson y Negbin, la función de enlace predeterminada es el enlace de , pero es una pista falsa; Si bien esto produce las mismas estimaciones de parámetros sin procesar , veremos a continuación que esta propiedad realmente no tiene nada que ver con la función de enlace.$\log$

Cuando estamos interesados ​​en una parametrización con más estructura, o que depende de predictores continuos, entonces la estructura de error asumida se vuelve relevante debido a la relación media-varianza de la distribución, ya que se relaciona con los parámetros y la función no lineal utilizada para modelar el condicional medio.

GLMs y familias de dispersión exponencial: curso intensivo

Una familia de dispersión exponencial en forma natural es tal que la densidad logarítmica es de la forma

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

Aquí es el parámetro natural y ν es el parámetro de dispersión . Si se conociera ν , esto sería solo una familia exponencial estándar de un parámetro. Todos los GLM considerados a continuación suponen un modelo de error de esta familia.$\theta$$\nu$$\nu$

Considere una muestra de una sola observación de esta familia. Si nos ajustamos a por máxima verosimilitud, obtenemos que y = b ' ( θ ) , independientemente del valor de ν . Esto se extiende fácilmente en el caso de una muestra iid ya que las probabilidades de registro añaden, produciendo ˉ y = b ' ( θ$\theta$$y = b'(\hat\theta)$$\nu$$\bar y = b'(\hat\theta)$ .

Pero, también sabemos, debido a la agradable regularidad de la densidad logarítmica en función de , que $\theta$ Entonces, de hecho 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$

Dado que las estimaciones de máxima verosimilitud son invariantes bajo transformaciones, esto significa que para esta familia de densidades.$\bar y = \hat\mu$

Ahora, en un GLM, modelamos como μ i = g - 1 ( x T i β ) donde g es la función de enlace. Pero si x i es un vector de todos los ceros a excepción de un solo 1 en la posición j , entonces μ i = g ( β j ) . La probabilidad de la GLM luego se factoriza de acuerdo con los β j 's y procedemos como se indica arriba. Este es precisamente el caso de los factores ortogonales.$\mu_i$$\mu_i = g^{-1}(\mathbf x_i^T \beta)$$g$$\mathbf x_i$$j$$\mu_i = g(\beta_j)$$\beta_j$

¿Qué tienen de diferentes los predictores continuos?

Cuando los predictores son continuos o categóricos, pero no pueden reducirse a una forma ortogonal, entonces la probabilidad ya no se factoriza en términos individuales con una media separada que depende de un parámetro separado. En este punto, la estructura y la función de enlace de error no entran en juego.

Si uno atraviesa el álgebra (tedioso), las ecuaciones de probabilidad se convierten en para todo j = 1 , ... , p donde λ i = x T i β . Aquí, losparámetros β y ν entran implícitamente a través de la relación de enlace μ i = g ( λ i ) = g ( x T i β ) y la varianza σ 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$ .

De esta manera, la función de enlace y el modelo de error asumido se vuelven relevantes para la estimación.

Ejemplo: el modelo de error (casi) no importa

En el siguiente ejemplo, generamos datos aleatorios binomiales negativos 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$ se utiliza ).

Luego ajustamos estos datos utilizando cinco GLM diferentes, cada uno con un enlace de : ( a ) binomio negativo, ( b ) Poisson, ( c ) gaussiano, ( d ) gaussiano inverso y ($\log$ e ) GLM gamma. Todos estos son ejemplos de familias de dispersión exponencial.

De la tabla, podemos ver que las estimaciones de los parámetros son idénticas , aunque algunos de estos GLM son para datos discretos y otros son para datos continuos, y algunos son para datos no negativos mientras que otros no.

      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 encabezado proviene del hecho de que el procedimiento de ajuste fallará si las observaciones no caen dentro del dominio de la densidad particular. Por ejemplo, si tuviéramos recuentos generados aleatoriamente en los datos anteriores, entonces el Gamma GLM no podría converger ya que los Gamma GLM requieren datos estrictamente positivos.$0$

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

$\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

La advertencia en el encabezado simplemente se refiere al hecho de que las estimaciones sin procesar variarán con la función de enlace, pero las estimaciones de parámetros medios implícitos no lo harán.

Código 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))

Para ver lo que está sucediendo aquí, es útil primero hacer la regresión sin la intercepción, ya que una intercepción en una regresión categórica con un solo predictor no tiene sentido:

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

Dado que Poisson y las regresiones binomiales negativas especifican el registro del parámetro medio, entonces, para la regresión categórica, exponiendo los coeficientes le dará el parámetro medio real 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 las medias reales sobre los diferentes valores de categoría:

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

$\lambda$

$\lambda$ que maximice la probabilidad de las observaciones en cada categoría particular. El estimador de máxima verosimilitud para la distribución de Poisson es simplemente la media muestral, razón por la cual los coeficientes de regresión son precisamente los (registros de) las medias muestrales para cada categoría.

$\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}$$ so the maximum is attained when $\lambda=\bar{x}$.

The reason that you don't get the same coefficients for continuous data is because in a continuous regression, $\text{log}(\lambda)$ is no longer going to be a piecewise constant function of the predictor variables, but a linear one. Maximizing the likelihood function in this case will not reduce to independently fitting a value $\lambda$ for disjoint subsets of the data, but will rather be a nontrivial problem that is solved numerically, and is likely to produce different results for different likelihood functions.

Similarly, if you have multiple categorical predictors, despite the fact that the fitted model will ultimately specify $\lambda$ as a piecewise constant function, in general there will not be enough degrees of freedom to allow $\lambda$ to be determined independently for each constant segment. For example, suppose that you have $2$ predictors with $5$ categories each. In this case, your model has $10$ degrees of freedom, whereas there are $5*5=25$ unique different combinations of the categories, each of which will have its own fitted value of $\lambda$. So, assuming that the intersections of these categories are non-empty (or at least that $11$ of them are nonempty), the likelihood maximization problem again becomes nontrivial and will generally produce different outcomes for Poisson versus negative binomial or any other distribution.