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

Logistisk regression

Logistisk regression laves med funktionen glm( y ~ x, data = data, family = "binomial"), altså en generaliseret linær model funktion hvor family-argumentet er defineret som "binomial". Resultatet kaldes derefter med summary()-funktionen:

model1 <- glm( mors1y ~ diagnosis,
    data = strokedata, family = "binomial")
summary(model1)

Call:
glm(formula = mors1y ~ diagnosis, family = "binomial", data = strokedata)

Deviance Residuals:
Min 1Q Median 3Q Max
-0.9658 -0.5853 -0.5853 -0.5853 1.9229
Coefficients:
EstimateStd. Error z valuePr(>|z|)
(Intercept) -0.5205 0.1972-2.6400.0083 **
diagdictIS -1.1569 0.2170 -5.332 9.69e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

  Null deviance: 973.29 on 1030 degrees of freedom
Residual deviance: 947.29 on 1029 degrees of freedom
AIC: 951.29

Number of Fisher Scoring iterations: 4

I summary() angives coefficienterne på logaritmisk skala og skal derfor transformeres med exp() for at man får odds ratioerne. En lidt lettere og overskuelig måde at få odds ratioerne er med funktionen publish fra pakken {Publish}:

Publish::publish(model1)

Variable Units OddsRatio CI.95 p-value
diagnosis ICH Ref
IS 0.31 [0.21;0.48] <1e-04

Ændring af referencegruppe

Som default sortere r factor variable i alfabetisk orden og det er ikke altid fortolkningsmæssigt det mest korrekte. Man kan derfor have brug for at vælge en anden referencekategori, hvilket kan gøres med funktionen relevel(var, ref = referenceværdi):

model2 <- glm(mors1y ~ relevel(diagnosis, ref = "IS"), data = strokedata, family = "binomial")
Publish::publish(model2)

Variable Units OddsRatio CI.95 p-value
relevel(diagnosis, ref = "IS") IS Ref
ICH 3.18 [2.08;4.87] <1e-04

Justering for kovariate

Modellen kan justeres for kovariate ved at tilføje dem på x-siden af formlen svarende til y = x + a. I nedenstående eksempel laves først en model med et års dødelighed som funktion af rygning:

model3 <- glm(mors1y ~ relevel(smoking, ref = "Never"), data = strokedata, family = "binomial")
Publish::publish(model3)

Variable Units OddsRatio CI.95 p-value
relevel(smoking, ref = "Never") Never Ref
Current 0.54 [0.35;0.83] 0.004904
Previous 0.87 [0.60;1.28] 0.486122

I modellen ser det ud til at rygning beskytter imod at dø, hvilket umiddelbart ikke giver mening. En mulig confounder kunne være alder (altså at forskellen i dødelighed er udtryk for at grupperne har forskellige aldre):

model4 <- glm(mors1y ~ relevel(smoking, ref = "Never") + age, data = strokedata, family = "binomial")
Publish::publish(model4)

Variable Units OddsRatio CI.95 p-value
relevel(smoking, ref = "Never") Never Ref
Current 1.27 [0.78;2.08] 0.3343
Previous 1.01 [0.67;1.53] 0.9634
age 1.11 [1.09;1.14] <1e-04

Resultatet her viser at der er en signifikant association mellem alder og dødelighed, men at associationen mellem rygning og dødelighed ikke længere er signifikant (for nogle af kategorierne). Alder var således en confounder og den tilsyneladende beskyttende effekt af at være nuværende ryger kan derfor forklares med at nuværende rygere er yngre.


Sammenligning af nestede modeller

Model 4 i ovenstående eksempel består af model 3 med en tilføjet kovariate (alder) - og man siger derfor at model 4 er nested i model 3. I modeludviklingen kan det være interessant at undersøge om tilføjelsen af en (eller flere) kovariate bidrager med mere information til modellen.

Summary() funktionen giver to mål man kan bruge til at sammenligne to nestede modeller: AIC og residual deviance.

AIC kan bruges til at lave en uformel sammenligning – idet en lavere AIC tyder på en "bedre" model - altså en model med mere information.


Hvis man vil kalde AIC uden at kalde hele summary()-outputtet, kan man nøjes med post-fixet $aic:

model3$aic
model4$aic
[1] 899.2056
[1] 759.162

AIC for model 4 er lavere end for model 3, og model 4 er derfor en mere præcis/informationsrig model.


Log-likelihood ratio test

Residual deviance kan enten bruges på samme måde som AIC, men kan også bruges til at lave en formel sammenligning ved hjælp af en log-likelihood ratio test.

Residual deviance kan kaldes fra modellen med post-fixet $deviance:

model3$deviance
model4$deviance
[1] 893.2056
[1] 751.162

Testet laves ved at man tager chi2-fordelingen af forskellen mellem deviancerne med funktionen pchisq( chi, df = df, lower.tail = FALSE):

c <- model3$deviance-model4$deviance
c
pchisq(c, df = 1, lower.tail = FALSE)
[1] 142.0435
[1] 9.514076e-33

Udregningen kan også laves med funktionen lrtest( m1, m2), som findes i pakken {lmtest}:

lmtest::lrtest(model3, model4)
Likelihood ratio test

Model 1: mors1y ~ relevel(smoking, ref = "Never")
Model 2: mors1y ~ relevel(smoking, ref = "Never") + age
#DfLogLikDfChisqPr(>Chisq)
1 3 -446.60
2 4-375.581142.04< 2.2e-16***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Bemærk at lmtest bruger en lidt anderledes fremgangmåde til at udregne Chi2-værdien, men at resultatet er det samme:

a <- model3$deviance * -0.5
a
b <- model4$deviance * -0.5
b
-2 * log( exp(a) / exp(b) )
[1] -446.6028
[1] -375.581
[1] 142.0435