Poissonregression

Introduktion till Poissonregression

Poissonregression används för att modellera antalet gånger en händelse inträffar under en given tid eller inom ett givet område. Utfallet är alltså ett antal (eng. count): 0, 1, 2, 3 och så vidare. För att använda metoden på rätt sätt behöver man förstå några centrala begrepp och den sannolikhetsfördelning som metoden bygger på, nämligen Poissonfördelningen. Kapitlet går igenom begreppen, Poissonfördelningen och modellen, och avslutas med ett fullständigt exempel i R.

Anta att vi undersöker hur många patienter som söker vård på en vårdcentral per dag. Vi vill veta vilka faktorer (prediktorer) som påverkar antalet besök, och vi vill kunna förutsäga (predicera) hur många besök som kan väntas en viss dag. Varje besök är en händelse (eng. event), och antalet besök per dag är ett antal. Sätts antalet i relation till den tid det gäller får man en rate (händelser per tidsenhet): om 10 patienter söker under en dag är raten 10 patienter per dag. Raten kan också gälla ett område, till exempel antal vårdsökande per kvadratkilometer.

Poissonregression kan användas så länge händelserna uppstår genom en Poissonprocess, det vill säga en process som ger upphov till en Poissonfördelning (se nedan). Vad händelsen är spelar ingen roll. Metoden används för att studera antal dödsfall, antal sjukhusinläggningar, antal infektioner, antal kunder eller antal e-postmeddelanden.

Poissonfördelningen

Poissonfördelningen beskrevs av den franske matematikern Siméon Denis Poisson. Den anger sannolikheten för att ett visst antal händelser ska inträffa inom ett givet intervall av tid eller rum, förutsatt att händelserna inträffar oberoende av varandra och med en konstant genomsnittlig takt. Fördelningen har bara en parameter, λ (lambda), som är det förväntade antalet händelser i intervallet. Sannolikheten att observera exakt $k$ händelser är

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

En utmärkande egenskap är att medelvärdet och variansen båda är lika med λ. Det är just denna egenskap som ofta inte håller i verkliga data (se avsnittet om overdispersion).

Figuren visar Poissonfördelningen för tre värden på λ. Fördelningen är diskret: sannolikhet finns bara för heltal. När λ är litet ligger det mesta av sannolikheten nära noll, och noll händelser är vanligt. När λ ökar förskjuts fördelningen åt höger, blir bredare och alltmer symmetrisk, och noll händelser blir osannolikt. Fördelningens topp ligger vid λ eller strax under (för heltal λ är λ och λ − 1 lika sannolika), medan medelvärdet alltid är exakt λ.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Poissonfördelningen för λ = 1, 4 och 10: sannolikheten för varje 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": "Sannolikhet 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": "Förväntat antal"}, {"field": "k", "title": "Antal händelser"}, {"field": "p", "title": "Sannolikhet", "format": ".4f"}]}}
  ]
}

Exempel. Du arbetar på en vårdcentral där i genomsnitt 10 patienter söker per dag, så λ = 10. Antalet varierar från dag till dag. Om patienterna kommer oberoende av varandra kan Poissonfördelningen användas för att beräkna hur sannolikt det är att exakt 5 patienter kommer en viss dag, eller att fler än 5 kommer. I R görs det med dpois() (sannolikheten för exakt k) och ppois() (sannolikheten för k eller färre):

# Sannolikheten att exakt 5 patienter kommer när det förväntade antalet är 10
dpois(5, lambda = 10)
#> [1] 0.03783327

# Sannolikheten att fler än 5 patienter kommer
ppois(5, lambda = 10, lower.tail = FALSE)
#> [1] 0.932914

Exakt 5 patienter inträffar alltså ungefär var 26:e dag (3,8 %), medan fler än 5 patienter kommer drygt 93 % av dagarna.

Antaganden för Poissonfördelningen

• Utfallet $k$ är ett antal händelser inom ett visst intervall och är ett icke-negativt heltal (0, 1, 2, 3 och så vidare).
• Händelserna inträffar oberoende av varandra: att en händelse inträffar påverkar inte när nästa händelse inträffar.
• Den genomsnittliga takten med vilken händelser inträffar är konstant inom intervallet.
• Två händelser kan inte inträffa exakt samtidigt.

Poissonregression är en generaliserad linjär modell

Poissonregression är, som namnet antyder, en typ av regressionsanalys. Den används när den beroende variabeln (utfallet, Y) är ett antal. Poissonregression tillhör de generaliserade linjära modellerna (eng. generalized linear models, GLM).

Vad är en GLM?

Generaliserade linjära modeller är ett ramverk som omfattar flera statistiska modeller, bland annat vanlig linjär regression. I linjär regression modelleras Y direkt som en linjär funktion av prediktorerna ($X_1, X_2, X_3, \ldots$) plus en normalfördelad felterm. Det fungerar dåligt för antal: ett antal kan inte vara negativt, och spridningen ökar ofta med medelvärdet. I en GLM tillåts utfallet i stället följa andra fördelningar ur den så kallade exponentialfamiljen (bland annat Poisson-, binomial-, gamma- och normalfördelningen), och sambandet mellan prediktorerna och utfallets förväntade värde går via en länkfunktion.

Poissonregression använder den naturliga logaritmen som länkfunktion. Det är logaritmen av det förväntade antalet, $\mu = E(Y)$, som modelleras som en linjär funktion av prediktorerna:

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

Till skillnad från linjär regression finns ingen separat felterm i ekvationen: slumpvariationen kring $\mu_i$ beskrivs i stället av Poissonfördelningen. Logaritmen gör att det förväntade antalet alltid blir positivt, $\mu_i = e^{\beta_0 + \beta_1 X_{1i} + \cdots}$, och att prediktorernas effekter blir multiplikativa.

Tolkning av koefficienterna: incidenskvoter (IRR)

Eftersom modellen är logaritmisk tolkas koefficienterna efter exponentiering. $e^{\beta}$ kallas incidenskvot (eng. incidence rate ratio, IRR) och anger med vilken faktor det förväntade antalet (eller raten) ändras när prediktorn ökar med en enhet, eller för en kategori jämfört med referenskategorin, med övriga prediktorer oförändrade. IRR = 1 innebär ingen skillnad, IRR = 2 innebär en fördubblad rate och IRR = 0,5 en halverad rate.

Exempel på frågor som besvaras med Poissonregression

Exempel 1. Hur påverkar utomhustemperaturen och vädret antalet förkylningar som diagnostiseras på en vårdcentral per vecka? Data kan se ut så här (de första 13 veckorna visas):

Vecka, temperatur och väder är prediktorer som används för att förklara och predicera antalet förkylningar. Modellen skrivs

$$\log(\text{förväntat antal förkylningar}_i) = \beta_0 + \beta_1,\text{vecka}_i + \beta_2,\text{temperatur}_i + \beta_3,\text{väder}_i$$

där $\beta_0$ är interceptet (logaritmen av det förväntade antalet när alla prediktorer är noll eller i sin referenskategori) och $\beta_1, \beta_2, \beta_3$ är prediktorernas koefficienter. Väder är en kategorisk variabel och får i praktiken en koefficient per kategori utom referenskategorin.

Exempel 2. Hur många kunder genomför ett köp i en butik per dag, och vad påverkar det?

Exempel 3. Hur många personer reser med en viss busslinje till arbetet, och hur påverkas antalet av väder och veckodag?

Offset: att analysera rater

Ofta är det inte antalet i sig som är intressant utan antalet i relation till en exponering: tid, area, distans eller befolkningsstorlek. Två exempel:

• Du studerar dödsfall och vet hur länge varje individ har följts upp (exponeringstiden). Då är det incidensen, antal händelser dividerat med exponeringstiden, som är intressant. Termen event rate används ofta för antal händelser per tidsenhet.
• Du studerar antalet djurarter i olika geografiska områden av olika storlek. Då är det antalet arter per area som är intressant.

Exponeringen förs in i modellen som en offset: logaritmen av exponeringen läggs till på högersidan med en koefficient som är fixerad till 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 då raten $\mu_i / t_i$, och $e^{\beta}$ blir en kvot mellan rater. I R anges offseten med funktionen offset() i formeln, eller med argumentet offset:

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

Observera att offseten ska anges på logaritmskalan (log(area)), eftersom modellen är logaritmisk.

Antaganden för Poissonregression

• Utfallet är ett antal (icke-negativa heltal).
• Observationerna är oberoende av varandra.
• Logaritmen av det förväntade antalet är en linjär funktion av prediktorerna. Kontinuerliga prediktorer kan behöva transformeras eller modelleras flexibelt (till exempel med splines) om sambandet inte är log-linjärt.
• Variansen är lika med medelvärdet (equidispersion). Detta antagande bryts ofta, se nästa avsnitt.
• Exponeringen (tid, area, befolkning) är korrekt angiven som offset när raten är det intressanta.

Overdispersion och nollinflation

Poissonmodellen förutsätter att variansen är lika med medelvärdet. I verkliga data är variansen ofta större än medelvärdet, vilket kallas overdispersion. Det beror ofta på att viktiga prediktorer saknas i modellen, att observationerna inte är oberoende (till exempel upprepade mätningar på samma individ) eller att individerna skiljer sig åt på sätt som modellen inte fångar. Punktskattningarna påverkas sällan nämnvärt, men standardfelen blir för små, konfidensintervallen för smala och p-värdena för låga.

Overdispersion kan bedömas genom att dividera Pearsons chi2-statistika (eller residualdeviansen) med antalet frihetsgrader. En kvot nära 1 talar för att Poissonantagandet håller. En kvot klart över 1 (tumregel: över 1,5) talar för overdispersion.

# Kvoten mellan Pearsons chi2 och antalet frihetsgrader
sum(residuals(model, type = "pearson")^2) / df.residual(model)

Vid overdispersion finns flera alternativ:

# 1. Kvasi-Poisson: samma skattningar, standardfelen skalas upp med dispersionen
model.qp <- glm(y ~ x1 + x2 + offset(log(tid)), family = quasipoisson, data = data)

# 2. Negativ binomialregression: modellerar extra variation med en egen parameter
library(MASS)
model.nb <- glm.nb(y ~ x1 + x2 + offset(log(tid)), data = data)

# 3. Robusta (sandwich) standardfel för en vanlig Poissonmodell
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model, type = "HC0"))

Ett annat vanligt problem är överskott av nollor (eng. zero inflation). Det uppstår när nollorna kommer från två olika processer. Tänk dig att vi studerar antalet cigaretter som röks per dag i en skolklass där halva klassen röker. Bland rökarna följer antalet cigaretter en Poissonprocess, medan icke-rökarna alltid röker noll cigaretter. Data innehåller då betydligt fler nollor än Poissonfördelningen förutsäger. I sådana situationer passar en nollinflaterad modell bättre, till exempel zeroinfl() i paketet pscl, som modellerar sannolikheten för en "strukturell" nolla och antalet bland övriga i två delar:

library(pscl)
# Före | : modellen för antalet. Efter | : modellen för de strukturella nollorna.
model.zip <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "poisson")
model.zinb <- zeroinfl(cigaretter ~ alder + kon | foraldrar_roker, data = data, dist = "negbin")

Exempel: hudcancer i Minneapolis och Dallas

Data kommer från en klassisk studie av icke-melanom hudcancer hos vita invånare i fyra områden 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). Här används två av städerna, Minneapolis och Dallas. Data består av antalet nya fall (cases) och befolkningens storlek (n) i åtta åldersgrupper i varje stad. Datamängden är liten och kan läsas in direkt i R med read.table(). Variablerna u1 till u7 är indikatorvariabler för åldersgrupperna och behövs inte i analysen nedan, där åldersgruppen i stället skapas som en faktor.

# Läser in data direkt 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
")

# Skapar variabler för åldersgrupp och stad.
# age.range är en faktor med 8 nivåer där den högsta åldersgruppen (85+) är referensnivå.
# within() används för att slippa skriva nonmel$ framför varje 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 är intresserade av incidensen av hudcancer, det vill säga antalet fall i förhållande till befolkningens storlek, och hur den skiljer sig mellan städerna och mellan åldersgrupperna. Antalet fall är beroende variabel och befolkningens storlek förs in som offset (på logaritmskalan):

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

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

Koefficienterna i summary() är på logaritmskalan. Interceptet (−5,4834) är logaritmen av incidensen i referensgruppen, alltså personer 85 år och äldre i Minneapolis: $e^{-5{,}4834} = 0{,}0042$ fall per person, eller 4,2 fall per 1 000 personer under studieperioden.

Incidenskvoter med konfidensintervall

För att få incidenskvoterna exponentierar man koefficienterna och deras konfidensintervall. Med paketet broom och funktionen tidy() blir utskriften en prydlig tabell:

library(broom)
# exponentiate = TRUE ger incidenskvoter (IRR) i stället för log-skattningar
tidy(model.1, exponentiate = TRUE, conf.int = TRUE)

# Samma sak i bas-R, med Wald-baserade konfidensintervall
exp(cbind(IRR = coef(model.1), confint.default(model.1)))

Incidenskvoterna med 95 % konfidensintervall (Wald) blir:

Hur tolkas resultatet?

• Stad. Incidenskvoten för Dallas är 2,23. Incidensen av hudcancer är alltså 2,23 gånger högre i Dallas än i Minneapolis (123 % högre) vid samma åldersfördelning. Det är en relativ skillnad i rate, inte 2,23 fler fall per person. Minneapolis är referenskategori och syns därför inte i tabellen. Skillnaden är förenlig med den högre UV-strålningen i Dallas, som ligger betydligt längre söderut.
• Ålder. Den äldsta gruppen (85+) är referenskategori. Alla yngre åldersgrupper har incidenskvoter under 1, det vill säga lägre incidens än de äldsta, och kvoten stiger stegvis med åldern. Personer 15–24 år har en incidens som är cirka 0,2 % av incidensen hos personer 85 år och äldre. För 75–84 år inkluderar konfidensintervallet 1 (p = 0,30): incidensen skiljer sig inte säkert från den i den äldsta gruppen.

Incidenskvoterna kan visas i ett forestdiagram. Skalan är logaritmisk, så att kvoter under och över 1 får samma avstånd från linjen för ingen skillnad. Med paketet sjPlot görs ett sådant diagram med plot_model(model.1); diagrammet nedan visar samma skattningar.

{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Incidenskvoter (IRR) med 95 % konfidensintervall från Poissonmodellen för hudcancer: stad och åldersgrupp.",
  "width": "container",
  "height": 300,
  "config": {"view": {"stroke": null}},
  "data": {
    "values": [
      {"term": "Dallas (jämfört 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": "Incidenskvot (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": "Jämförelse"}, {"field": "irr", "title": "IRR", "format": ".3f"}, {"field": "lo", "title": "Nedre gräns", "format": ".3f"}, {"field": "hi", "title": "Övre gräns", "format": ".3f"}]}}
  ]
}

Passar modellen?

Två kontroller bör alltid göras efter en Poissonregression.

Anpassning (goodness of fit). Residualdeviansen (8,26 på 7 frihetsgrader) kan jämföras med en chi2-fördelning. Ett högt p-värde betyder att det inte finns tecken på att modellen passar dåligt:

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

Overdispersion. Kvoten mellan residualdeviansen och frihetsgraderna är 8,26 / 7 = 1,18, alltså nära 1. Poissonantagandet om lika medelvärde och varians verkar rimligt, och det finns inget behov av kvasi-Poisson eller negativ binomialregression i detta exempel.

Prediktioner

Med predict() kan modellen förutsäga incidensen för valfria kombinationer av stad och åldersgrupp. Eftersom modellen är logaritmisk exponentierar man prediktionen (eller anger type = "response"). Genom att sätta befolkningen n till 1 får man det förväntade antalet fall per person, och med n = 1 000 antalet per 1 000 personer. Ett antal per 1 000 personer är lättare att tolka än ett litet decimaltal.

# Förväntat antal fall per person, Minneapolis, 85 år och ä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 äldsta åldersgruppen i Minneapolis väntas alltså 4,2 fall per 1 000 personer under studieperioden. Motsvarande siffra i Dallas är 4,2 × 2,23 = 9,3 per 1 000.

Prediktioner för alla rader i datamängden görs på samma sätt:

# Nya datamängder med samma stad och åldersgrupp som originalet
newdat1 <- nonmel
newdat2 <- nonmel

# Befolkningen sätts till 1 för att få antal fall per person
newdat1$n <- 1
nonmel$pred.cases.per.one <- predict(model.1, newdat1, type = "response")

# Befolkningen sätts till 1 000 för att få antal fall per 1 000 personer
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- predict(model.1, newdat2, type = "response")

# Förväntat antal fall i den faktiska 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

Kolumnen pred.cases kan jämföras direkt med det observerade antalet fall (cases). Överensstämmelsen är god i alla grupper, vilket stämmer med den låga residualdeviansen.

Sammanfattning

• Poissonregression används när utfallet är ett antal händelser, och med en offset när utfallet är en rate (händelser per tid, area eller befolkning).
• Modellen är en GLM med logaritmisk länkfunktion. Exponentierade koefficienter är incidenskvoter (IRR) som tolkas multiplikativt.
• Poissonfördelningen förutsätter att variansen är lika med medelvärdet. Kontrollera alltid modellens anpassning och dispersion.
• Vid overdispersion används kvasi-Poisson, negativ binomialregression eller robusta standardfel. Vid överskott av nollor används nollinflaterade modeller.