Quand les régressions de Poisson et binomiales négatives correspondent-elles aux mêmes coefficients?

38

J'ai remarqué que dans R, les régressions de Poisson et binomiales négatives (NB) semblent toujours correspondre aux mêmes coefficients pour les prédicteurs catégoriels, mais non continus.

Par exemple, voici une régression avec un prédicteur catégorique:

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

entrez la description de l'image ici

Voici un exemple avec un prédicteur continu, où les lois de Poisson et NB correspondent à des coefficients différents:

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

entrez la description de l'image ici

(Bien sûr, ce ne sont pas des données de comptage, et les modèles ne sont pas significatifs ...)

Ensuite, je recodifie le prédicteur en un facteur, et les deux modèles correspondent à nouveau aux mêmes coefficients:

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

entrez la description de l'image ici

Cependant, la régression binomiale négative de Joseph Hilbe donne un exemple (6.3, pages 118-119) dans lequel un prédicteur catégorique, le sexe, est ajusté avec des coefficients légèrement différents selon les lois de Poisson ( b=0.883 ) et NB ( b=0.881 ). Il déclare: «Les ratios de taux d'incidence entre les modèles de Poisson et NB sont assez similaires. Ceci est pas surprenant étant donné la proximité de α [correspondant à 1/θ dans R] à zéro « .

Je comprends cela, mais dans les exemples ci-dessus, summary(rs2)nous dit que est estimé à 9.16 et 7.93 respectivement.θ

Alors, pourquoi les coefficients sont-ils exactement les mêmes? Et pourquoi uniquement pour les prédicteurs catégoriques?


Modifier # 1

Voici un exemple avec deux prédicteurs non orthogonaux. En effet, les coefficients ne sont plus les mêmes:

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

entrez la description de l'image ici

De plus, l'ajout d'un autre prédicteur force les modèles à s'adapter à différents coefficients, même lorsque le nouveau prédicteur est continu. Donc, cela a quelque chose à voir avec l'orthogonalité des variables factices que j'ai créées dans mon exemple initial?

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

entrez la description de l'image ici

demi-passe
la source
6
(+1) Belle question. Je vais essayer d'écrire quelque chose pour vous dans quelques heures. En attendant, vous pouvez essayer de voir ce qui se passe avec plusieurs prédicteurs catégoriels qui ne sont pas orthogonaux (indice!).
cardinal
1
Intrigue! Je vais certainement essayer ça. Et merci, j'attends avec impatience votre réponse.
demi-passage du
Désolé @ demi-passe; Je n’ai pas oublié cela et je vais essayer de faire quelque chose en moins d’une journée. (J'ai la moitié d'une réponse, mais je me suis laissé entraîner par d'autres tâches.) J'espère que la prime attirera également d'autres réponses. À votre santé. :-)
cardinal
Pas de soucis, cardinal! Je sais que vous et tous les autres gourous étonnants ici avez des vies en dehors de CV :)
demi-réussite du

Réponses:

29

Vous avez découvert une propriété intime, mais générique, de GLM qui s'adapte parfaitement . Le résultat disparaît dès que l'on considère le cas le plus simple: ajuster un seul paramètre à une seule observation!

Une réponse de la phrase : Si tout ce que nous préoccupons convient des moyens séparés pour des sous - ensembles disjoints de notre échantillon, puis GLM produira toujours μ j = ˉ y pour chaque sousensemble j ,sortela structure d'erreur réelle et paramétrisation de la densitéfois devenu sans objet à l'estimation (point)!μ^j=y¯jj

Un peu plus : ajuster les facteurs catégoriels orthogonaux au maximum de vraisemblance équivaut à ajuster des moyens distincts à des sous-ensembles disjoints de notre échantillon, ce qui explique pourquoi les GLM de Poisson et binomiaux négatifs génèrent les mêmes estimations de paramètres. En fait, il en va de même si nous utilisons une régression de Poisson, Negbin, Gaussienne, Gaussienne Inverse ou Gamma (voir ci-dessous). Dans le cas de Poisson et negbin, la fonction de liaison par défaut est la liaison de , mais il s’agit d’un faux fil rouge. Bien que cela produise les mêmes estimations de paramètres bruts , nous verrons ci-dessous que cette propriété n'a vraiment rien à voir avec la fonction de liaison.bûche

Lorsque nous nous intéressons à une paramétrisation avec plus de structure, ou qui dépend de prédicteurs continus, alors la structure d'erreur supposée devient pertinente en raison de la relation moyenne-variance de la distribution en relation avec les paramètres et de la fonction non linéaire utilisée pour modéliser le conditionnel veux dire.

GLM et familles de dispersion exponentielle: cours intensif

Une famille de dispersion exponentielle sous forme naturelle est telle que la densité de log est de la forme

logf(y;θ,ν)=θyb(θ)ν+a(y,ν).

Ici est le paramètre naturel et ν le paramètre de dispersionθν . Si était connu, il ne s'agirait que d'une famille exponentielle standard à un paramètre. Tous les GLM considérés ci-dessous supposent un modèle d'erreur de cette famille.ν

Prenons un exemple d’une seule observation de cette famille. Si nous nous situons par maximum de vraisemblance, on obtient que y = b ' ( θ ) , quelle que soit la valeur de ν . Cela va facilement au cas d'un échantillon IID depuis les vraisemblances log ajouter, ce qui donne ˉ y = b ' ( θ ) .θy=b(θ^)νy¯=b(θ^)

Mais, nous savons aussi, en raison de la belle régularité de la densité du journal en fonction de , que θ Donc, en fait, b ( θ ) = E Y

θElogf(Y;θ,ν)=Eθlogf(Y;θ,ν)=0.
.b(θ)=EY=μ

Étant donné que les estimations du maximum de vraisemblance sont invariantes par transformation, cela signifie que y¯=μ^ pour cette famille de densités.

Maintenant, dans un GLM, nous modélisons tant que μ i = g - 1 (μigest la fonction de liaison. Mais si x i est un vecteur de tous les zéros à l'exception d'un seul 1 en positionj, alors µ i =g( β j ). La vraisemblance du GLM est alors factorisée en fonction de la βμi=g1(xiTβ)gxijμi=g(βj) et on procède comme ci-dessus. C'est précisément le cas des facteurs orthogonaux.βj

En quoi les prédicteurs continus sont-ils si différents?

Lorsque les prédicteurs sont continus ou catégoriques, mais ne peuvent pas être réduits à une forme orthogonale, la probabilité ne tient plus compte de termes individuels avec une moyenne distincte dépendant d'un paramètre distinct. A ce stade, la structure d'erreur et la fonction de lien n'entrent en jeu.

Si une manivelle par le biais de l'algèbre (fastidieux), les équations de vraisemblance deviennent Pour tout j = 1 , ... , p λ i = x T i β . Ici, les β et ν

i=1n(yiμi)xijσi2μiλi=0,
j=1,,pλje=XjeTββν paramètres entrent implicitement via la relation de liaison et la variance σ 2 iμje=g(λje)=g(XjeTβ)σje2 .

De cette manière, la fonction de lien et le modèle d'erreur supposé deviennent pertinents pour l'estimation.

Exemple: le modèle d'erreur (presque) n'a pas d'importance

Dans l'exemple ci-dessous, nous générons des données aléatoires binomiales négatives en fonction de trois facteurs catégoriels. Chaque observation provient d’une seule catégorie et du même paramètre de dispersion ( k=6 utilise ).

Nous nous sommes ensuite adaptés à ces données en utilisant cinq GLM différents, chacun avec un lien de : ( a ) binôme négatif, ( b ) Poisson, ( c ) gaussien, ( d ) gaussien inverse et (bûche e ) GLM. Ce sont tous des exemples de familles de dispersion exponentielle.

Dans le tableau, nous pouvons constater que les estimations de paramètres sont identiques , même si certains de ces GLM concernent des données discrètes et d’autres des données continues, et d’autres des données non négatives, d’autres non.

      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 mise en garde dans l'en-tête vient du fait que la procédure d'ajustement échouera si les observations ne relèvent pas du domaine de la densité particulière. Par exemple, si nous avions 0 comptes générés aléatoirement dans les données ci-dessus, le GLM Gamma ne pourra pas converger car les GLM Gamma nécessitent des données strictement positives.

Exemple: la fonction de lien (presque) n'a pas d'importance

En utilisant les mêmes données, nous répétons la procédure en ajustant les données avec un GLM de Poisson avec trois fonctions de liaison différentes: ( a )bûchebûche(β^)bûche(β^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

La mise en garde dans l'en-tête fait simplement référence au fait que les estimations brutes varieront avec la fonction de lien, mais pas les estimations implicites en paramètres moyens.

Code 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))
cardinal
la source
+1, clair et complet, une réponse meilleure que celle que j'ai donnée.
Mpr
@mpr, je préférerais de beaucoup le considérer comme complémentaire du vôtre. J'ai été très heureux quand j'ai vu votre message. vous avez décrit de manière claire et concise (et montré) ce qui se passe. Mes publications sont parfois un peu longues à lire; Je crains que ce soit l'un d'entre eux.
cardinal
Vous êtes tous les deux incroyables. Merci beaucoup de m'avoir expliqué cela avec tant de clarté et de rigueur. J'ai besoin de passer plus de temps à digérer maintenant :)
demi-passage du
@ cardinal D'où avez-vous obtenu le 2 " family=negative.binomial(theta=2)"?
timothy.s.lau
23

Pour voir ce qui se passe ici, il est utile de commencer par effectuer la régression sans intercept, car une interception dans une régression catégorique avec un seul prédicteur n'a pas de sens:

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

Étant donné que les régressions de Poisson et binomiales négatives spécifient le journal du paramètre moyen, alors pour une régression catégorielle, une exponentiation des coefficients vous donnera le paramètre moyen réel pour chaque catégorie:

>  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

Ces paramètres correspondent aux moyennes réelles sur les différentes valeurs de catégorie:

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

λ

λ qui maximise la vraisemblance pour les observations de chaque catégorie. L'estimateur de vraisemblance maximum pour la distribution de Poisson est simplement la moyenne de l'échantillon, raison pour laquelle les coefficients de régression sont précisément les moyennes (journaux des) échantillons de chaque catégorie.

λθθλ

L(X,λ,θ)=Π(θλ+θ)θΓ(θ+Xje)Xje!Γ(θ)(λλ+θ)XjebûcheL(X,λ,θ)=Σθ(bûcheθ-bûche(λ+θ))+Xje(bûcheλ-bûche(λ+θ))+bûche(Γ(θ+Xje)Xje!Γ(θ))λbûcheL(X,λ,θ)=ΣXjeλ-θ+Xjeλ+θ=n(X¯λ-X¯+θλ+θ),
de sorte que le maximum est atteint lorsque λ=X¯.

La raison pour laquelle vous n'obtenez pas les mêmes coefficients pour les données continues est parce que, dans une régression continue, bûche(λ)ne sera plus une fonction constante par morceaux des variables prédictives, mais linéaire. Maximiser la fonction de vraisemblance dans ce cas ne se réduira pas à l'ajustement indépendant d'une valeurλ pour des sous-ensembles disjoints des données, mais sera plutôt un problème non trivial qui est résolu numériquement, et est susceptible de produire des résultats différents pour différentes fonctions de vraisemblance.

De même, si vous avez plusieurs prédicteurs catégoriels, malgré le fait que le modèle ajusté spécifiera finalement λ fonction constante par morceaux, en général, il n’y aura pas assez de degrés de liberté pour permettre λà déterminer indépendamment pour chaque segment constant. Par exemple, supposons que vous ayez2 prédicteurs avec 5catégories chacune. Dans ce cas, votre modèle adix degrés de liberté, alors qu'il y a 5*5=25 différentes combinaisons uniques de catégories, chacune ayant sa propre valeur ajustée de λ. Donc, en supposant que les intersections de ces catégories ne soient pas vides (ou du moins que11 d’entre eux sont non vides), le problème de la maximisation de la probabilité devient non trivial et produira généralement des résultats différents pour la distribution de Poisson par rapport au binôme négatif ou à n’importe quelle autre distribution.

mpr
la source
5
(+1) Bonne réponse. Une chose que je pense pouvoir expliquer plus explicitement ici est que cela n’a vraiment rien à voir avec une relation entre le Poisson et le binôme négatif, mais plutôt avec des faits plus fondamentaux sur l’ajustement des GLM avec le maximum de vraisemblance.
cardinal
Bon point. La seule chose réelle que Poisson et le binôme négatif ont à voir avec cela est qu'ils sont spécifiés par le log du paramètre moyen. Par exemple, si vous utilisiez une régression des moindres carrés ordinaires, les coefficients seraient nominalement différents car le paramètre correspondrait alors à la moyenne réelle plutôt qu'au journal.
Mpr
2
C'est vrai, mais je pense que cela va un peu plus loin: Ajustez n'importe quel GLM en utilisant n'importe quelle fonction de lien (par ML, et avec des mises en garde très mineures) et, comme les moyennes ajustées correspondent aux moyennes d'échantillon, les estimations de paramètres sont identiques jusqu'à la méthode non linéaire. transformation entre différents liens. Le modèle d'erreur spécifique n'est pas pertinent, outre le fait qu'il provient d'une famille de dispersion exponentielle.
cardinal
C'est un bon point que je n'ai pas couvert. J'ai abordé cette question du point de vue plus général de l'estimation de la valeur maximale, plutôt que des GLM en particulier. Il existe de nombreux modèles naturels dans lesquels ML produira des moyennes ajustées différentes de la moyenne de l'échantillon (par exemple, lognormales). Mais pour les GLM, votre observation conduit à une explication plus concise et plus générale.
Mpr
2
@ half-pass: adaptez tous vos modèles catégoriques orthogonaux y~X+0et réessayez. :-)
cardinal