Warum wird das Quasi-Poisson in GLM nicht als Sonderfall eines negativen Binomials behandelt?

21

Ich versuche, verallgemeinerte lineare Modelle an einige Sätze von Zähldaten anzupassen, die möglicherweise überdispers sind oder nicht. Die beiden hier geltenden kanonischen Verteilungen sind das Poisson- und das Negative Binomial (Negbin) mit EV und Varianzμ

VarP=μ

VarNB=μ+μ2θ

in denen R montiert werden unter Verwendung von, glm(..,family=poisson)und glm.nb(...)sind. Es gibt auch die quasipoissonFamilie, die nach meinem Verständnis ein angepasster Poisson mit dem gleichen EV und der gleichen Varianz ist

VarQP=ϕμ ,

dh irgendwo zwischen Poisson und Negbin fallen. Das Hauptproblem bei der Quasipoisson-Familie ist, dass es keine entsprechende Wahrscheinlichkeit dafür gibt und daher viele äußerst nützliche statistische Tests und Anpassungsmaße (AIC, LR usw.) nicht verfügbar sind.

Wenn Sie die QP- und Negbin-Varianzen vergleichen, stellen Sie möglicherweise fest, dass Sie sie gleichsetzen können, indem Sie . Wenn Sie diese Logik fortsetzen, können Sie versuchen, die Quasipoisson-Verteilung als Sonderfall des Negbin auszudrücken:ϕ=1+μθ

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

dh ein Negbin mit thgr; linear abhängig von mgr ; . Ich habe versucht, diese Idee zu verifizieren, indem ich eine zufällige Folge von Zahlen gemäß der obigen Formel generiert und sie mit Folgendem ausstattet :μθμ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

Beide Anpassungen reproduzieren die Parameter, und die Quasipoisson liefert eine 'vernünftige' Schätzung für . Wir können jetzt auch einen AIC-Wert für die Quasipoisson definieren:ϕ

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

(Ich musste den angepassten Wert von manuell kopieren , da ich ihn im Objekt nicht finden konnte. )ϕsummary(glmQP)glmQP

Da , würde dies anzeigen, dass die Quasipoisson, nicht überraschend, die bessere Anpassung ist; also zumindest , was es tun sollte, und daher könnte es eine vernünftige Definition für den AIC (und damit auch die Wahrscheinlichkeit) einer Quasipoisson sein. Die großen Fragen, die ich noch habe, sind dann A I C Q PEINichCQ.P<EINichCNBEINichCQ.P

  1. Ist diese Idee sinnvoll? Basiert meine Überprüfung auf Zirkelschluss?
  2. Die Hauptfrage für alle, die etwas „erfinden“, was einem etablierten Thema zu fehlen scheint: Wenn diese Idee Sinn macht, warum ist sie nicht bereits in implementiert glm?

Bearbeiten: Figur hinzugefügt

glm passt und + -1 Sigma-Bänder

user28400
quelle
1
(+1) Willkommen bei Cross Validated! Und danke für eine ausgezeichnete Frage (obwohl ein paar Kommentare im Code für Leute, die R nicht benutzen, hilfreich sein könnten). Ich glaube, Sie haben das NB1-Modell möglicherweise neu erfunden (obwohl ich es noch nicht im Detail verfolgt habe). Beachten Sie auch, dass es keine Quasi-Poisson- Verteilung gibt - weshalb es keine Wahrscheinlichkeit oder AIC gibt - es bezieht sich nur auf eine Art der Anpassung von Mitteln und Varianzen.
Scortchi - wieder einzusetzen Monica
2
Vielen Dank! Ich habe in der Zwischenzeit einige Kommentare hinzugefügt, ich hoffe das klärt die Dinge auf. Ich verstehe, dass die Quasi-Poisson-Verteilung per se nicht existiert - was ich wirklich herausfinden wollte, ist, warum QP überhaupt eine Sache ist, wenn man bedenkt, dass die NB1-Verteilung existiert und keines der Quasi-Probleme des QP hat (siehe Achims Antwort für eine offensichtliche Lösung).
User28400
1
@Scortchi --- tatsächlich gibt es eine solche Verteilung ... Wenn und , dann ist eine Exponentialfamilie mit dem Mittelwert und der Varianz . Wenn . Es ist nicht unbedingt für Zähldaten geeignet (außer als Annäherung), da es auf . Y = K X Y μ = k λ k μ k 1 0 , k , 2 k , . . .XPois(λ)Y.=kXY.μ=kλkμk10,k,2k,...
Glen_b
1
@ Glen_b: Nennen die Leute das wirklich Quasi-Poisson? Auf jeden Fall ist es eine gute Illustration - wenn Sie ein "QuasiPoisson" -Modell verwenden, nehmen Sie nicht wirklich diese Verteilung oder die NB1 oder irgendeine andere Beziehung zwischen Mittelwert und Varianz an, die Ihre Schätzungen der Koeffizienten und ihrer Standardfehler ermöglicht besser, wenn die Stichprobe größer wird.
Scortchi
1
@Scortchi Es ist die einzige exponentielle Familienverteilung, die die Annahmen des Quasi-Poisson erfüllt, also irgendwie - gelegentlich habe ich gesehen, dass Leute darauf hinweisen, dass es die Verteilung ist, die die Annahme impliziert. Wenn Leute es benutzen, haben sie natürlich fast nie die Absicht, dass ihre Daten von dieser spezifischen Verteilung stammen - es ist nur als grobe Beschreibung der Beziehung zwischen Mittelwert und Varianz gedacht. (Bei einigen Versicherungsanwendungen kann es unter sehr einfachen Voraussetzungen sinnvoll sein - die Gesamtkosten für Versicherungsfälle, bei denen die Anzahl der Versicherungsfälle Poisson beträgt und die Kosten pro
Versicherungsfall praktisch

Antworten:

24

Das Quasi-Poisson ist kein Full Maximum Likelihood (ML) -Modell, sondern ein Quasi-ML-Modell. Verwenden Sie einfach die Schätzfunktion (oder Bewertungsfunktion) aus dem Poisson-Modell, um die Koeffizienten zu schätzen, und verwenden Sie dann eine bestimmte Varianzfunktion, um geeignete Standardfehler (oder vielmehr eine vollständige Kovarianzmatrix) zu erhalten, um eine Inferenz durchzuführen. Daher glm()nicht liefern und logLik()oder AIC()hier usw.

Wie Sie richtig hervorheben, kann ein Modell mit derselben Erwartungs- und in das negative Binomial-Framework (NB) eingebettet werden, wenn der sizeParameter zusammen mit der Erwartung variiert . In der Literatur wird dies typischerweise als NB1-Parametrisierung bezeichnet. Siehe beispielsweise das Cameron & Trivedi-Buch (Regressionsanalyse von Zähldaten) oder "Analysis of Microdata" von Winkelmann & Boes.μ iθichμi

Wenn es keine Regressoren gibt (nur einen Abschnitt), stimmen die NB1-Parametrisierung und die NB2-Parametrisierung, die von verwendet werden MASS, glm.nb()überein. Bei Regressoren unterscheiden sie sich. In der statistischen Literatur wird die NB2-Parametrisierung häufiger verwendet, einige Softwarepakete bieten jedoch auch die NB1-Version an. Beispielsweise können Sie in R das gamlssPaket verwenden, um zu tun gamlss(y ~ x, family = NBII). Beachten Sie, dass etwas verwirrend gamlssverwendet NBIfür die NB2 Parametrisierung und NBIIfür NB1. (Fachsprache und Terminologie sind jedoch nicht in allen Communities einheitlich.)

Dann könnten Sie sich natürlich fragen, warum Sie Quasi-Poisson verwenden, wenn NB1 verfügbar ist. Es gibt noch einen subtilen Unterschied: Der erstere verwendet Quasi-ML und erhält die Schätzung aus der Dispersion der quadratischen Abweichungs- (oder Pearson-) Residuen. Letzterer verwendet volle ML. In der Praxis ist der Unterschied oft nicht groß, aber die Gründe für die Verwendung beider Modelle sind leicht unterschiedlich.

Achim Zeileis
quelle
1
Vielen Dank! Sehr hilfreiche Antwort, ich experimentiere gamlssjetzt damit und es sieht so aus, als ob es genau das ist, was ich brauchte. Könnten Sie die Gründe für die Verwendung von Quasi-Likelihood im Vergleich zu Full ML erläutern?
User28400
2
Sie nehmen weniger an: Sie nehmen nur (1) eine logarithmisch-lineare Beziehung zwischen der Erwartung und den Regressoren an (2) eine lineare Beziehung zwischen Varianz und Erwartung. Der Rest der Wahrscheinlichkeit ist völlig unbestimmt. Als Alternative zu (2) setzen die Praktiker manchmal sogenannte "robuste" Sandwich-Standardfehler ein, die allgemeinere Heteroskedastizitätsmuster ermöglichen würden. Natürlich könnte man den NB1 auch mit Sandwich-Standardfehlern einsetzen ... Ein paar weitere Kommentare sind in unserem vignette("countreg", package = "pscl").
Achim Zeileis,