Avsnitt 1, Ämne 1
Startad

Poisson regression

Poisson-regression (Poissonregression)

Introduktion till Poisson-regression

Poisson-regression används för att modellera antalet gånger en händelse inträffar inom en given tidsperiod eller plats. För att förstå och använda Poisson-regression är det viktigt att begripa de termer som används i dessa sammanhang. Det är också viktigt att ha förståelse för den sannolikhetsfördelning som ligger till grund för Poisson-regression, nämligen Poissonfördelning (Poissondistribution). Vi inleder denna diskussionen med viktiga termer. Därnäst diskuteras Poissonfördelning med flera exempel och slutligen illustreras ett praktiskt exempel med beräkningar i R.

Ponera att vi undersöker hur många patienter som söker vård på en vårdcentral per dag. Vi är intresserade av att undersöka vilka prediktorer som påverkar antalet besök per dag, liksom vill vi kunna predicera (förutsäga) hur många besök som komma skall under en dag då vissa omständigheter råder. Vi är alltså intresserade av antal patienter som söker per dag. Detta antal kallas på engelska count. Vi skulle kunna betrakta varje besök på vårdcentralen som en händelse (eng. event). Om vi undersöker antalet patienter som söker vård under en given tidsperiod (exempelvis under en dag) så kan antalet specificeras som en hastighet (eng. rate); om 10 patienter söker vård under en dag så är rate 10 patienter/dag. Man skulle dock även kunna analysera antalet patienter som söker vård per kvadratkilometer i ett land. Om 100 patienter söker vård i en stad som är 100 kvadratkilometer stor, så är rate 100 patienter/100 kvadratkilometer. Poisson-regression används alltså för att undersöka antal gånger en händelse inträffar i ett givet tidsintervall eller plats.

Man kan använda Poisson-regression så länge händelsen följer en Poissonprocess. Det innebär att den processen (oavsett vilken den är) som ger upphov till en händelse måste resultera i en Poissonfördelning (se nedan). Själva händelsens natur har ingen nämnvärd betydelse, så länge processen resulterar i en Poissonfördelning och antagandena för Poisson-regression är uppfyllda (se nedan). Det innebär att Poisson-regression kan användas för att studera antalet dödsfall, antalet e-post, antal kunder eller antal datorer.

Regression (regressionsanalys) med Poisson

Poisson-regression är, som namnet antyder, en typ av regressionsanalys. Denna modell används för att predicera utfall (beroende variabel) som är ett antal (eng. count). Poisson-regression lämpar sig alltså då den beroende variabeln utgörs av ett antal som har en Poissonfördelning (eng. Poisson distribution). En förutsättning för att använda Poisson-regression är att logaritmen av det förväntade värdet på den beroende variabeln (Y) kan modelleras med en linjär funktion. Poissonregression är en typ av generaliserad linjär modell (eng. generalized linear modell, GLM).

Vad är GLM (Generalized Linear Model)?

Generaliserade linjära modeller (GLM) är ett matematiskt ramverk som omfattar flera statistiska modeller och där klassiska modeller, exempelvis linjär regression, ingår. GLM, som är en matematisk generalisering av linjär regression, skapades för att erbjuda en allmän och flexibel statistisk modell som inte är lika begränsad (av svåruppfyllda matematiska antaganden) som enskilda modeller. I vanlig linjär regression antar man att den beroende variabeln (Y) kan modelleras med hjälp av ett antal prediktorer (X1, X2, X3, ...) samt en felterm (ε, residualer). Linjär regression kräver att residualerna är oberoende av varandra och normalfördelade. I generaliserade linjära modeller tillåts andra fördelningar än normalfördelningen. De fördelningar som tillåts ingår i exponentialfamiljen, som innehåller Poissonfördelning, binomialfördelning, normalfördelning och gamma-fördelning. I GLM används en vanlig linjär funktion av prediktorerna (X1, X2, X3, ...) för att skapa en funktion av Y, men sambandet mellan prediktorerna och Y går via en länkfunktion som kan vara icke-linjär. Syftet med en icke-linjär länkfunktion är att öka flexibiliteten (användbarheten) i modellen. Poisson-regression nyttjar logaritmen som länkfunktion, vilket innebär att en linjär kombination av prediktorerna (X1, X2, X3, ...) används för att predicera logaritmen på Y.

Exempel då Poisson-regression används

Poisson-regression användas för att studera data där utfallet (Y) är ett antal.

Exempel 1. Man kan undersöka hur temperaturen utomhus samt vädret påverkar antalet förkylningar som diagnostiseras på en vårdcentral under årets veckor. Data kan se ut som följer (första 13 veckorna visas).

Vecka TemperaturVäderAntal förkylningar
1-10A45
2-8B34
3-9B18
4-6A30
5-1C33
6-3A29
7-2D45
81C40
90A22
102E18
115E19
126A15
139B10

Vecka, temperatur och väder är prediktorer (X) som används för att predicera antal förkylningar (Y). Matematiskt kan detta beskrivas som följer:

Antal förkylningari = β0 + β1Xi + β2Xi + β3Xi + ei

β0 är interceptet (där regressionslinjen skär y-axeln).

β1β2 och βär koefficienter för prediktorerna.

ei är modellens felmarginal (residual error [residualer]).

Exempel 2. Man kan använda Poisson-regression för att studera antal personer som kommer genomföra ett köp i en butik.

Exempel 3. Man kan använda Poisson-regression för att studera hur många personer som kommer ta bussen till jobbet.

Poissonfördelning (Poisson-distribution)

Poissonfördelning beskrevs av den franske matematikern Siméon Denis Poisson. Poissonfördelning beskriver sannolikheten för att ett visst antal händelser skall inträffa inom ett givet intervall av tid eller rum, förutsatt att händelserna inträffar med konstant hastighet och varje händelse skall vara oberoende av tiden som förlöpt sedan den föregående händelsen. Oftast illustreras Poissonfördelningen genom att använda exempel med antal händelser inom en given tidsperiod, men Poissonfördelningen fungerar likväl för antal händelser inom en given area, volym eller avstånd.

Ponera att vi analyserar antal mejl (e-post) som ankommer varje dag till universitetet. I genomsnitt ankommer 12000 mejl dagligen. Om ankomst av ett mejl inte påverkar ankomsttiden för nästa mejl, samt om mejlen ankommer oberoende av varandra, så kan man anta att Poissonfördelningen är lämplig för att analysera mejlen.
Ponera att vi analyserar antal patienter som söker på vårdcentralen varje dag. I genomsnitt kommer 10 patienter mer dag. Det är rimligt att utgå ifrån att patienterna inte påverkar varandra, så till vida att den ena patientens ankomsttid inte påverkar nästkommande patient. Även i det fallet är Poissonfördelningen lämplig.

I denna figuren presenteras 9 olika Poissonfördelningar baserat på λ, som är det förväntade antalet händelser (dvs genomsnittligt antal händelser inom en viss tidsperiod/plats). λ = 1 innebär att vi förväntar oss att 1 händelse skall inträffa. λ = 9 innebär att vi förväntar oss att 9 händelser skall inträffa. En händelse kan vara ett dödsfall, eller ankomst av e-post, eller antal köp som genomförs i en butik. X-axeln visar k, som är antalet händelser som inträffar. Y-axeln visas sannolikheten för att ett givet antal händelser skall inträffa.
I denna figuren presenteras 9 olika Poissonfördelningar baserat på λ, som är det förväntade antalet händelser (dvs genomsnittligt antal händelser inom en viss tidsperiod/plats). λ = 1 innebär att vi förväntar oss att 1 händelse skall inträffa. λ = 9 innebär att vi förväntar oss att 9 händelser skall inträffa. En händelse kan vara ett dödsfall, eller ankomst av e-post, eller antal köp som genomförs i en butik. X-axeln visar k, som är antalet händelser som inträffar. Y-axeln visas sannolikheten för att ett givet antal händelser skall inträffa.

För att förklara Poissonfördelningen utgår vi från ovanstående figur. I figuren visas 9 olika Poissonfördelningar, vilka skiljer sig åt i λ, som är det genomsnittliga antalet händelser som inträffar inom ett givet intervall (i tid eller rum). λ kallas även event rate och kan ses som det förväntade antalet fall. Om vi förväntar oss att ett fåtal händelser skall inträffa, t ex om λ=1, så är det sannolikt att antalet händelser kommer infinna sig mellan 0 och 5. Om λ ökas till 10 (dvs vi förväntar oss 10 händelser), så blir noll osannolikt och istället blir ett större intervall sannolikt (värden från 1 till 18). Detta betyder att när λ är lågt så föreligger hög sannolikhet för att noll händelser skall inträffa, och allteftersom λ ökar (dvs händelsen blir vanligare) så förskjuts kurvan åt höger och noll blir mindre sannolikt.

Exempel: Du jobbar på en vårdcentral (alternativt i en klädbutik) och i genomsnitt kommer 10 patienter (alternativt kunder) per dag. Det förväntade värdet (λ) blir då 10. Det kommer dock finnas en viss spridning, så till vida att antalet patienter kan vara lägre eller högre än 10. Om patienternas ankomst till vårdcentralen är oberoende av varandra så kan Poissonfördelningen användas för att exempelvis beräkna hur sannolikt det är att fler än 5 patienter, eller exakt 5 patienter skall komma. I detta fallet är alltså medelvärdet (λ) ett värde som förtäljer det genomsnittliga antalet händelser. Medelvärdet i varje Poissonfördelning är helt enkelt kurvans apex (högsta punkten).

Antaganden för Poissonfördelningen och Poissonregression

  • Poissonfördelning är lämpligt om följande antaganden är uppfyllda.
  • k är antalet gånger en händelse inträffar inom ett visst intervall och k är ett heltal (0, 1, 2, 3, ...)
  • En händelse (att en händelse inträffar) påverkar inte när nästa händelse inträffar. Händelserna inträffar oberoende av varandra.
  • Takten/hastigheten med vilken händelser inträffar är konstant. Det innebär att hastigheten med vilken händelser inträffar får inte vara långsammare eller snabbare inom givna intervallet.
  • Två, eller flera händelser, får inte inträffa exakt samtidigt.

Användning av offset för att analysera en rate

Poisson-regression används alltså för att studera ett utfall som är ett antal, t ex antal dödsfall inom en viss period. I vissa situationer är det viktigt att antalet sätts i relation till en tidsperiod, plats, distans eller volym. Två exempel följer: Exempel 1. Du studerar en händelse, som exempelvis död. För varje individ har du uppgift om hur länge du följer upp individen (tiden då den är exponerad). I detta fall kan Poissonregression användas för att studera en incidens (antal händelser dividerat med exponeringstiden). Termen event rate används ofta för att beskriva antalet händelser per tidsperiod. Exempel 2: Du studerar hur många djurarter som lever inom ett visst geografiskt området. Det geografiska området kan mätas i form av area (t ex hektar). I detta fall kan Poissonregression användas för att studera antalet djurarter som kan förväntas per area. För att på ett korrekt sätt hantera exponeringstid (exempel 1) och area (exempel 2), vilka kallas exponeringar, så inkluderar man en offset i modellen. I R är det med funktionen offset() som man anger vilken variabel som är modellens offset (förutsatt att en sådan finns). Se nedan för exempel.

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

Overdispersion och zero inflation (zero inflated model)

Poissonregression har en särskilt egenskap, nämligen att medelvärdet (λ) skall vara ekvivalent med variansen (λ). I vissa fall kan variansen överstiga medelvärdet, varvid fenomenet overdispersion uppstår. Ofta brukar detta förklaras av att man saknar viktiga prediktorer men andra faktorer kan också vara orsaken. I dessa fall kan man istället använda negativ binomial distribution, som omtalas annorstädes. Ett annat vanligt problem med Poissonregression är överflöd av nollor (zero inflation). I situationer där det finns överflöd av nollor anses två processer göra sig gällande, så till vida att överflödet av nollor orsakas av en process som är separat från processen som ger upphov till Poissonfördelningen. Föreställ dig att vi studerar antalet cigaretter som röks under en dag i en skolklass där halva klassen röker. Då kommer gruppen som röker generera en Poissonprocess medan gruppen som inte röker kommer generera en separat process, vilken genererar ett överflöd av nollor (inga rökta cigaretter). Situationer där det finns två separata processer och där den ena genererar överflöd av nollor är också olämpliga för Poissonregression. Istället kan man då använda zero inflated models, som omtalas annorstädes.

Poissonregression: exempel på hudcancer i två amerikanska städer

Data kommer från “Non-melanoma skin cancer among Caucasians in four areas of the United States” (www.pubmed.org/4422100). Exemplet återfinns även i Applied Regression Analysis and Multivariable Methods, 4th Edition. Data är tämligen litet, varför det kan läsas in direkt i R via funktionen read.table().

# Skapar data manuellt

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 åldersgrupper och städer
# Först skapas variabeln age.range som är en faktor med 8 nivåer, där högsta ålder sätts till referensnivå (med funktionen "relevel()"").
# funktionen "within" används för att inte behöva använda "$" för att specificera variablerna.
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"))
})

Skapar Poissonmodeller

Ovanstående dataset skall nu användas för att göra en Poissonregression. Vi är intresserade av antalet cancerfall i olika ålderskategorier i två amerikanska städer (Minneapolis och Dallas). Vi är alltså intresserade av en "rate", vilken beskriver antalet cancerfall/invånarantal. För att detta skall fungera i R så använder vi antal cancerfall som beroende variabel (Y) och invånarantal som offset i modellen, enligt följande:

# offset kan anges på följande vis
model.1 <- glm(cases ~ city + age.range + offset(log(n)), family = poisson(link = "log"), data = nonmel)

# offset kan också anges på detta sättet
model.1 <- glm(cases ~ city + age.range, offset = log(n), family = poisson(link = "log"), data = nonmel)

# Resultat
summary(model.1)

Resultat

## 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

Alternativ till “summary()”

För att få en mer städad utskrift att resultaten så kan du använda paketet “broom” och funktionen “tidy()”, enligt följande

# installera "broom"
install.packages("broom")

# aktivera broom
library(broom)

tidy(model.1, exp=T) # exp=T för att exponentiera resultatet, eftersom Poissonregression användar logaritm som länkfunktion

Resultat

##             term    estimate  std.error  statistic      p.value
## 1    (Intercept) 0.004155093 0.10367655 -52.889690 0.000000e+00
## 2     cityDallas 2.234233723 0.05220475  15.398950 1.663497e-53
## 3 age.range15_24 0.002082493 0.45774051 -13.488406 1.830153e-41
## 4 age.range25_34 0.028897843 0.16748387 -21.160176 2.223891e-99
## 5 age.range35_44 0.097605351 0.12746646 -18.254394 1.909191e-74
## 6 age.range45_54 0.206182619 0.11383370 -13.871051 9.487693e-44
## 7 age.range55_64 0.337260039 0.11090787  -9.800035 1.125469e-22
## 8 age.range65_74 0.589325065 0.10861972  -4.868153 1.126463e-06
## 9 age.range75_84 0.890782074 0.11094746  -1.042435 2.972102e-01

Hur tolkas resultatet?

När city = Dallas så har vi 2.23 i estimat, vilket betyder 2.23 fler fall per person i Dallas, jämfört med Minneapolis. Eftersom Minneapolis är referenskategori för “city” så redovisas den inte i ovanstående tabell. För “age.range” så är den högsta ålderskategorin referensgrupp (varför den inte syns i utskriften). Vi ser att alla åldersgrupper har lägre än 1.0, vilket innebär att de har lägre rate än den högsta ålderskategorin. Vill vi få ut en graf direkt med alla incidence rate ratios, så kan vi använda funktionen sjp.glm(), som följer:

#install.packages("sjPlot")
library(sjPlot)
sjp.glm(model.1)

Nu ska vi predicera antal fall per person, per 1000 personer, per riktig populationsstorlek. Detta görs genom att använda “predict()” i R. Vi börjar med att predikta antalet fall per person i Minneapolis i åldrarna 85+. Resultatet blir 0.004155093 vilket innebär att sjukdomen inträffar cirka 0.004155093 gånger per person. Detta är inte särskitl pedagogiskt, och därför är det bättre att predikta antalet fall per 1000 personer, vilket görs i nästa steg.

exp(predict(model.1, newdata = data.frame(city = "Minneapolis", age.range = "85+", n = 1)))

Resultat

## 0.004155093

Skapar ett nytt data set för att predikta

# nya datasettet innehåller city och age.range från ordinarie fil
newdat1 <- nonmel

# kopierar filen
newdat2 <- newdat1

# Vi sätter antal personer till 1, för att predikta antalet fall per person
newdat1$n <- 1
nonmel$pred.cases.per.one <- exp(predict(model.1, newdat1))

# Nu prediktar vi antal fall per 1000 personer
newdat2$n <- 1000
nonmel$pred.cases.per.thousand <- exp(predict(model.1, newdat2))

# Nu prediktar vi antal fall i den faktiska populationen
nonmel$pred.cases <- exp(predict(model.1))

Resultaten från ovanstående prediktioner

##    cases        city u1 u2 u3 u4 u5 u6 u7      n age.range
## 1      1 Minneapolis  1  0  0  0  0  0  0 172675     15_24
## 2     16 Minneapolis  0  1  0  0  0  0  0 123065     25_34
## 3     30 Minneapolis  0  0  1  0  0  0  0  96216     35_44
## 4     71 Minneapolis  0  0  0  1  0  0  0  92051     45_54
## 5    102 Minneapolis  0  0  0  0  1  0  0  72159     55_64
## 6    130 Minneapolis  0  0  0  0  0  1  0  54722     65_74
## 7    133 Minneapolis  0  0  0  0  0  0  1  32185     75_84
## 8     40 Minneapolis  0  0  0  0  0  0  0   8328       85+
## 9      4      Dallas  1  0  0  0  0  0  0 181343     15_24
## 10    38      Dallas  0  1  0  0  0  0  0 146207     25_34
## 11   119      Dallas  0  0  1  0  0  0  0 121374     35_44
## 12   221      Dallas  0  0  0  1  0  0  0 111353     45_54
## 13   259      Dallas  0  0  0  0  1  0  0  83004     55_64
## 14   310      Dallas  0  0  0  0  0  1  0  55932     65_74
## 15   226      Dallas  0  0  0  0  0  0  1  29007     75_84
## 16    65      Dallas  0  0  0  0  0  0  0   7583       85+
##    pred.cases.per.one pred.cases.per.thousand pred.cases
## 1        8.652950e-06              0.00865295   1.494148
## 2        1.200732e-04              0.12007322  14.776810
## 3        4.055593e-04              0.40555928  39.021292
## 4        8.567079e-04              0.85670789  78.860818
## 5        1.401347e-03              1.40134673 101.119779
## 6        2.448700e-03              2.44870028 133.997777
## 7        3.701282e-03              3.70128210 119.125764
## 8        4.155093e-03              4.15509270  34.603612
## 9        1.933271e-05              0.01933271   3.505852
## 10       2.682716e-04              0.26827163  39.223190
## 11       9.061142e-04              0.90611423 109.978708
## 12       1.914086e-03              1.91408567 213.139182
## 13       3.130936e-03              3.13093612 259.880221
## 14       5.470969e-03              5.47096874 306.002223
## 15       8.269529e-03              8.26952928 239.874236
## 16       9.283448e-03              9.28344824  70.396388