Quando Poisson e regressões binomiais negativas se ajustam aos mesmos coeficientes?

38

Notei que em R, Poisson e regressões binomiais negativas (RN) sempre parecem se encaixar nos mesmos coeficientes para preditores categóricos, mas não contínuos.

Por exemplo, aqui está uma regressão com um preditor 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))

insira a descrição da imagem aqui

Aqui está um exemplo com um preditor contínuo, em que Poisson e NB se ajustam 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))

insira a descrição da imagem aqui

(É claro que esses dados não são contados e os modelos não são significativos ...)

Então, recodifico o preditor em um fator e os dois modelos novamente se encaixam nos mesmos 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))

insira a descrição da imagem aqui

No entanto, a regressão binomial negativa de Joseph Hilbe dá um exemplo (6.3, pág. 118-119) em que um preditor categórico, sexo, se encaixa com coeficientes ligeiramente diferentes pelos Poisson ( ) e NB ( ). Ele diz: “As taxas de incidência entre os modelos de Poisson e RN são bastante semelhantes. Isso não é surpreendente, dada a proximidade de [correspondente a em R] a zero. ”b = 0,881 α 1 / θb=0,883b=0,881α1/θ

Entendo isso, mas nos exemplos acima, summary(rs2)diz-nos que é estimado em 9,16 e 7,93, respectivamente.θ

Então, por que os coeficientes são exatamente os mesmos? E por que apenas para preditores categóricos?


Editar # 1

Aqui está um exemplo com dois preditores não ortogonais. De fato, os coeficientes não são mais os mesmos:

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

insira a descrição da imagem aqui

Além disso, a inclusão de outro preditor faz com que os modelos se ajustem a diferentes coeficientes, mesmo quando o novo preditor é contínuo. Então, isso tem algo a ver com a ortogonalidade das variáveis ​​fictícias que criei no meu exemplo 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))

insira a descrição da imagem aqui

meia passagem
fonte
6
(+1) Boa pergunta. Vou tentar escrever algo para você em algumas horas. Enquanto isso, você pode tentar ver o que acontece com vários preditores categóricos que não são ortogonais (dica!).
cardeal
1
Intriga! Eu certamente tentarei aquilo. E obrigado, aguardo ansiosamente sua resposta.
half-pass
Desculpe @ meia passagem; Não esqueci disso e tentarei fazer algo dentro de um dia. (Tenho meia resposta em conjunto, mas fui afastado por outras tarefas.) Espero que a recompensa também atraia outras respostas. Felicidades. :-)
cardeal
Não se preocupe, @ cardinal! Eu sei que você e todos os outros gurus incríveis aqui têm vidas fora CV :)
meia-pass

Respostas:

29

Você descobriu uma propriedade íntima, mas genérica, dos GLMs, ajustada pela máxima probabilidade . O resultado desaparece quando se considera o caso mais simples de todos: ajustar um único parâmetro a uma única observação!

Uma resposta frase : Se tudo o que interessa é encaixe meios separados para subconjuntos disjuntos de amostra, em seguida, MLG irá sempre originar μ J = ˉ y j para cada subconjunto j , de modo que a estrutura de erro real e parametrização da densidade tanto tornar irrelevante para a estimativa (ponto)!μ^j=y¯jj

Um pouco mais : Ajustar fatores categóricos ortogonais à máxima probabilidade é equivalente a ajustar meios separados para separar subconjuntos de nossa amostra, portanto, isso explica por que Poisson e GLMs binomiais negativos produzem as mesmas estimativas de parâmetros. De fato, o mesmo acontece se usamos regressão de Poisson, negbin, Gaussian, Gaussian inversa ou Gamma (veja abaixo). No caso de Poisson e negbin, a função de link padrão é o link de , mas esse é um arenque vermelho; Embora isso produza as mesmas estimativas brutas de parâmetro, veremos abaixo que essa propriedade realmente não tem nada a ver com a função de link.registro

Quando estamos interessados ​​em uma parametrização com mais estrutura, ou que depende de preditores contínuos, a estrutura de erro assumida se torna relevante devido à relação de variância média da distribuição no que se refere aos parâmetros e à função não linear usada para modelar a condicional significa.

GLMs e famílias de dispersão exponencial: curso intensivo

Uma família de dispersão exponencial na forma natural é aquela que a densidade logarítmica é da forma

registrof(y;θ,ν)=θy-b(θ)ν+uma(y,ν).

Aqui é o parâmetro natural e ν é o parâmetro de dispersão . Se v fosse conhecido, seria apenas uma família exponencial de um parâmetro padrão. Todos os GLMs considerados abaixo assumem um modelo de erro dessa família.θνν

Considere uma amostra de uma única observação dessa família. Se nos encaixamos por máxima verossimilhança, temos que y = b ' ( θ ) , independentemente do valor de ν . Este estende-se facilmente ao caso de uma amostra iid uma vez que as probabilidades de log adicionar, obtendo-se ˉ y = b ' ( θ ) .θy=b(θ^)νy¯=b(θ^)

Mas também sabemos, devido à boa regularidade da densidade do log em função de , que θ Então, de fato b ( θ ) = E Y = μ .

θEregistrof(Y;θ,ν)=Eθregistrof(Y;θ,ν)=0 0.
b(θ)=EY=μ

Uma vez que as estimativas de probabilidade máxima são invariante sob transformações, isto significa que para esta família de densidades.y¯=μ^

μEuμEu=g-1(xEuTβ)gxEujμEu=g(βj)βj

O que há de tão diferente nos preditores contínuos?

Quando os preditores são contínuos ou são categóricos, mas não podem ser reduzidos a uma forma ortogonal, a probabilidade não é mais fatorada em termos individuais com uma média separada, dependendo de um parâmetro separado. Neste ponto, a estrutura de erro e função de ligação não entram em jogo.

Eu=1n(yEu-μEu)xEujσEu2μEuλEu=0 0,
j=1,...,pλEu=xEuTββνμEu=g(λEu)=g(xEuTβ)σEu2

Dessa forma, a função de link e o modelo de erro assumido tornam-se relevantes para a estimativa.

Exemplo: o modelo de erro (quase) não importa

k=6

registro

A partir da tabela, podemos ver que as estimativas dos parâmetros são idênticas , embora alguns desses GLMs sejam para dados discretos e outros sejam contínuos, e alguns sejam para dados não negativos, enquanto outros não.

      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

0 0

Exemplo: a função de link (quase) não importa

registroregistro(β^)registro(β^2)β^

> 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

A ressalva no cabeçalho refere-se simplesmente ao fato de que as estimativas brutas variarão com a função de link, mas as estimativas implícitas dos parâmetros médios não.

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))
cardeal
fonte
+1, clara e abrangente, uma resposta melhor do que a que eu dei.
Mpr5
@mpr, eu preferiria pensar nisso como complementar ao seu. Fiquei muito satisfeito quando vi seu post; você descreveu clara e concisamente (e mostrou) o que está acontecendo. Às vezes, minhas postagens demoram um pouco para serem lidas; Temo que este seja um deles.
cardeal
Vocês são incríveis. Muito obrigado por me explicar isso com tanta clareza e rigor. Eu preciso gastar um pouco mais de tempo digerindo agora :)
half-pass
@cardinal De onde você tirou os 2 " family=negative.binomial(theta=2)"?
timothy.s.lau
23

Para ver o que está acontecendo aqui, é útil primeiro fazer a regressão sem a interceptação, pois uma interceptação em uma regressão categórica com apenas um preditor não tem sentido:

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

Como Poisson e regressões binomiais negativas especificam o log do parâmetro da média, para a regressão categórica, exponenciando os coeficientes, você obterá o parâmetro da média real para cada categoria:

>  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

Estes parâmetros correspondem às médias reais sobre os diferentes valores de categoria:

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

λ

λ

λθθλ

eu(X,λ,θ)=(θλ+θ)θΓ(θ+xEu)xEu!Γ(θ)(λλ+θ)xEuregistroeu(X,λ,θ)=θ(registroθ-registro(λ+θ))+xEu(registroλ-registro(λ+θ))+registro(Γ(θ+xEu)xEu!Γ(θ))ddλregistroeu(X,λ,θ)=xEuλ-θ+xEuλ+θ=n(x¯λ-x¯+θλ+θ),
λ=x¯

registro(λ)λ

λλ251055=25λ11

mpr
fonte
5
(+1) Boa resposta. Uma coisa que eu acho que poderia ser explicada aqui de maneira mais explícita é que isso realmente não tem nada a ver com qualquer relação entre o Poisson e o binômio negativo, mas com fatos mais básicos sobre a adaptação de GLMs por máxima probabilidade.
cardeal
Bom ponto. A única coisa real que Poisson e o binômio negativo têm a ver com isso são que eles são especificados pelo log do parâmetro mean. Por exemplo, se você fizesse uma regressão ordinária de mínimos quadrados, os coeficientes seriam nominalmente diferentes, pois o parâmetro teria a média real e não o log.
Mpr5 #
2
É verdade, mas acho que vai um pouco além disso: ajuste qualquer GLM usando qualquer função de link (por ML e com ressalvas muito pequenas) e, como os meios ajustados corresponderão aos da amostra, as estimativas dos parâmetros serão idênticas às não-lineares transformação entre diferentes links. O modelo de erro específico é irrelevante além do fato de ser proveniente de uma família de dispersão exponencial.
cardeal
Este é um bom ponto que não cobri. Abordei esta questão do ponto de vista mais geral da estimativa de ML, em vez de especificamente GLMs. Existem muitos modelos de ocorrência natural em que o ML produzirá meios ajustados que diferem dos meios da amostra (por exemplo, lognormal). Mas para GLMs, sua observação leva a uma explicação mais concisa e mais geral.
Mpr5 #
2
@ Meia-passagem: ajuste todos os seus modelos categóricos ortogonais y~X+0e tente novamente. :-)
cardeal