Poisson-Regression

Einführung in die Poisson-Regression

Die Poisson-Regression dient dazu, zu modellieren, wie oft ein Ereignis innerhalb einer bestimmten Zeit oder eines bestimmten Gebiets eintritt. Die Zielgröße ist also eine Anzahl (engl. count): 0, 1, 2, 3 und so weiter. Um die Methode richtig anzuwenden, muss man einige zentrale Begriffe und die Wahrscheinlichkeitsverteilung verstehen, auf der die Methode beruht, nämlich die Poisson-Verteilung. Das Kapitel behandelt die Begriffe, die Poisson-Verteilung und das Modell und schließt mit einem vollständigen Beispiel in R.

Angenommen, wir untersuchen, wie viele Patienten pro Tag eine hausärztliche Praxis aufsuchen. Wir möchten wissen, welche Faktoren (Prädiktoren) die Anzahl der Besuche beeinflussen, und wir möchten vorhersagen (prädizieren) können, wie viele Besuche an einem bestimmten Tag zu erwarten sind. Jeder Besuch ist ein Ereignis (engl. event), und die Anzahl der Besuche pro Tag ist eine Anzahl. Setzt man die Anzahl in Beziehung zu der Zeit, auf die sie sich bezieht, erhält man eine Rate (Ereignisse pro Zeiteinheit): Kommen an einem Tag 10 Patienten, beträgt die Rate 10 Patienten pro Tag. Die Rate kann sich auch auf ein Gebiet beziehen, etwa die Anzahl der Personen, die ärztliche Hilfe suchen, pro Quadratkilometer.

Die Poisson-Regression kann verwendet werden, solange die Ereignisse durch einen Poisson-Prozess entstehen, also durch einen Prozess, der eine Poisson-Verteilung hervorbringt (siehe unten). Welcher Art das Ereignis ist, spielt keine Rolle. Die Methode wird eingesetzt, um die Anzahl von Todesfällen, Krankenhauseinweisungen, Infektionen, Kunden oder E-Mails zu untersuchen.

Die Poisson-Verteilung

Die Poisson-Verteilung wurde von dem französischen Mathematiker Siméon Denis Poisson beschrieben. Sie gibt die Wahrscheinlichkeit an, dass eine bestimmte Anzahl von Ereignissen innerhalb eines gegebenen Zeit- oder Raumintervalls eintritt, vorausgesetzt, die Ereignisse treten unabhängig voneinander und mit einer konstanten durchschnittlichen Rate auf. Die Verteilung hat nur einen Parameter, λ (Lambda), die erwartete Anzahl von Ereignissen im Intervall. Die Wahrscheinlichkeit, genau $k$ Ereignisse zu beobachten, beträgt

$$P(Y = k) = \frac{\lambda^{k} e^{-\lambda}}{k!}, \qquad k = 0, 1, 2, \ldots$$

Eine charakteristische Eigenschaft ist, dass Erwartungswert und Varianz beide gleich λ sind. Gerade diese Eigenschaft ist in realen Daten häufig nicht erfüllt (siehe den Abschnitt zur Überdispersion).

Die Abbildung zeigt die Poisson-Verteilung für drei Werte von λ. Die Verteilung ist diskret: Wahrscheinlichkeit gibt es nur für ganze Zahlen. Ist λ klein, liegt der Großteil der Wahrscheinlichkeit nahe null, und null Ereignisse sind häufig. Mit zunehmendem λ verschiebt sich die Verteilung nach rechts, wird breiter und zunehmend symmetrisch, und null Ereignisse werden unwahrscheinlich. Der Gipfel der Verteilung liegt bei λ oder knapp darunter (für ganzzahliges λ sind λ und λ − 1 gleich wahrscheinlich), während der Erwartungswert stets genau λ beträgt.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Poisson-Verteilung für λ = 1, 4 und 10: Wahrscheinlichkeit für jede Anzahl von Ereignissen k.",
  "width": "container",
  "height": 280,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"k": 0, "lambda": "λ = 1", "p": 0.3679}, {"k": 1, "lambda": "λ = 1", "p": 0.3679}, {"k": 2, "lambda": "λ = 1", "p": 0.1839}, {"k": 3, "lambda": "λ = 1", "p": 0.0613}, {"k": 4, "lambda": "λ = 1", "p": 0.0153}, {"k": 5, "lambda": "λ = 1", "p": 0.0031}, {"k": 6, "lambda": "λ = 1", "p": 0.0005}, {"k": 7, "lambda": "λ = 1", "p": 0.0001},
      {"k": 0, "lambda": "λ = 4", "p": 0.0183}, {"k": 1, "lambda": "λ = 4", "p": 0.0733}, {"k": 2, "lambda": "λ = 4", "p": 0.1465}, {"k": 3, "lambda": "λ = 4", "p": 0.1954}, {"k": 4, "lambda": "λ = 4", "p": 0.1954}, {"k": 5, "lambda": "λ = 4", "p": 0.1563}, {"k": 6, "lambda": "λ = 4", "p": 0.1042}, {"k": 7, "lambda": "λ = 4", "p": 0.0595}, {"k": 8, "lambda": "λ = 4", "p": 0.0298}, {"k": 9, "lambda": "λ = 4", "p": 0.0132}, {"k": 10, "lambda": "λ = 4", "p": 0.0053}, {"k": 11, "lambda": "λ = 4", "p": 0.0019}, {"k": 12, "lambda": "λ = 4", "p": 0.0006},
      {"k": 1, "lambda": "λ = 10", "p": 0.0005}, {"k": 2, "lambda": "λ = 10", "p": 0.0023}, {"k": 3, "lambda": "λ = 10", "p": 0.0076}, {"k": 4, "lambda": "λ = 10", "p": 0.0189}, {"k": 5, "lambda": "λ = 10", "p": 0.0378}, {"k": 6, "lambda": "λ = 10", "p": 0.0631}, {"k": 7, "lambda": "λ = 10", "p": 0.0901}, {"k": 8, "lambda": "λ = 10", "p": 0.1126}, {"k": 9, "lambda": "λ = 10", "p": 0.1251}, {"k": 10, "lambda": "λ = 10", "p": 0.1251}, {"k": 11, "lambda": "λ = 10", "p": 0.1137}, {"k": 12, "lambda": "λ = 10", "p": 0.0948}, {"k": 13, "lambda": "λ = 10", "p": 0.0729}, {"k": 14, "lambda": "λ = 10", "p": 0.0521}, {"k": 15, "lambda": "λ = 10", "p": 0.0347}, {"k": 16, "lambda": "λ = 10", "p": 0.0217}, {"k": 17, "lambda": "λ = 10", "p": 0.0128}, {"k": 18, "lambda": "λ = 10", "p": 0.0071}, {"k": 19, "lambda": "λ = 10", "p": 0.0037}, {"k": 20, "lambda": "λ = 10", "p": 0.0019}
    ]
  },
  "encoding": {
    "x": {"field": "k", "type": "quantitative", "title": "Anzahl der Ereignisse (k)", "scale": {"domain": [0, 20]}, "axis": {"tickMinStep": 1, "labelFontSize": 12, "titleFontSize": 13, "grid": false}},
    "y": {"field": "p", "type": "quantitative", "title": "Wahrscheinlichkeit P(Y = k)", "axis": {"labelFontSize": 12, "titleFontSize": 13, "format": ".2f"}},
    "color": {"field": "lambda", "type": "nominal", "sort": ["λ = 1", "λ = 4", "λ = 10"], "scale": {"range": ["#0e7490", "#b0442c", "#149e61"]}, "legend": {"title": null, "orient": "top", "labelFontSize": 13}}
  },
  "layer": [
    {"mark": {"type": "line", "strokeWidth": 1, "opacity": 0.35}},
    {"mark": {"type": "point", "filled": true, "size": 60}, "encoding": {"tooltip": [{"field": "lambda", "title": "Erwartete Anzahl"}, {"field": "k", "title": "Anzahl der Ereignisse"}, {"field": "p", "title": "Wahrscheinlichkeit", "format": ".4f"}]}}
  ]
}

Beispiel. Sie arbeiten in einer hausärztlichen Praxis, in die durchschnittlich 10 Patienten pro Tag kommen, also λ = 10. Die Anzahl schwankt von Tag zu Tag. Kommen die Patienten unabhängig voneinander, lässt sich mit der Poisson-Verteilung berechnen, wie wahrscheinlich es ist, dass an einem bestimmten Tag genau 5 Patienten kommen oder mehr als 5 kommen. In R geschieht das mit dpois() (Wahrscheinlichkeit für genau k) und ppois() (Wahrscheinlichkeit für k oder weniger):

# Wahrscheinlichkeit, dass genau 5 Patienten kommen, wenn die erwartete Anzahl 10 beträgt
dpois(5, lambda = 10)
#> [1] 0.03783327

# Wahrscheinlichkeit, dass mehr als 5 Patienten kommen
ppois(5, lambda = 10, lower.tail = FALSE)
#> [1] 0.932914

Genau 5 Patienten kommen also etwa an jedem 26. Tag (3,8 %), während an gut 93 % der Tage mehr als 5 Patienten kommen.

Annahmen der Poisson-Verteilung

• Die Zielgröße $k$ ist eine Anzahl von Ereignissen innerhalb eines bestimmten Intervalls und eine nichtnegative ganze Zahl (0, 1, 2, 3 und so weiter).
• Die Ereignisse treten unabhängig voneinander auf: Das Eintreten eines Ereignisses beeinflusst nicht, wann das nächste Ereignis eintritt.
• Die durchschnittliche Rate, mit der Ereignisse auftreten, ist innerhalb des Intervalls konstant.
• Zwei Ereignisse können nicht exakt gleichzeitig eintreten.

Die Poisson-Regression ist ein verallgemeinertes lineares Modell

Die Poisson-Regression ist, wie der Name schon sagt, eine Form der Regressionsanalyse. Sie wird verwendet, wenn die abhängige Variable (die Zielgröße, Y) eine Anzahl ist. Die Poisson-Regression gehört zu den verallgemeinerten linearen Modellen (engl. generalized linear models, GLM).

Was ist ein GLM?

Verallgemeinerte lineare Modelle bilden einen Rahmen, der mehrere statistische Modelle umfasst, darunter die gewöhnliche lineare Regression. In der linearen Regression wird Y direkt als lineare Funktion der Prädiktoren ($X_1, X_2, X_3, \ldots$) plus einem normalverteilten Fehlerterm modelliert. Für Anzahlen funktioniert das schlecht: Eine Anzahl kann nicht negativ sein, und die Streuung nimmt häufig mit dem Mittelwert zu. In einem GLM darf die Zielgröße stattdessen anderen Verteilungen aus der sogenannten Exponentialfamilie folgen (unter anderem der Poisson-, Binomial-, Gamma- und Normalverteilung), und der Zusammenhang zwischen den Prädiktoren und dem Erwartungswert der Zielgröße wird über eine Linkfunktion hergestellt.

Die Poisson-Regression verwendet den natürlichen Logarithmus als Linkfunktion. Modelliert wird der Logarithmus der erwarteten Anzahl, $\mu = E(Y)$, als lineare Funktion der Prädiktoren:

$$\log(\mu_i) = \beta_0 + \beta_1 X_{1i} + \beta_2 X_{2i} + \cdots + \beta_p X_{pi}$$

Anders als in der linearen Regression enthält die Gleichung keinen separaten Fehlerterm: Die zufällige Variation um $\mu_i$ wird stattdessen durch die Poisson-Verteilung beschrieben. Durch den Logarithmus ist die erwartete Anzahl stets positiv, $\mu_i = e^{\beta_0 + \beta_1 X_{1i} + \cdots}$, und die Effekte der Prädiktoren werden multiplikativ.

Interpretation der Koeffizienten: Inzidenzratenverhältnisse (IRR)

Da das Modell logarithmisch ist, werden die Koeffizienten nach Exponentiation interpretiert. $e^{\beta}$ heißt Inzidenzratenverhältnis (engl. incidence rate ratio, IRR) und gibt an, um welchen Faktor sich die erwartete Anzahl (oder die Rate) ändert, wenn der Prädiktor um eine Einheit zunimmt, beziehungsweise für eine Kategorie im Vergleich zur Referenzkategorie, bei unveränderten übrigen Prädiktoren. Ein IRR = 1 bedeutet keinen Unterschied, ein IRR = 2 eine verdoppelte Rate und ein IRR = 0,5 eine halbierte Rate.

Beispiele für Fragen, die mit der Poisson-Regression beantwortet werden

Beispiel 1. Wie beeinflussen Außentemperatur und Wetter die Anzahl der Erkältungen, die in einer hausärztlichen Praxis pro Woche diagnostiziert werden? Die Daten könnten so aussehen (die ersten 13 Wochen sind dargestellt):

Woche, Temperatur und Wetter sind Prädiktoren, mit denen die Anzahl der Erkältungen erklärt und vorhergesagt wird. Das Modell lautet

$$\log(\text{erwartete Anzahl Erkältungen}_i) = \beta_0 + \beta_1,\text{Woche}_i + \beta_2,\text{Temperatur}_i + \beta_3,\text{Wetter}_i$$

wobei $\beta_0$ der Achsenabschnitt ist (der Logarithmus der erwarteten Anzahl, wenn alle Prädiktoren null sind oder in ihrer Referenzkategorie liegen) und $\beta_1, \beta_2, \beta_3$ die Koeffizienten der Prädiktoren sind. Das Wetter ist eine kategoriale Variable und erhält in der Praxis einen Koeffizienten pro Kategorie mit Ausnahme der Referenzkategorie.

Beispiel 2. Wie viele Kunden tätigen pro Tag einen Kauf in einem Geschäft, und was beeinflusst das?

Beispiel 3. Wie viele Personen fahren mit einer bestimmten Buslinie zur Arbeit, und wie wird diese Anzahl von Wetter und Wochentag beeinflusst?

Offset: Analyse von Raten

Oft ist nicht die Anzahl an sich von Interesse, sondern die Anzahl im Verhältnis zu einer Exposition: Zeit, Fläche, Strecke oder Bevölkerungsgröße. Zwei Beispiele:

• Sie untersuchen Todesfälle und wissen, wie lange jede Person nachbeobachtet wurde (die Expositionszeit). Dann ist die Inzidenz von Interesse, also die Anzahl der Ereignisse geteilt durch die Expositionszeit. Für die Anzahl der Ereignisse pro Zeiteinheit wird häufig der Begriff event rate verwendet.
• Sie untersuchen die Anzahl der Tierarten in verschiedenen geografischen Gebieten unterschiedlicher Größe. Dann ist die Anzahl der Arten pro Fläche von Interesse.

Die Exposition geht als Offset in das Modell ein: Der Logarithmus der Exposition wird auf der rechten Seite mit einem auf 1 fixierten Koeffizienten hinzugefügt,

$$\log(\mu_i) = \log(t_i) + \beta_0 + \beta_1 X_{1i} + \cdots \quad\Longleftrightarrow\quad \log!\left(\frac{\mu_i}{t_i}\right) = \beta_0 + \beta_1 X_{1i} + \cdots$$

Das Modell beschreibt dann die Rate $\mu_i / t_i$, und $e^{\beta}$ wird zu einem Verhältnis von Raten. In R wird der Offset mit der Funktion offset() in der Formel oder mit dem Argument offset angegeben:

# Poisson-Modell mit der Fläche als Offset
resultat <- glm(events ~ väder + temperatur + vecka + offset(log(area)),
                family = poisson(link = "log"), data = data)

Beachten Sie, dass der Offset auf der logarithmischen Skala anzugeben ist (log(area)), da das Modell logarithmisch ist.

Annahmen der Poisson-Regression

• Die Zielgröße ist eine Anzahl (nichtnegative ganze Zahlen).
• Die Beobachtungen sind voneinander unabhängig.
• Der Logarithmus der erwarteten Anzahl ist eine lineare Funktion der Prädiktoren. Stetige Prädiktoren müssen gegebenenfalls transformiert oder flexibel modelliert werden (etwa mit Splines), wenn der Zusammenhang nicht log-linear ist.
• Die Varianz ist gleich dem Erwartungswert (Äquidispersion). Diese Annahme ist häufig verletzt, siehe nächster Abschnitt.
• Die Exposition (Zeit, Fläche, Bevölkerung) ist korrekt als Offset angegeben, wenn die Rate von Interesse ist.

Überdispersion und Nullinflation

Das Poisson-Modell setzt voraus, dass die Varianz gleich dem Erwartungswert ist. In realen Daten ist die Varianz oft größer als der Erwartungswert, was als Überdispersion bezeichnet wird. Ursache ist häufig, dass wichtige Prädiktoren im Modell fehlen, dass die Beobachtungen nicht unabhängig sind (zum Beispiel wiederholte Messungen an derselben Person) oder dass sich die Personen auf eine Weise unterscheiden, die das Modell nicht erfasst. Die Punktschätzer werden selten nennenswert beeinflusst, doch die Standardfehler werden zu klein, die Konfidenzintervalle zu schmal und die p-Werte zu niedrig.

Überdispersion lässt sich beurteilen, indem man die Pearson-Chi2-Statistik (oder die Residualdevianz) durch die Anzahl der Freiheitsgrade teilt. Ein Quotient nahe 1 spricht dafür, dass die Poisson-Annahme erfüllt ist. Ein Quotient deutlich über 1 (Faustregel: über 1,5) spricht für Überdispersion.

# Quotient aus Pearson-Chi2 und Anzahl der Freiheitsgrade
sum(residuals(model, type = "pearson")^2) / df.residual(model)

Bei Überdispersion gibt es mehrere Möglichkeiten:

# 1. Quasi-Poisson: gleiche Schätzer, die Standardfehler werden mit der Dispersion hochskaliert
model.qp <- glm(y ~ x1 + x2 + offset(log(tid)), family = quasipoisson, data = data)

# 2. Negative Binomialregression: modelliert zusätzliche Variation mit einem eigenen Parameter
library(MASS)
model.nb <- glm.nb(y ~ x1 + x2 + offset(log(tid)), data = data)

# 3. Robuste (Sandwich-)Standardfehler für ein gewöhnliches Poisson-Modell
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

Ein weiteres häufiges Problem ist ein Überschuss an Nullen (engl. zero inflation). Er entsteht, wenn die Nullen aus zwei verschiedenen Prozessen stammen. Stellen Sie sich vor, wir untersuchen die Anzahl der pro Tag gerauchten Zigaretten in einer Schulklasse, in der die Hälfte der Klasse raucht. Bei den Rauchern folgt die Anzahl der Zigaretten einem Poisson-Prozess, während die Nichtraucher stets null Zigaretten rauchen. Die Daten enthalten dann deutlich mehr Nullen, als die Poisson-Verteilung vorhersagt. In solchen Situationen eignet sich ein nullinflationiertes Modell besser, zum Beispiel zeroinfl() aus dem Package pscl, das die Wahrscheinlichkeit einer "strukturellen" Null und die Anzahl bei den Übrigen in zwei Teilen modelliert:

library(pscl)
# Vor | : das Modell für die Anzahl. Nach | : das Modell für die strukturellen Nullen.
model.zip <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "poisson")
model.zinb <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "negbin")

Beispiel: Hautkrebs in Minneapolis und Dallas

Die Daten stammen aus einer klassischen Studie zum nicht-melanozytären Hautkrebs bei weißen Einwohnern in vier Regionen der USA (Scotto J, Kopf AW, Urbach F. Non-melanoma skin cancer among Caucasians in four areas of the United States. Cancer. 1974;34(4):1333-8. PMID 4422100). Hier werden zwei der Städte verwendet, Minneapolis und Dallas. Die Daten bestehen aus der Anzahl der Neuerkrankungen (cases) und der Bevölkerungsgröße (n) in acht Altersgruppen je Stadt. Der Datensatz ist klein und kann mit read.table() direkt in R eingelesen werden. Die Variablen u1 bis u7 sind Indikatorvariablen für die Altersgruppen und werden in der folgenden Analyse nicht benötigt, da die Altersgruppe dort stattdessen als Faktor erzeugt wird.

# Liest die Daten direkt in R ein
nonmel <- read.table(header = TRUE,
text = "
cases city u1 u2 u3 u4 u5 u6 u7 n
1 1 0 1 0 0 0 0 0 0 172675
2 16 0 0 1 0 0 0 0 0 123065
3 30 0 0 0 1 0 0 0 0 96216
4 71 0 0 0 0 1 0 0 0 92051
5 102 0 0 0 0 0 1 0 0 72159
6 130 0 0 0 0 0 0 1 0 54722
7 133 0 0 0 0 0 0 0 1 32185
8 40 0 0 0 0 0 0 0 0 8328
9 4 1 1 0 0 0 0 0 0 181343
10 38 1 0 1 0 0 0 0 0 146207
11 119 1 0 0 1 0 0 0 0 121374
12 221 1 0 0 0 1 0 0 0 111353
13 259 1 0 0 0 0 1 0 0 83004
14 310 1 0 0 0 0 0 1 0 55932
15 226 1 0 0 0 0 0 0 1 29007
16 65 1 0 0 0 0 0 0 0 7583
")

# Erzeugt Variablen für Altersgruppe und Stadt.
# age.range ist ein Faktor mit 8 Stufen, wobei die höchste Altersgruppe (85+) die Referenzstufe ist.
# within() wird verwendet, damit man nicht vor jede Variable nonmel$ schreiben muss.
nonmel <- within(nonmel, {
  age.range <- rep(c("15_24", "25_34", "35_44", "45_54", "55_64", "65_74", "75_84", "85+"), 2)
  age.range <- factor(age.range)
  age.range <- relevel(age.range, ref = "85+")
  city <- factor(city, 0:1, c("Minneapolis", "Dallas"))
})

Das Modell

Uns interessiert die Inzidenz von Hautkrebs, also die Anzahl der Fälle im Verhältnis zur Bevölkerungsgröße, und wie sie sich zwischen den Städten und zwischen den Altersgruppen unterscheidet. Die Anzahl der Fälle ist die abhängige Variable, und die Bevölkerungsgröße geht als Offset (auf der logarithmischen Skala) ein:

# Der Offset kann in der Formel angegeben werden
model.1 <- glm(cases ~ city + age.range + offset(log(n)),
               family = poisson(link = "log"), data = nonmel)

# oder mit dem Argument offset (dasselbe Modell)
model.1 <- glm(cases ~ city + age.range, offset = log(n),
               family = poisson(link = "log"), data = nonmel)

summary(model.1)

Call:
glm(formula = cases ~ city + age.range, family = poisson(link = "log"),
    data = nonmel, offset = log(n))

Deviance Residuals:
     Min        1Q    Median        3Q       Max
-1.50598  -0.48566   0.01639   0.36926   1.24763

Coefficients:
               Estimate Std. Error z value Pr(>|z|)
(Intercept)     -5.4834     0.1037 -52.890  < 2e-16 ***
cityDallas       0.8039     0.0522  15.399  < 2e-16 ***
age.range15_24  -6.1742     0.4577 -13.488  < 2e-16 ***
age.range25_34  -3.5440     0.1675 -21.160  < 2e-16 ***
age.range35_44  -2.3268     0.1275 -18.254  < 2e-16 ***
age.range45_54  -1.5790     0.1138 -13.871  < 2e-16 ***
age.range55_64  -1.0869     0.1109  -9.800  < 2e-16 ***
age.range65_74  -0.5288     0.1086  -4.868 1.13e-06 ***
age.range75_84  -0.1157     0.1109  -1.042    0.297
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 2789.6810  on 15  degrees of freedom
Residual deviance:    8.2585  on  7  degrees of freedom
AIC: 120.5

Number of Fisher Scoring iterations: 4

Die Koeffizienten in summary() liegen auf der logarithmischen Skala. Der Achsenabschnitt (−5,4834) ist der Logarithmus der Inzidenz in der Referenzgruppe, also bei Personen ab 85 Jahren in Minneapolis: $e^{-5{,}4834} = 0{,}0042$ Fälle pro Person oder 4,2 Fälle pro 1 000 Personen im Studienzeitraum.

Inzidenzratenverhältnisse mit Konfidenzintervallen

Um die Inzidenzratenverhältnisse zu erhalten, exponentiert man die Koeffizienten und ihre Konfidenzintervalle. Mit dem Package broom und der Funktion tidy() wird die Ausgabe zu einer übersichtlichen Tabelle:

library(broom)
# exponentiate = TRUE liefert Inzidenzratenverhältnisse (IRR) statt Log-Schätzern
tidy(model.1, exponentiate = TRUE, conf.int = TRUE)

# Dasselbe in Basis-R, mit Wald-basierten Konfidenzintervallen
exp(cbind(IRR = coef(model.1), confint.default(model.1)))

Die Inzidenzratenverhältnisse mit 95-%-Konfidenzintervall (Wald) lauten:

Wie wird das Ergebnis interpretiert?

• Stadt. Das Inzidenzratenverhältnis für Dallas beträgt 2,23. Die Inzidenz von Hautkrebs ist in Dallas also bei gleicher Altersverteilung 2,23-mal so hoch wie in Minneapolis (123 % höher). Es handelt sich um einen relativen Unterschied der Rate, nicht um 2,23 zusätzliche Fälle pro Person. Minneapolis ist die Referenzkategorie und erscheint daher nicht in der Tabelle. Der Unterschied passt zur höheren UV-Strahlung in Dallas, das deutlich weiter südlich liegt.
• Alter. Die älteste Gruppe (85+) ist die Referenzkategorie. Alle jüngeren Altersgruppen haben Inzidenzratenverhältnisse unter 1, also eine niedrigere Inzidenz als die Ältesten, und das Verhältnis steigt mit dem Alter stufenweise an. Personen im Alter von 15–24 Jahren haben eine Inzidenz von etwa 0,2 % der Inzidenz bei Personen ab 85 Jahren. Für 75–84 Jahre schließt das Konfidenzintervall die 1 ein (p = 0,30): Die Inzidenz unterscheidet sich nicht sicher von derjenigen in der ältesten Gruppe.

Die Inzidenzratenverhältnisse lassen sich in einem Forest-Plot darstellen. Die Skala ist logarithmisch, sodass Verhältnisse unter und über 1 denselben Abstand zur Linie für keinen Unterschied haben. Mit dem Package sjPlot erstellt man ein solches Diagramm mit plot_model(model.1); das folgende Diagramm zeigt dieselben Schätzer.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Inzidenzratenverhältnisse (IRR) mit 95-%-Konfidenzintervall aus dem Poisson-Modell für Hautkrebs: Stadt und Altersgruppe.",
  "width": "container",
  "height": 300,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"term": "Dallas (im Vergleich zu Minneapolis)", "irr": 2.2342, "lo": 2.017, "hi": 2.4749},
      {"term": "75–84 Jahre", "irr": 0.8907, "lo": 0.7167, "hi": 1.107},
      {"term": "65–74 Jahre", "irr": 0.5893, "lo": 0.4763, "hi": 0.7291},
      {"term": "55–64 Jahre", "irr": 0.3373, "lo": 0.2714, "hi": 0.4191},
      {"term": "45–54 Jahre", "irr": 0.2062, "lo": 0.165, "hi": 0.2577},
      {"term": "35–44 Jahre", "irr": 0.0976, "lo": 0.076, "hi": 0.1253},
      {"term": "25–34 Jahre", "irr": 0.0289, "lo": 0.0208, "hi": 0.0401},
      {"term": "15–24 Jahre", "irr": 0.0021, "lo": 0.0008, "hi": 0.0051}
    ]
  },
  "encoding": {
    "y": {"field": "term", "type": "nominal", "sort": null, "title": null, "axis": {"labelFontSize": 13, "labelLimit": 260}},
    "x": {"type": "quantitative", "scale": {"type": "log", "domain": [0.0005, 5]}, "title": "Inzidenzratenverhältnis (IRR), logarithmische Skala", "axis": {"labelFontSize": 12, "titleFontSize": 13, "values": [0.001, 0.01, 0.1, 0.5, 1, 2, 5]}}
  },
  "layer": [
    {"mark": {"type": "rule", "strokeDash": [4, 4], "color": "#8a8a96"}, "encoding": {"x": {"datum": 1}, "y": null}},
    {"mark": {"type": "rule", "strokeWidth": 2, "color": "#0e7490"}, "encoding": {"x": {"field": "lo"}, "x2": {"field": "hi"}}},
    {"mark": {"type": "point", "filled": true, "size": 80, "color": "#0e7490"}, "encoding": {"x": {"field": "irr"}, "tooltip": [{"field": "term", "title": "Vergleich"}, {"field": "irr", "title": "IRR", "format": ".3f"}, {"field": "lo", "title": "Untere Grenze", "format": ".3f"}, {"field": "hi", "title": "Obere Grenze", "format": ".3f"}]}}
  ]
}

Passt das Modell?

Nach einer Poisson-Regression sollten stets zwei Prüfungen durchgeführt werden.

Anpassungsgüte (goodness of fit). Die Residualdevianz (8,26 bei 7 Freiheitsgraden) kann mit einer Chi2-Verteilung verglichen werden. Ein hoher p-Wert bedeutet, dass es keine Hinweise auf eine schlechte Anpassung des Modells gibt:

pchisq(deviance(model.1), df.residual(model.1), lower.tail = FALSE)
#> [1] 0.3103526

Überdispersion. Der Quotient aus Residualdevianz und Freiheitsgraden beträgt 8,26 / 7 = 1,18, liegt also nahe 1. Die Poisson-Annahme, dass Erwartungswert und Varianz gleich sind, erscheint plausibel, und in diesem Beispiel besteht kein Bedarf an Quasi-Poisson- oder negativer Binomialregression.

Vorhersagen

Mit predict() kann das Modell die Inzidenz für beliebige Kombinationen aus Stadt und Altersgruppe vorhersagen. Da das Modell logarithmisch ist, exponentiert man die Vorhersage (oder gibt type = "response" an). Setzt man die Bevölkerung n auf 1, erhält man die erwartete Anzahl der Fälle pro Person, mit n = 1 000 die Anzahl pro 1 000 Personen. Eine Anzahl pro 1 000 Personen ist leichter zu interpretieren als eine kleine Dezimalzahl.

# Erwartete Anzahl der Fälle pro Person, Minneapolis, 85 Jahre und älter
predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1),
        type = "response")
#>           1
#> 0.004155093

# Pro 1 000 Personen
predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1000),
        type = "response")
#>        1
#> 4.155093

In der ältesten Altersgruppe in Minneapolis sind im Studienzeitraum also 4,2 Fälle pro 1 000 Personen zu erwarten. Der entsprechende Wert für Dallas beträgt 4,2 × 2,23 = 9,3 pro 1 000.

Vorhersagen für alle Zeilen des Datensatzes werden auf dieselbe Weise erstellt:

# Neue Datensätze mit derselben Stadt und Altersgruppe wie im Original
newdat1 <- nonmel
newdat2 <- nonmel

# Die Bevölkerung wird auf 1 gesetzt, um die Anzahl der Fälle pro Person zu erhalten
newdat1$n <- 1
nonmel$pred.cases.per.one <- predict(model.1, newdat1, type = "response")

# Die Bevölkerung wird auf 1 000 gesetzt, um die Anzahl der Fälle pro 1 000 Personen zu erhalten
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- predict(model.1, newdat2, type = "response")

# Erwartete Anzahl der Fälle in der tatsächlichen Bevölkerung
nonmel$pred.cases <- predict(model.1, type = "response")

nonmel[, c("city", "age.range", "cases", "n", "pred.cases.per.thousand", "pred.cases")]

          city age.range cases      n pred.cases.per.thousand pred.cases
1  Minneapolis     15_24     1 172675              0.00865295   1.494148
2  Minneapolis     25_34    16 123065              0.12007322  14.776810
3  Minneapolis     35_44    30  96216              0.40555928  39.021292
4  Minneapolis     45_54    71  92051              0.85670789  78.860818
5  Minneapolis     55_64   102  72159              1.40134673 101.119779
6  Minneapolis     65_74   130  54722              2.44870028 133.997777
7  Minneapolis     75_84   133  32185              3.70128210 119.125764
8  Minneapolis       85+    40   8328              4.15509270  34.603612
9       Dallas     15_24     4 181343              0.01933271   3.505852
10      Dallas     25_34    38 146207              0.26827163  39.223190
11      Dallas     35_44   119 121374              0.90611423 109.978708
12      Dallas     45_54   221 111353              1.91408567 213.139182
13      Dallas     55_64   259  83004              3.13093612 259.880221
14      Dallas     65_74   310  55932              5.47096874 306.002223
15      Dallas     75_84   226  29007              8.26952928 239.874236
16      Dallas       85+    65   7583              9.28344824  70.396388

Die Spalte pred.cases lässt sich direkt mit der beobachteten Anzahl der Fälle (cases) vergleichen. Die Übereinstimmung ist in allen Gruppen gut, was zur niedrigen Residualdevianz passt.

Zusammenfassung

• Die Poisson-Regression wird verwendet, wenn die Zielgröße eine Anzahl von Ereignissen ist, und mit einem Offset, wenn die Zielgröße eine Rate ist (Ereignisse pro Zeit, Fläche oder Bevölkerung).
• Das Modell ist ein GLM mit logarithmischer Linkfunktion. Exponentierte Koeffizienten sind Inzidenzratenverhältnisse (IRR), die multiplikativ interpretiert werden.
• Die Poisson-Verteilung setzt voraus, dass die Varianz gleich dem Erwartungswert ist. Prüfen Sie stets die Anpassung und die Dispersion des Modells.
• Bei Überdispersion verwendet man Quasi-Poisson, negative Binomialregression oder robuste Standardfehler. Bei einem Überschuss an Nullen verwendet man nullinflationierte Modelle.