Pourquoi le quasi-Poisson en GLM n'est-il pas traité comme un cas particulier de binôme négatif?

21

J'essaie d'adapter les modèles linéaires généralisés à certains ensembles de données de comptage qui pourraient ou non être sur-dispersés. Les deux distributions canoniques qui s'appliquent ici sont le binôme de Poisson et négatif (Negbin), avec EV et varianceμ

VunerP=μ

VunerNB=μ+μ2θ

qui peuvent être montés dans R en utilisant glm(..,family=poisson)et glm.nb(...), respectivement. Il y a aussi la quasipoissonfamille, qui selon moi est un Poisson ajusté avec le même EV et la même variance

VunerQP=ϕμ ,

c'est-à-dire se situant quelque part entre Poisson et Negbin. Le principal problème avec la famille des quasipoisson est qu'il n'y a pas de probabilité correspondante pour elle, et donc de nombreux tests statistiques et mesures d'ajustement extrêmement utiles (AIC, LR, etc.) ne sont pas disponibles.

Si vous comparez les variances QP et Negbin, vous remarquerez peut-être que vous pouvez les égaliser en mettant . En poursuivant dans cette logique, vous pourriez essayer d'exprimer la distribution quasipoisson comme un cas particulier du Negbin:ϕ=1+μθ

QP(μ,ϕ)=NB(μ,θ=μϕ-1) ,

c'est-à-dire un Negbin avec dépendant linéairement de . J'ai essayé de vérifier cette idée en générant une séquence aléatoire de nombres selon la formule ci-dessus et en l'adaptant avec :μθμglm

#fix parameters

phi = 3
a = 1/50
b = 3
x = 1:100

#generating points according to an exp-linear curve
#this way the default log-link recovers the same parameters for comparison

mu = exp(a*x+b) 
y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator

#fit a generalized linear model y = f(x)  
glmQP = glm(y~x, family=quasipoisson) #quasipoisson
glmNB = glm.nb(y~x) #negative binomial

> glmQP

Call:  glm(formula = y ~ x, family = quasipoisson)

Coefficients:
(Intercept)            x  
    3.11257      0.01854  
(Dispersion parameter for quasipoisson family taken to be 3.613573)

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      2097 
Residual Deviance: 356.8    AIC: NA

> glmNB

Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)

Coefficients:
(Intercept)            x  
    3.10182      0.01873  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      578.1 
Residual Deviance: 107.8    AIC: 824.7

Les deux ajustements reproduisent les paramètres et le quasipoisson donne une estimation «raisonnable» pour . On peut maintenant aussi définir une valeur AIC pour le quasipoisson:ϕ

df = 3 # three model parameters: a,b, and phi
phi.fit = 3.613573 #fitted phi value copied from summary(glmQP)
mu.fit = glmQP$fitted.values 

#dnbinom = negbin density, log=T returns log probabilities
AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T))
> AIC
[1] 819.329

(J'ai dû copier manuellement la valeur ajustée depuis , car je ne pouvais pas la trouver dans l' objet)ϕsummary(glmQP)glmQP

Puisque , cela indiquerait que le quasipoisson est, sans surprise, le meilleur ajustement; donc au moins fait ce qu'il devrait faire, et donc ce pourrait être une définition raisonnable pour l'AIC (et par extension, la vraisemblance) d'un quasipoisson. Les grandes questions qui me restent sont alors A I C Q PUNEjeCQP<UNEjeCNBUNEjeCQP

  1. Cette idée a-t-elle un sens? Ma vérification est-elle basée sur un raisonnement circulaire?
  2. La principale question pour quiconque «invente» quelque chose qui semble manquer dans un sujet bien établi: si cette idée a du sens, pourquoi n'est-elle pas déjà mise en œuvre glm?

Modifier: figure ajoutée

ajustements glm et bandes +1 sigma

user28400
la source
1
(+1) Bienvenue sur Cross Validated! Et merci pour une excellente question (même si quelques commentaires dans le code peuvent être bien pour les personnes qui n'utilisent pas R). Je pense que vous avez peut-être réinventé le modèle NB1 (même si je ne l'ai pas encore suivi en détail). Notez également qu'il n'y a pas de distribution quasi-Poisson - c'est pourquoi il n'y a pas de probabilité ou d'AIC - cela se réfère simplement à un moyen d'ajuster les moyennes et les variances.
Scortchi - Réintégrer Monica
2
Merci! J'ai ajouté quelques commentaires en attendant, j'espère que cela clarifie les choses. Je comprends que la distribution quasi-Poisson n'existe pas en soi - ce que j'essayais vraiment de comprendre, c'est pourquoi QP est même une chose du tout, étant donné que la distribution NB1 existe et n'a aucun des quasi-problèmes de QP (voir la réponse d'Achims pour une résolution apparente).
user28400
1
XPois(λ)Oui=kXOuiμ=kλkμk10,k,2k,...
1
@Glen_b: Les gens appellent-ils vraiment ça le quasi-Poisson? Dans tous les cas, c'est une belle illustration - lorsque vous utilisez un modèle "quasiPoisson", vous ne supposez pas vraiment que la distribution, ou le NB1, ou tout autre, juste une relation entre la moyenne et la variance qui fait vos estimations des coefficients et leurs erreurs standard mieux que l'échantillon devient plus grand.
Scortchi - Réintégrer Monica
1
@Scortchi C'est la seule distribution exponentielle de la famille qui satisfait les hypothèses du quasi-Poisson, donc en quelque sorte - à l'occasion, j'ai vu des gens souligner que c'est la distribution qu'implique l'hypothèse. Bien sûr, lorsque les gens l'utilisent, ils n'ont presque * jamais l'intention que leurs données proviennent de cette distribution spécifique - il s'agit simplement d'une description approximative de la façon dont leur moyenne et leur variance sont liées. (Cela peut avoir du sens dans des hypothèses très simples dans certaines applications d'assurance - coût total des sinistres, où le nombre de sinistres est Poisson et le coût par sinistre est effectivement constant.)
Glen_b -Reinstate Monica

Réponses:

24

Le quasi-Poisson n'est pas un modèle à maximum de vraisemblance maximale (ML) mais un modèle quasi-ML. Vous utilisez simplement la fonction d'estimation (ou la fonction de score) du modèle de Poisson pour estimer les coefficients, puis utilisez une certaine fonction de variance pour obtenir des erreurs standard appropriées (ou plutôt une matrice de covariance complète) pour effectuer l'inférence. Par conséquent, glm()ne fournit pas et logLik()ou AIC()ici etc.

sizeθjeμje

S'il n'y a pas régresseurs (juste une interception) la paramétrisation NB1 et la paramétrisation NB2 employé par MASSde » glm.nb()la coïncidence. Avec les régresseurs, ils diffèrent. Dans la littérature statistique, le paramétrage NB2 est plus fréquemment utilisé mais certains progiciels proposent également la version NB1. Par exemple, dans R, vous pouvez utiliser le gamlsspackage pour le faire gamlss(y ~ x, family = NBII). Notez que quelque peu confuse gamlssutilise NBIpour le paramétrage NB2 et NBIIpour NB1. (Mais le jargon et la terminologie ne sont pas unifiés dans toutes les communautés.)

Vous pourriez alors vous demander, bien sûr, pourquoi utiliser du quasi-Poisson si NB1 est disponible? Il y a encore une différence subtile: le premier utilise du quasi-ML et obtient l'estimation à partir de la dispersion à partir des résidus de déviance au carré (ou Pearson). Ce dernier utilise le ML complet. Dans la pratique, la différence n'est souvent pas grande, mais les motivations à utiliser l'un ou l'autre modèle sont légèrement différentes.

Achim Zeileis
la source
1
Merci! Réponse très utile, j'expérimente gamlssmaintenant et il semble que c'est exactement ce dont j'avais besoin. Pourriez-vous nous expliquer les motivations de l'utilisation de la quasi-vraisemblance par rapport au ML complet?
user28400
2
Vous supposez moins: vous supposez simplement (1) une relation log-linéaire entre l'espérance et les régresseurs (2) une relation linéaire entre la variance et l'espérance. Le reste de la probabilité n'est pas précisé. Comme alternative à (2), les praticiens emploient parfois des erreurs standard dites "robustes" en sandwich qui permettraient des schémas d'hétéroskédasticité plus généraux. Bien sûr, on pourrait également utiliser le NB1 avec des erreurs standard sandwich ... Quelques autres commentaires sont dans notre vignette("countreg", package = "pscl").
Achim Zeileis