Poissonregression

Introduktion til Poissonregression

Poissonregression bruges til at modellere antallet af gange, en hændelse indtræffer inden for en given tid eller et givet område. Udfaldet er altså et antal (eng. count): 0, 1, 2, 3 og så videre. For at bruge metoden korrekt skal man forstå nogle centrale begreber og den sandsynlighedsfordeling, som metoden bygger på, nemlig Poissonfordelingen. Kapitlet gennemgår begreberne, Poissonfordelingen og modellen og afsluttes med et komplet eksempel i R.

Antag, at vi undersøger, hvor mange patienter der henvender sig i et lægehus per dag. Vi vil vide, hvilke faktorer (prædiktorer) der påvirker antallet af besøg, og vi vil kunne forudsige (prædiktere), hvor mange besøg der kan forventes en bestemt dag. Hvert besøg er en hændelse (eng. event), og antallet af besøg per dag er et antal. Sættes antallet i relation til den tid, det gælder, får man en rate (hændelser per tidsenhed): hvis 10 patienter henvender sig i løbet af en dag, er raten 10 patienter per dag. Raten kan også gælde et område, for eksempel antal personer, der søger behandling, per kvadratkilometer.

Poissonregression kan bruges, så længe hændelserne opstår gennem en Poissonproces, det vil sige en proces, der giver anledning til en Poissonfordeling (se nedenfor). Hvad hændelsen er, spiller ingen rolle. Metoden bruges til at studere antal dødsfald, antal hospitalsindlæggelser, antal infektioner, antal kunder eller antal e-mails.

Poissonfordelingen

Poissonfordelingen blev beskrevet af den franske matematiker Siméon Denis Poisson. Den angiver sandsynligheden for, at et bestemt antal hændelser indtræffer inden for et givet interval i tid eller rum, forudsat at hændelserne indtræffer uafhængigt af hinanden og med en konstant gennemsnitlig hyppighed. Fordelingen har kun én parameter, λ (lambda), som er det forventede antal hændelser i intervallet. Sandsynligheden for at observere præcis $k$ hændelser er

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

En karakteristisk egenskab er, at middelværdien og variansen begge er lig med λ. Det er netop denne egenskab, der ofte ikke holder i virkelige data (se afsnittet om overdispersion).

Figuren viser Poissonfordelingen for tre værdier af λ. Fordelingen er diskret: der er kun sandsynlighed for heltal. Når λ er lille, ligger det meste af sandsynligheden tæt på nul, og nul hændelser er almindeligt. Når λ stiger, forskydes fordelingen mod højre, bliver bredere og mere og mere symmetrisk, og nul hændelser bliver usandsynligt. Fordelingens top ligger ved λ eller lige under (for heltallige λ er λ og λ − 1 lige sandsynlige), mens middelværdien altid er præcis λ.

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

Eksempel. Du arbejder i et lægehus, hvor der i gennemsnit henvender sig 10 patienter per dag, så λ = 10. Antallet varierer fra dag til dag. Hvis patienterne kommer uafhængigt af hinanden, kan Poissonfordelingen bruges til at beregne, hvor sandsynligt det er, at præcis 5 patienter kommer en bestemt dag, eller at flere end 5 kommer. I R gøres det med dpois() (sandsynligheden for præcis k) og ppois() (sandsynligheden for k eller færre):

# Sandsynligheden for, at præcis 5 patienter kommer, når det forventede antal er 10
dpois(5, lambda = 10)
#> [1] 0.03783327

# Sandsynligheden for, at flere end 5 patienter kommer
ppois(5, lambda = 10, lower.tail = FALSE)
#> [1] 0.932914

Præcis 5 patienter forekommer altså cirka hver 26. dag (3,8 %), mens flere end 5 patienter kommer på godt 93 % af dagene.

Antagelser for Poissonfordelingen

• Udfaldet $k$ er et antal hændelser inden for et bestemt interval og er et ikke-negativt heltal (0, 1, 2, 3 og så videre).
• Hændelserne indtræffer uafhængigt af hinanden: at en hændelse indtræffer, påvirker ikke, hvornår den næste hændelse indtræffer.
• Den gennemsnitlige hyppighed, hvormed hændelser indtræffer, er konstant inden for intervallet.
• To hændelser kan ikke indtræffe præcis samtidig.

Poissonregression er en generaliseret lineær model

Poissonregression er, som navnet antyder, en type regressionsanalyse. Den bruges, når den afhængige variabel (udfaldet, Y) er et antal. Poissonregression hører til de generaliserede lineære modeller (eng. generalized linear models, GLM).

Hvad er en GLM?

Generaliserede lineære modeller er en ramme, der omfatter flere statistiske modeller, blandt andet almindelig lineær regression. I lineær regression modelleres Y direkte som en lineær funktion af prædiktorerne ($X_1, X_2, X_3, \ldots$) plus et normalfordelt fejlled. Det fungerer dårligt for antal: et antal kan ikke være negativt, og spredningen stiger ofte med middelværdien. I en GLM tillades udfaldet i stedet at følge andre fordelinger fra den såkaldte eksponentielle familie (blandt andet Poisson-, binomial-, gamma- og normalfordelingen), og sammenhængen mellem prædiktorerne og udfaldets forventede værdi går via en linkfunktion.

Poissonregression bruger den naturlige logaritme som linkfunktion. Det er logaritmen af det forventede antal, $\mu = E(Y)$, der modelleres som en lineær funktion af prædiktorerne:

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

I modsætning til lineær regression er der intet separat fejlled i ligningen: den tilfældige variation omkring $\mu_i$ beskrives i stedet af Poissonfordelingen. Logaritmen gør, at det forventede antal altid bliver positivt, $\mu_i = e^{\beta_0 + \beta_1 X_{1i} + \cdots}$, og at prædiktorernes effekter bliver multiplikative.

Fortolkning af koefficienterne: incidensrate-ratioer (IRR)

Da modellen er logaritmisk, fortolkes koefficienterne efter eksponentiering. $e^{\beta}$ kaldes incidensrate-ratio (eng. incidence rate ratio, IRR) og angiver, med hvilken faktor det forventede antal (eller raten) ændres, når prædiktoren øges med én enhed, eller for en kategori sammenlignet med referencekategorien, med de øvrige prædiktorer uændrede. IRR = 1 betyder ingen forskel, IRR = 2 betyder en fordoblet rate, og IRR = 0,5 en halveret rate.

Eksempler på spørgsmål, der besvares med Poissonregression

Eksempel 1. Hvordan påvirker udetemperaturen og vejret antallet af forkølelser, der diagnosticeres i et lægehus per uge? Data kan se sådan ud (de første 13 uger vises):

Uge, temperatur og vejr er prædiktorer, der bruges til at forklare og prædiktere antallet af forkølelser. Modellen skrives

$$\log(\text{forventet antal forkølelser}_i) = \beta_0 + \beta_1,\text{uge}_i + \beta_2,\text{temperatur}_i + \beta_3,\text{vejr}_i$$

hvor $\beta_0$ er interceptet (logaritmen af det forventede antal, når alle prædiktorer er nul eller i deres referencekategori), og $\beta_1, \beta_2, \beta_3$ er prædiktorernes koefficienter. Vejr er en kategorisk variabel og får i praksis én koefficient per kategori undtagen referencekategorien.

Eksempel 2. Hvor mange kunder gennemfører et køb i en butik per dag, og hvad påvirker det?

Eksempel 3. Hvor mange personer rejser med en bestemt buslinje til arbejde, og hvordan påvirkes antallet af vejr og ugedag?

Offset: at analysere rater

Ofte er det ikke antallet i sig selv, der er interessant, men antallet i relation til en eksponering: tid, areal, afstand eller befolkningsstørrelse. To eksempler:

• Du studerer dødsfald og ved, hvor længe hvert individ er fulgt op (eksponeringstiden). Så er det incidensen, antal hændelser divideret med eksponeringstiden, der er interessant. Termen event rate bruges ofte om antal hændelser per tidsenhed.
• Du studerer antallet af dyrearter i forskellige geografiske områder af forskellig størrelse. Så er det antallet af arter per areal, der er interessant.

Eksponeringen føres ind i modellen som en offset: logaritmen af eksponeringen lægges til på højresiden med en koefficient, der er fikseret 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}$ bliver et forhold mellem rater. I R angives offsetten med funktionen offset() i formlen eller med argumentet offset:

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

Bemærk, at offsetten skal angives på logaritmisk skala (log(area)), da modellen er logaritmisk.

Antagelser for Poissonregression

• Udfaldet er et antal (ikke-negative heltal).
• Observationerne er uafhængige af hinanden.
• Logaritmen af det forventede antal er en lineær funktion af prædiktorerne. Kontinuerte prædiktorer kan have behov for at blive transformeret eller modelleret fleksibelt (for eksempel med splines), hvis sammenhængen ikke er log-lineær.
• Variansen er lig med middelværdien (equidispersion). Denne antagelse brydes ofte, se næste afsnit.
• Eksponeringen (tid, areal, befolkning) er korrekt angivet som offset, når raten er det interessante.

Overdispersion og nulinflation

Poissonmodellen forudsætter, at variansen er lig med middelværdien. I virkelige data er variansen ofte større end middelværdien, hvilket kaldes overdispersion. Det skyldes ofte, at vigtige prædiktorer mangler i modellen, at observationerne ikke er uafhængige (for eksempel gentagne målinger på samme individ), eller at individerne adskiller sig fra hinanden på måder, som modellen ikke indfanger. Punktestimaterne påvirkes sjældent nævneværdigt, men standardfejlene bliver for små, konfidensintervallerne for smalle og p-værdierne for lave.

Overdispersion kan vurderes ved at dividere Pearsons chi2-statistik (eller residualdeviansen) med antallet af frihedsgrader. En ratio tæt på 1 taler for, at Poissonantagelsen holder. En ratio klart over 1 (tommelfingerregel: over 1,5) taler for overdispersion.

# Forholdet mellem Pearsons chi2 og antallet af frihedsgrader
sum(residuals(model, type = "pearson")^2) / df.residual(model)

Ved overdispersion er der flere muligheder:

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

# 2. Negativ binomialregression: modellerer ekstra variation med en særskilt parameter
library(MASS)
model.nb <- glm.nb(y ~ x1 + x2 + offset(log(tid)), data = data)

# 3. Robuste (sandwich-)standardfejl for en almindelig Poissonmodel
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

Et andet almindeligt problem er overskud af nuller (eng. zero inflation). Det opstår, når nullerne kommer fra to forskellige processer. Forestil dig, at vi studerer antallet af cigaretter, der ryges per dag i en skoleklasse, hvor halvdelen af klassen ryger. Blandt rygerne følger antallet af cigaretter en Poissonproces, mens ikke-rygerne altid ryger nul cigaretter. Data indeholder så betydeligt flere nuller, end Poissonfordelingen forudsiger. I sådanne situationer passer en nulinflateret model bedre, for eksempel zeroinfl() i pakken pscl, som modellerer sandsynligheden for et "strukturelt" nul og antallet blandt de øvrige i to dele:

library(pscl)
# Før | : modellen for antallet. Efter | : modellen for de strukturelle nuller.
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: hudkræft i Minneapolis og Dallas

Data stammer fra et klassisk studie af non-melanom hudkræft hos hvide indbyggere 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 bruges to af byerne, Minneapolis og Dallas. Data består af antallet af nye tilfælde (cases) og befolkningens størrelse (n) i otte aldersgrupper i hver by. Datasættet er lille og kan indlæses direkte i R med read.table(). Variablerne u1 til u7 er indikatorvariabler for aldersgrupperne og er ikke nødvendige i analysen nedenfor, hvor aldersgruppen i stedet oprettes som en faktor.

# Indlæser 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
")

# Opretter variabler for aldersgruppe og by.
# age.range er en faktor med 8 niveauer, hvor den højeste aldersgruppe (85+) er referenceniveau.
# within() bruges for at undgå at 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 interesserede i incidensen af hudkræft, det vil sige antallet af tilfælde i forhold til befolkningens størrelse, og hvordan den varierer mellem byerne og mellem aldersgrupperne. Antallet af tilfælde er afhængig variabel, og befolkningens størrelse føres ind som offset (på logaritmisk skala):

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

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

Koefficienterne i summary() er på logaritmisk skala. Interceptet (−5,4834) er logaritmen af incidensen i referencegruppen, altså personer på 85 år og ældre i Minneapolis: $e^{-5{,}4834} = 0{,}0042$ tilfælde per person, eller 4,2 tilfælde per 1 000 personer i løbet af studieperioden.

Incidensrate-ratioer med konfidensintervaller

For at få incidensrate-ratioerne eksponentierer man koefficienterne og deres konfidensintervaller. Med pakken broom og funktionen tidy() bliver udskriften en overskuelig tabel:

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

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

Incidensrate-ratioerne med 95 % konfidensinterval (Wald) bliver:

Hvordan fortolkes resultatet?

• By. Incidensrate-ratioen for Dallas er 2,23. Incidensen af hudkræft er altså 2,23 gange højere i Dallas end i Minneapolis (123 % højere) ved samme aldersfordeling. Det er en relativ forskel i rate, ikke 2,23 flere tilfælde per person. Minneapolis er referencekategori og ses derfor ikke i tabellen. Forskellen er forenelig med den højere UV-stråling i Dallas, som ligger betydeligt længere mod syd.
• Alder. Den ældste gruppe (85+) er referencekategori. Alle yngre aldersgrupper har incidensrate-ratioer under 1, det vil sige lavere incidens end de ældste, og ratioen stiger trinvis med alderen. Personer på 15–24 år har en incidens, der er cirka 0,2 % af incidensen hos personer på 85 år og ældre. For 75–84 år omfatter konfidensintervallet 1 (p = 0,30): incidensen adskiller sig ikke sikkert fra den i den ældste gruppe.

Incidensrate-ratioerne kan vises i et forest plot. Skalaen er logaritmisk, så ratioer under og over 1 får samme afstand fra linjen for ingen forskel. Med pakken sjPlot laves et sådant diagram med plot_model(model.1); diagrammet nedenfor viser de samme estimater.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Incidensrate-ratioer (IRR) med 95 % konfidensinterval fra Poissonmodellen for hudkræft: 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": "Incidensrate-ratio (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 grænse", "format": ".3f"}, {"field": "hi", "title": "Øvre grænse", "format": ".3f"}]}}
  ]
}

Passer modellen?

To kontroller bør altid udføres efter en Poissonregression.

Modeltilpasning (goodness of fit). Residualdeviansen (8,26 på 7 frihedsgrader) kan sammenlignes med en chi2-fordeling. En høj p-værdi betyder, at der ikke er tegn på, at modellen passer dårligt:

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

Overdispersion. Forholdet mellem residualdeviansen og frihedsgraderne er 8,26 / 7 = 1,18, altså tæt på 1. Poissonantagelsen om ens middelværdi og varians virker rimelig, og der er ikke behov for kvasi-Poisson eller negativ binomialregression i dette eksempel.

Prædiktioner

Med predict() kan modellen forudsige incidensen for vilkårlige kombinationer af by og aldersgruppe. Da modellen er logaritmisk, eksponentierer man prædiktionen (eller angiver type = "response"). Ved at sætte befolkningen n til 1 får man det forventede antal tilfælde per person, og med n = 1 000 antallet per 1 000 personer. Et antal per 1 000 personer er lettere at fortolke end et lille decimaltal.

# Forventet antal tilfælde per person, Minneapolis, 85 år og ældre
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 ældste aldersgruppe i Minneapolis forventes der altså 4,2 tilfælde per 1 000 personer i løbet af studieperioden. Det tilsvarende tal i Dallas er 4,2 × 2,23 = 9,3 per 1 000.

Prædiktioner for alle rækker i datasættet laves på samme måde:

# Nye datasæt med samme by og aldersgruppe som originalen
newdat1 <- nonmel
newdat2 <- nonmel

# Befolkningen sættes til 1 for at få antal tilfælde per person
newdat1$n <- 1
nonmel$pred.cases.per.one <- predict(model.1, newdat1, type = "response")

# Befolkningen sættes til 1 000 for at få antal tilfælde per 1 000 personer
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- predict(model.1, newdat2, type = "response")

# Forventet antal tilfælde i den faktiske befolkning
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 observerede antal tilfælde (cases). Overensstemmelsen er god i alle grupper, hvilket stemmer overens med den lave residualdevians.

Opsummering

• Poissonregression bruges, når udfaldet er et antal hændelser, og med en offset, når udfaldet er en rate (hændelser per tid, areal eller befolkning).
• Modellen er en GLM med logaritmisk linkfunktion. Eksponentierede koefficienter er incidensrate-ratioer (IRR), der fortolkes multiplikativt.
• Poissonfordelingen forudsætter, at variansen er lig med middelværdien. Kontrollér altid modellens tilpasning og dispersion.
• Ved overdispersion bruges kvasi-Poisson, negativ binomialregression eller robuste standardfejl. Ved overskud af nuller bruges nulinflaterede modeller.