Politikwissenschaft

 

Webseite durchsuchen

Poisson, Quasi Poisson oder Negativ Binomial?

Benjamin Schlegel | 10. Oktober 2016

PDF

Dieser Artikel zeigt, wie man in R herausfinden kann, ob ein Poisson, Quasi Poisson oder Negativ Binomial Modell das geeignetste ist.

Zuerst werden die Daten aus der Bibliothek COUNT und anschliessend recodiert.


library(COUNT)
data("affairs")
affairs$kids = factor(affairs$kids)
relig.matrix = as.matrix(affairs[,8:12])
affairs$religion = factor(relig.matrix%*%(1:ncol(relig.matrix)), 
    labels = c("anti religious","not religious","slightly religious",
               " somewhat religious","very religious")) 

Nun werden die drei Modelle geschätzt.


model.poisson = glm(naffairs ~ kids + religion, data=affairs, family=poisson)
model.quasi = glm(naffairs ~ kids + religion, data=affairs, family=quasipoisson)
model.nb = glm.nb(naffairs ~ kids + religion, data=affairs)

Anschliessend kann mit dem likelihood ratio test für Überdispersion getestet werden, ob der Durchschnitt genug Nahe bei der Varianz liegt, um ein Poisson Modell zu schätzen. Der Test kann mit der Funktion odTest aus der Bibliothek pscl durchgeführt werden. Als Parameter wird das Negativ Binomial Modell angegeben.


library(pscl)
odTest(model.nb)

H0 besagt, dass keine Überdispersion vorhanden ist. Der Test zeigt uns, dass H0 verworfen werden kann (p<0.05). Deshalb sollte das Poisson Modell nicht verwendet werden.


Likelihood ratio test of H0: Poisson, as restricted NB model:
n.b., the distribution of the test-statistic under H0 is non-standard
e.g., see help(odTest) for details/references

Critical value of test statistic at the alpha= 0.05 level: 2.7055 
Chi-Square Test Statistic =  1731.7461 p-value = < 2.2e-16 

Nun ist klar, dass ein Quasi Poisson oder Negativ Binomial geschätzt werden muss. Welches der beiden besser ist, lässt sich aber nicht mit einem Test herausfinden. Die beste Möglichkeit zwischen den beiden Modellen zu unterscheiden ist ein Durchschnitt-Varianz-Beziehungsgrafik. Damit kann von Auge entscheiden werden, welches Modell die Daten besser abbildet. Beim Quasi Poisson ist die Beziehung zwischen Varianz und Durchschnitt linear, beim Negativ Binomial hingegen hängt die Beziehung von einem weiteren Parameter ab. Mit folgender Funktion kann die Grafik gezeichnet werden. Als Parameter werden die Modelle Poisson und Negativ Binomial gebraucht.


mean.var.plot = function(model.poisson,model.nb){
  xb = predict(model.nb)
  g = cut(xb, breaks=unique(quantile(xb,seq(0,1,0.1))))
  m = tapply(model.poisson$y, g, mean)
  v = tapply(model.poisson$y, g, var)
  pr <- residuals(model.poisson,"pearson")
  phi <- sum(pr^2)/df.residual(model.poisson)
  x = seq(min(m),max(m),length.out = 500)
  line.data = data.frame(x=rep(x,2),y=c(x*phi,x*(1+x/model.nb$theta)),
    model=c(rep("Q. Poisson",length(x)),rep("Neg. Binomial",length(x))))
  library(ggplot2)
  ggplot() + geom_point(aes(x=m,y=v)) + 
    geom_line(aes(x=x,y=y,linetype=model),data=line.data) + 
    theme_bw() + theme(panel.background = element_rect(rgb(.95,.95,.95))) +
    ylab("variance") + xlab("mean") +
    scale_linetype_manual(values = c("solid","dashed")) +
    ggtitle("Mean-Variance Relationship") 
}

Nachdem die Funktion eingelesen ist, kann sie verwendet werden.


mean.var.plot(model.poisson, model.nb)

Es ist klar erkennbar, dass das Quasi Poisson in diesem Fall dem Negativ Binomial Modell vorgezogen werden soll.

Durchschitt-Varianz-Beziehung