Statistik
eksempler i R

Neurologi Neuroanatomi Statistik Home

Deskriptiv


Kvantitative data
Kategoriske data
Intervaller

Analytisk


Sandsynligheder

Kategoriske udfald

Kategoriske eksponeringer
Logistisk regression

Kvantitative udfald

Kvantitative udfald
Linær regression
Korrelationer
Overlevelse
Poisson regression

Tilfældighed


Randomisering

Forskning


PhD thesis



Jacob Liljehult
Klinisk sygeplejespecialist
cand.scient.san, Ph.d.

Neurologisk afdeling
Nordsjællands Hospital

Poisson regression

Poisson regression bruges til at lave intensitets-baserede, marginale modeller, til at sammenligne antallet af tilfælde i to grupper i forhold til mængden af risiko-tid.

Regressionen udføres ved hjælp af en generaliseret linær model med funktionn glm(), med syntaxen:

glm( y ~ x + offset(log(z)), data = data, family="poisson")

Offsetet er logaritmen af den tidsstruktur der er i data - det kan fx være person-risiko-år, altså summen af år hvor alle personer i gruppen er i risiko for udfaldet. Hvis en person bliver fulgt op efter et år har personen bidraget med 1 person-år; mens en person der dør efter 60 dage ikke længere er i risiko efter sin død og derfor kun bidrager med 60/365 person-år til modellen. Events tælles som antallet af events over risiko-perioden og samme person kan derfor godt bidrage med flere events over risiko-personen. Hvis udfaldet er antal epileptiske anfald kan samme patient godt have flere epileptiske anfald over perioden - og hvert anfald tælles derfor som ét event. family-argumentet defineres som "poisson"

Eksempel

Ud fra datasættet strokedata kan man undersøge om det øger risikoen for en ny apopleksi, hvis man tidligere har haft apopleksi. Der bruges følgende variable: prevapo = tidligere apopleksi; recurrence = ny apopleksi; daystorcurrence = antal dage indtil ny apopleksi, dødsfald eller afslutningen af opfølgningen.

Fordi det er en marginal model skal data omstruktureres til et tabel-format, så man får de samlede antal for hver gruppe. Det kan fx gøres med {dplyr}

poisson_df <- strokedata %>%
select(age, prevapo, recurrence, daystorecurrence) %>%
mutate(tidligereapo = prevapo) %>%
group_by(tidligereapo) %>%
summarise(
events = sum(recurrence == "Event"),
personaar = sum(daystorecurrence/365),
rate = events/personaar)
poisson_df
tidligereapo events personaar rate
<fctr> <int><dbl><dbl>
N 73 2133.6164 0.03421421
Y 43 595.2466 0.07223897

Alternativt kan man manuelt indtastning data:

poisson_df <- data.frame(
tidligereapo = c(0,1),
events = c(73,43),
personaar = c(2133.6164,595.2466) )

Her kan vi altså se at blandt personer der tidligere har haft en apopleksi er der 43 nye indlæggelser med apopleksi indenfor en samlet risiko-periode på 595,2466 person-år; hvilket giver en rate på 0,072 per person-år. Blandt personer der ikke tidligere havde haft apopleksi var der 73 nye indlæggelser indenfor en samlet risiko-periode på 2133,2466 person-år; svarende til en rate på 0,034 per person-år.

Poisson-modellen opstilles således:

model1_1 = glm(events ~ factor(tidligereapo) + offset(log(personaar)),
data = poisson_df, family="poisson")

Det 'rå' resultat af modellen kan kaldes med summary(), men resultatet er her på en logaritmisk skala og kan umiddelbart være svært at fortolke:

summary(model1_1)
Call:
glm(formula = events ~ factor(tidligereapo) + offset(log(personaar)),
family = "poisson", data = poisson_df)

Deviance Residuals:
[1] 0 0

Coefficients:
EstimateStd. Errorz valuePr(>|z|)
(Intercept)-3.37510.1170-28.837< 2e-16***
factor(tidligereapo)Y0.74730.19223.8880.000101***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 1.3912e+01 on 1 degrees of freedom
Residual deviance: 5.7732e-15 on 0 degrees of freedom
AIC: 15.734

Number of Fisher Scoring iterations: 2

For at gøre fortolkningen lettere kan man bruge funktionen publish() fra pakken {Publish} til at omregne estimat og konfidens-interval til Hazard-ratio skala:

Publish::publish(model1_1)
VariableUnitsHazardRatioCI.95p-value
factor(tidligereapo)NRef
Y2.11[1.45;3.08]0.0001012

Modellen viser altså at risikoen for at få en ny apopleksi efter udskrivelsen er omkring dobbelt så høj blandt de patienter der tidligere havde haft en apopleksi sammenlignet med de patienter, der ikke tidligere havde haft det.

Justering for alder

Hvis man vil justere modellen for en mulig confounder kan der gøres ved at stratificere tabellen for denne variabel, så man får antal events og person-tid for hver strata. I dette eksempel er der stratificeret for alder ved at opdele alders variablen (som i strokedata er kontinuert) i tre grupper (henholdsvis <60 år/ 60-75 år/ >75 år) - herefter er tabellen er stratificeret for både tidligere apopleksi og alder:

poisson_df <- strokedata %>%
select(age, prevapo, recurrence, daystorecurrence) %>%
mutate(aldersgruppe = factor(case_when(
age < 60 ~ 1,
age >= 60 & age <= 75 ~ 2,
age > 75 ~ 3) ),
tidligereapo = prevapo) %>%
group_by(tidligereapo, aldersgruppe) %>%
summarise(
events = sum(recurrence == "Event"),
personaar = sum(daystorecurrence/365),
rate = events/personaar)
poisson_df
tidligereapoaldersgruppeeventspersonaarrate
<fctr><fctr><int><dbl><dbl>
N113449.808220.02890121
N2351032.468490.03389934
N325651.339730.03838243
Y1592.956160.05378879
Y223273.726030.08402562
Y315228.564380.06562702

For at justere for alder indsætte alders variablen efter x-argumentet (men før offset argumentet):

model1_2 = glm(events ~ factor(tidligereapo) + factor(aldersgruppe) + offset(log(personaar)),
data = poisson_df, family="poisson")
Publish::publish(model1_2)
VariableUnitsHazardRatio CI.95p-value
factor(tidligereapo)NRef
Y2.09[1.43;3.05]0.0001369
factor(aldersgruppe)1 Ref
2 1.29 [0.76;2.20]0.3405091
31.27[0.73;2.22]0.4045718

Den justerede model viser at Hazard-ratioen for tidligere apopleksi ikke ændre sig nævneværdigt efter justering for alder og at der ikke er noget der tyder på en alders gradient.


Goodness-of-fit

Man kan lave en goodness-of-fit analyse for at sammenligne to nestede modeller (altså om tilføjelsen af en justering forbedre modellen signifikant) ved hjælp af en log-likelihood test. Den laves lettest med funktionen lrtest() fra pakken {lmtest}. For at testen virker er det vigtigt at de to modeller, der sammenlignes er nestede og er udregnet på samme dataset.

model1_1 = glm(events ~ factor(tidligereapo) + offset(log(personaar)),
data = poisson_df, family="poisson")
model1_2 = glm(events ~ factor(tidligereapo) + factor(aldersgruppe) + offset(log(personaar)),
data = poisson_df, family="poisson")
lmtest::lrtest(model1_1, model1_2)
Likelihood ratio test

Model 1: events ~ factor(tidligereapo) + offset(log(personaar))
Model 2: events ~ factor(tidligereapo) + factor(aldersgruppe) + offset(log(personaar))
#Df LogLik Df Chisq Pr(>Chisq)
1 2 -14.866
2 4 -14.3642 1.005 0.605

Testet viser en p-værdi på 0,60, hvilket vil sige at tilføjelse af alder ikke har forbedret modellen.