Poissonregresjon

Introduksjon til Poissonregresjon

Poissonregresjon brukes til å modellere antall ganger en hendelse inntreffer i løpet av en gitt tid eller innenfor et gitt område. Utfallet er altså et antall (eng. count): 0, 1, 2, 3 og så videre. For å bruke metoden riktig må man forstå noen sentrale begreper og den sannsynlighetsfordelingen metoden bygger på, nemlig Poissonfordelingen. Kapitlet gjennomgår begrepene, Poissonfordelingen og modellen, og avsluttes med et fullstendig eksempel i R.

Anta at vi undersøker hvor mange pasienter som oppsøker et legesenter per dag. Vi vil vite hvilke faktorer (prediktorer) som påvirker antall besøk, og vi vil kunne forutsi (predikere) hvor mange besøk som kan forventes en bestemt dag. Hvert besøk er en hendelse (eng. event), og antall besøk per dag er et antall. Settes antallet i forhold til tiden det gjelder, får man en rate (hendelser per tidsenhet): hvis 10 pasienter kommer i løpet av en dag, er raten 10 pasienter per dag. Raten kan også gjelde et område, for eksempel antall som søker helsehjelp per kvadratkilometer.

Poissonregresjon kan brukes så lenge hendelsene oppstår gjennom en Poissonprosess, det vil si en prosess som gir opphav til en Poissonfordeling (se nedenfor). Hva hendelsen er, spiller ingen rolle. Metoden brukes til å studere antall dødsfall, antall sykehusinnleggelser, antall infeksjoner, antall kunder eller antall e-postmeldinger.

Poissonfordelingen

Poissonfordelingen ble beskrevet av den franske matematikeren Siméon Denis Poisson. Den angir sannsynligheten for at et bestemt antall hendelser skal inntreffe innenfor et gitt intervall i tid eller rom, forutsatt at hendelsene inntreffer uavhengig av hverandre og med en konstant gjennomsnittlig takt. Fordelingen har bare én parameter, λ (lambda), som er det forventede antallet hendelser i intervallet. Sannsynligheten for å observere nøyaktig $k$ hendelser er

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

En karakteristisk egenskap er at gjennomsnittet og variansen begge er lik λ. Det er nettopp denne egenskapen som ofte ikke holder i virkelige data (se avsnittet om overdispersjon).

Figuren viser Poissonfordelingen for tre verdier av λ. Fordelingen er diskret: sannsynlighet finnes bare for heltall. Når λ er liten, ligger det meste av sannsynligheten nær null, og null hendelser er vanlig. Når λ øker, forskyves fordelingen mot høyre, blir bredere og stadig mer symmetrisk, og null hendelser blir usannsynlig. Fordelingens topp ligger ved λ eller like under (for heltallige λ er λ og λ − 1 like sannsynlige), mens gjennomsnittet alltid er nøyaktig λ.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Poissonfordelingen for λ = 1, 4 og 10: sannsynligheten for hvert antall hendelser 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": "Antall hendelser (k)", "scale": {"domain": [0, 20]}, "axis": {"tickMinStep": 1, "labelFontSize": 12, "titleFontSize": 13, "grid": false}},
    "y": {"field": "p", "type": "quantitative", "title": "Sannsynlighet 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": "Forventet antall"}, {"field": "k", "title": "Antall hendelser"}, {"field": "p", "title": "Sannsynlighet", "format": ".4f"}]}}
  ]
}

Eksempel. Du arbeider på et legesenter der i gjennomsnitt 10 pasienter kommer per dag, så λ = 10. Antallet varierer fra dag til dag. Hvis pasientene kommer uavhengig av hverandre, kan Poissonfordelingen brukes til å beregne hvor sannsynlig det er at nøyaktig 5 pasienter kommer en bestemt dag, eller at flere enn 5 kommer. I R gjøres det med dpois() (sannsynligheten for nøyaktig k) og ppois() (sannsynligheten for k eller færre):

# Sannsynligheten for at nøyaktig 5 pasienter kommer når forventet antall er 10
dpois(5, lambda = 10)
#> [1] 0.03783327

# Sannsynligheten for at flere enn 5 pasienter kommer
ppois(5, lambda = 10, lower.tail = FALSE)
#> [1] 0.932914

Nøyaktig 5 pasienter inntreffer altså omtrent hver 26. dag (3,8 %), mens flere enn 5 pasienter kommer på drøyt 93 % av dagene.

Forutsetninger for Poissonfordelingen

• Utfallet $k$ er et antall hendelser innenfor et bestemt intervall og er et ikke-negativt heltall (0, 1, 2, 3 og så videre).
• Hendelsene inntreffer uavhengig av hverandre: at én hendelse inntreffer, påvirker ikke når neste hendelse inntreffer.
• Den gjennomsnittlige takten som hendelsene inntreffer med, er konstant innenfor intervallet.
• To hendelser kan ikke inntreffe nøyaktig samtidig.

Poissonregresjon er en generalisert lineær modell

Poissonregresjon er, som navnet antyder, en type regresjonsanalyse. Den brukes når den avhengige variabelen (utfallet, Y) er et antall. Poissonregresjon hører til de generaliserte lineære modellene (eng. generalized linear models, GLM).

Hva er en GLM?

Generaliserte lineære modeller er et rammeverk som omfatter flere statistiske modeller, blant annet vanlig lineær regresjon. I lineær regresjon modelleres Y direkte som en lineær funksjon av prediktorene ($X_1, X_2, X_3, \ldots$) pluss et normalfordelt feilledd. Det fungerer dårlig for antall: et antall kan ikke være negativt, og spredningen øker ofte med gjennomsnittet. I en GLM tillates utfallet i stedet å følge andre fordelinger fra den såkalte eksponentialfamilien (blant annet Poisson-, binomial-, gamma- og normalfordelingen), og sammenhengen mellom prediktorene og utfallets forventede verdi går via en linkfunksjon.

Poissonregresjon bruker den naturlige logaritmen som linkfunksjon. Det er logaritmen av det forventede antallet, $\mu = E(Y)$, som modelleres som en lineær funksjon av prediktorene:

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

I motsetning til lineær regresjon finnes det ikke noe separat feilledd i ligningen: den tilfeldige variasjonen rundt $\mu_i$ beskrives i stedet av Poissonfordelingen. Logaritmen gjør at det forventede antallet alltid blir positivt, $\mu_i = e^{\beta_0 + \beta_1 X_{1i} + \cdots}$, og at prediktorenes effekter blir multiplikative.

Tolkning av koeffisientene: insidensrateratio (IRR)

Fordi modellen er logaritmisk, tolkes koeffisientene etter eksponensiering. $e^{\beta}$ kalles insidensrateratio (eng. incidence rate ratio, IRR) og angir med hvilken faktor det forventede antallet (eller raten) endres når prediktoren øker med én enhet, eller for en kategori sammenlignet med referansekategorien, med de øvrige prediktorene uendret. IRR = 1 betyr ingen forskjell, IRR = 2 betyr en fordoblet rate og IRR = 0,5 en halvert rate.

Eksempler på spørsmål som besvares med Poissonregresjon

Eksempel 1. Hvordan påvirker utetemperaturen og været antall forkjølelser som diagnostiseres på et legesenter per uke? Dataene kan se slik ut (de første 13 ukene vises):

Uke, temperatur og vær er prediktorer som brukes til å forklare og predikere antall forkjølelser. Modellen skrives

$$\log(\text{forventet antall forkjølelser}_i) = \beta_0 + \beta_1,\text{uke}_i + \beta_2,\text{temperatur}_i + \beta_3,\text{vær}_i$$

der $\beta_0$ er konstantleddet (logaritmen av det forventede antallet når alle prediktorer er null eller i sin referansekategori) og $\beta_1, \beta_2, \beta_3$ er prediktorenes koeffisienter. Vær er en kategorisk variabel og får i praksis én koeffisient per kategori unntatt referansekategorien.

Eksempel 2. Hvor mange kunder gjennomfører et kjøp i en butikk per dag, og hva påvirker det?

Eksempel 3. Hvor mange personer reiser med en bestemt busslinje til jobb, og hvordan påvirkes antallet av vær og ukedag?

Offset: å analysere rater

Ofte er det ikke antallet i seg selv som er interessant, men antallet i forhold til en eksponering: tid, areal, distanse eller befolkningsstørrelse. To eksempler:

• Du studerer dødsfall og vet hvor lenge hvert individ har vært fulgt opp (eksponeringstiden). Da er det insidensen, antall hendelser delt på eksponeringstiden, som er interessant. Begrepet event rate brukes ofte om antall hendelser per tidsenhet.
• Du studerer antall dyrearter i ulike geografiske områder av ulik størrelse. Da er det antall arter per areal som er interessant.

Eksponeringen føres inn i modellen som en offset: logaritmen av eksponeringen legges til på høyre side med en koeffisient som er fiksert til 1,

$$\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$$

Modellen beskriver da raten $\mu_i / t_i$, og $e^{\beta}$ blir en ratio mellom rater. I R angis offseten med funksjonen offset() i formelen, eller med argumentet offset:

# Poissonmodell der areal brukes som offset
resultat <- glm(events ~ väder + temperatur + vecka + offset(log(area)),
                family = poisson(link = "log"), data = data)

Merk at offseten skal angis på logaritmisk skala (log(area)), fordi modellen er logaritmisk.

Forutsetninger for Poissonregresjon

• Utfallet er et antall (ikke-negative heltall).
• Observasjonene er uavhengige av hverandre.
• Logaritmen av det forventede antallet er en lineær funksjon av prediktorene. Kontinuerlige prediktorer kan måtte transformeres eller modelleres fleksibelt (for eksempel med splines) hvis sammenhengen ikke er log-lineær.
• Variansen er lik gjennomsnittet (ekvidispersjon). Denne forutsetningen brytes ofte, se neste avsnitt.
• Eksponeringen (tid, areal, befolkning) er korrekt angitt som offset når raten er det interessante.

Overdispersjon og nullinflasjon

Poissonmodellen forutsetter at variansen er lik gjennomsnittet. I virkelige data er variansen ofte større enn gjennomsnittet, noe som kalles overdispersjon. Det skyldes ofte at viktige prediktorer mangler i modellen, at observasjonene ikke er uavhengige (for eksempel gjentatte målinger på samme individ) eller at individene skiller seg fra hverandre på måter modellen ikke fanger opp. Punktestimatene påvirkes sjelden nevneverdig, men standardfeilene blir for små, konfidensintervallene for smale og p-verdiene for lave.

Overdispersjon kan vurderes ved å dele Pearsons kjikvadratstatistikk (eller residualdeviansen) på antall frihetsgrader. En ratio nær 1 taler for at Poissonforutsetningen holder. En ratio klart over 1 (tommelfingerregel: over 1,5) taler for overdispersjon.

# Ratioen mellom Pearsons kjikvadrat og antall frihetsgrader
sum(residuals(model, type = "pearson")^2) / df.residual(model)

Ved overdispersjon finnes det flere alternativer:

# 1. Kvasi-Poisson: samme estimater, standardfeilene skaleres opp med dispersjonen
model.qp <- glm(y ~ x1 + x2 + offset(log(tid)), family = quasipoisson, data = data)

# 2. Negativ binomial regresjon: modellerer ekstra variasjon med en egen parameter
library(MASS)
model.nb <- glm.nb(y ~ x1 + x2 + offset(log(tid)), data = data)

# 3. Robuste (sandwich-) standardfeil for en vanlig Poissonmodell
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

Et annet vanlig problem er overskudd av nuller (eng. zero inflation). Det oppstår når nullene kommer fra to ulike prosesser. Tenk deg at vi studerer antall sigaretter som røykes per dag i en skoleklasse der halve klassen røyker. Blant røykerne følger antall sigaretter en Poissonprosess, mens ikke-røykerne alltid røyker null sigaretter. Dataene inneholder da betydelig flere nuller enn Poissonfordelingen forutsier. I slike situasjoner passer en nullinflatert modell bedre, for eksempel zeroinfl() i pakken pscl, som modellerer sannsynligheten for en "strukturell" null og antallet blant de øvrige i to deler:

library(pscl)
# Før | : modellen for antallet. Etter | : modellen for de strukturelle nullene.
model.zip <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "poisson")
model.zinb <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "negbin")

Eksempel: hudkreft i Minneapolis og Dallas

Dataene kommer fra en klassisk studie av ikke-melanom hudkreft hos hvite innbyggere i fire områder i 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). Her brukes to av byene, Minneapolis og Dallas. Dataene består av antall nye tilfeller (cases) og befolkningens størrelse (n) i åtte aldersgrupper i hver by. Datasettet er lite og kan leses inn direkte i R med read.table(). Variablene u1 til u7 er indikatorvariabler for aldersgruppene og trengs ikke i analysen nedenfor, der aldersgruppen i stedet lages som en faktor.

# Leser inn data direkte i R
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
")

# Lager variabler for aldersgruppe og by.
# age.range er en faktor med 8 nivåer der den høyeste aldersgruppen (85+) er referansenivå.
# within() brukes for å slippe å skrive nonmel$ foran hver variabel.
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"))
})

Modellen

Vi er interessert i insidensen av hudkreft, det vil si antall tilfeller i forhold til befolkningens størrelse, og hvordan den varierer mellom byene og mellom aldersgruppene. Antall tilfeller er avhengig variabel, og befolkningens størrelse føres inn som offset (på logaritmisk skala):

# Offset kan angis i formelen
model.1 <- glm(cases ~ city + age.range + offset(log(n)),
               family = poisson(link = "log"), data = nonmel)

# eller med argumentet offset (samme 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

Koeffisientene i summary() er på logaritmisk skala. Konstantleddet (−5,4834) er logaritmen av insidensen i referansegruppen, altså personer 85 år og eldre i Minneapolis: $e^{-5{,}4834} = 0{,}0042$ tilfeller per person, eller 4,2 tilfeller per 1 000 personer i løpet av studieperioden.

Insidensrateratio med konfidensintervall

For å få insidensrateratioene eksponensierer man koeffisientene og konfidensintervallene deres. Med pakken broom og funksjonen tidy() blir utskriften en ryddig tabell:

library(broom)
# exponentiate = TRUE gir insidensrateratio (IRR) i stedet for log-estimater
tidy(model.1, exponentiate = TRUE, conf.int = TRUE)

# Det samme i base-R, med Wald-baserte konfidensintervaller
exp(cbind(IRR = coef(model.1), confint.default(model.1)))

Insidensrateratioene med 95 % konfidensintervall (Wald) blir:

Hvordan tolkes resultatet?

• By. Insidensrateratioen for Dallas er 2,23. Insidensen av hudkreft er altså 2,23 ganger høyere i Dallas enn i Minneapolis (123 % høyere) ved samme aldersfordeling. Det er en relativ forskjell i rate, ikke 2,23 flere tilfeller per person. Minneapolis er referansekategori og vises derfor ikke i tabellen. Forskjellen er forenlig med den høyere UV-strålingen i Dallas, som ligger betydelig lenger sør.
• Alder. Den eldste gruppen (85+) er referansekategori. Alle yngre aldersgrupper har insidensrateratio under 1, det vil si lavere insidens enn de eldste, og ratioen stiger trinnvis med alderen. Personer 15–24 år har en insidens som er omtrent 0,2 % av insidensen hos personer 85 år og eldre. For 75–84 år inkluderer konfidensintervallet 1 (p = 0,30): insidensen skiller seg ikke sikkert fra insidensen i den eldste gruppen.

Insidensrateratioene kan vises i et forestplot. Skalaen er logaritmisk, slik at ratioer under og over 1 får samme avstand fra linjen for ingen forskjell. Med pakken sjPlot lages et slikt diagram med plot_model(model.1); diagrammet nedenfor viser de samme estimatene.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Insidensrateratio (IRR) med 95 % konfidensintervall fra Poissonmodellen for hudkreft: by og aldersgruppe.",
  "width": "container",
  "height": 300,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"term": "Dallas (sammenlignet med Minneapolis)", "irr": 2.2342, "lo": 2.017, "hi": 2.4749},
      {"term": "75–84 år", "irr": 0.8907, "lo": 0.7167, "hi": 1.107},
      {"term": "65–74 år", "irr": 0.5893, "lo": 0.4763, "hi": 0.7291},
      {"term": "55–64 år", "irr": 0.3373, "lo": 0.2714, "hi": 0.4191},
      {"term": "45–54 år", "irr": 0.2062, "lo": 0.165, "hi": 0.2577},
      {"term": "35–44 år", "irr": 0.0976, "lo": 0.076, "hi": 0.1253},
      {"term": "25–34 år", "irr": 0.0289, "lo": 0.0208, "hi": 0.0401},
      {"term": "15–24 år", "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": "Insidensrateratio (IRR), logaritmisk 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": "Sammenligning"}, {"field": "irr", "title": "IRR", "format": ".3f"}, {"field": "lo", "title": "Nedre grense", "format": ".3f"}, {"field": "hi", "title": "Øvre grense", "format": ".3f"}]}}
  ]
}

Passer modellen?

To kontroller bør alltid gjøres etter en Poissonregresjon.

Tilpasning (goodness of fit). Residualdeviansen (8,26 på 7 frihetsgrader) kan sammenlignes med en kjikvadratfordeling. En høy p-verdi betyr at det ikke finnes tegn på at modellen passer dårlig:

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

Overdispersjon. Ratioen mellom residualdeviansen og frihetsgradene er 8,26 / 7 = 1,18, altså nær 1. Poissonforutsetningen om likt gjennomsnitt og lik varians virker rimelig, og det er ikke behov for kvasi-Poisson eller negativ binomial regresjon i dette eksemplet.

Prediksjoner

Med predict() kan modellen forutsi insidensen for valgfrie kombinasjoner av by og aldersgruppe. Fordi modellen er logaritmisk, eksponensierer man prediksjonen (eller angir type = "response"). Ved å sette befolkningen n til 1 får man forventet antall tilfeller per person, og med n = 1 000 antallet per 1 000 personer. Et antall per 1 000 personer er lettere å tolke enn et lite desimaltall.

# Forventet antall tilfeller per person, Minneapolis, 85 år og eldre
predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1),
        type = "response")
#>           1
#> 0.004155093

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

I den eldste aldersgruppen i Minneapolis forventes altså 4,2 tilfeller per 1 000 personer i løpet av studieperioden. Tilsvarende tall i Dallas er 4,2 × 2,23 = 9,3 per 1 000.

Prediksjoner for alle rader i datasettet lages på samme måte:

# Nye datasett med samme by og aldersgruppe som originalen
newdat1 <- nonmel
newdat2 <- nonmel

# Befolkningen settes til 1 for å få antall tilfeller per person
newdat1$n <- 1
nonmel$pred.cases.per.one <- predict(model.1, newdat1, type = "response")

# Befolkningen settes til 1 000 for å få antall tilfeller per 1 000 personer
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- predict(model.1, newdat2, type = "response")

# Forventet antall tilfeller i den faktiske befolkningen
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

Kolonnen pred.cases kan sammenlignes direkte med det observerte antallet tilfeller (cases). Samsvaret er godt i alle grupper, noe som stemmer med den lave residualdeviansen.

Sammendrag

• Poissonregresjon brukes når utfallet er et antall hendelser, og med en offset når utfallet er en rate (hendelser per tid, areal eller befolkning).
• Modellen er en GLM med logaritmisk linkfunksjon. Eksponensierte koeffisienter er insidensrateratioer (IRR) som tolkes multiplikativt.
• Poissonfordelingen forutsetter at variansen er lik gjennomsnittet. Kontroller alltid modellens tilpasning og dispersjon.
• Ved overdispersjon brukes kvasi-Poisson, negativ binomial regresjon eller robuste standardfeil. Ved overskudd av nuller brukes nullinflaterte modeller.