Die binäre logistische Regression ist immer dann zu rechnen, wenn die abhängige Variable nur zwei Ausprägungen hat, also binär bzw. dichotom ist. Es wird dann die Wahrscheinlichkeit des Eintritts bei Ändern der unabhängigen Variable geschätzt. Die Schätzung der Wahrscheinlichkeit ist neben der binären Codierung der wesentliche Unterschied zur einfachen Regression.
Im Vorfeld der Regressionsanalyse kann zudem eine Filterung vorgenommen werden, um nur einen gewissen Teil der Stichprobe zu untersuchen, bei dem man am ehesten einen Effekt erwartet. Ist die abhängige Variable metrisch oder quasi-metrisch, kann eine lineare Regression (einfach oder multipel) gerechnet werden.
Inhaltsverzeichnis
1 Voraussetzungen der binär logistischen Regression
Die wichtigsten Voraussetzungen sind:
- die y-Variable ist binär codiert (0 und 1 als Ausprägungen),
- metrisch skalierte x-Variable(n), ordinal funktioniert auch, nominal skalierte x-Variablen sind als “Faktor” in R zu definieren
- keine hoch korrelierten x-Variablen
2 Durchführung der binär logistischen Regression in R
Es gibt verschiedene Möglichkeiten und Pakete hierfür. Ich werde aber zum besseren Verständnis für das Vorgehen mit Bordmitteln arbeiten.
2.1 Testmodell definieren
Das Testmodell beinhaltet die zu prüfenden Hypothesen bzw. die Variablen. In die glm()-Funktion wird zuerst die y-Variable gegeben, gefolgt von ~ und den x-Variablen.
Das Argument family = binomial() definiert das Verwenden einer logit-Funktion, also der natürliche Logarithmus einer Chance (“Odd’s”).
# Testmodell
model1 <- glm(Krankheit ~ Geschlecht + BMI, data = df, family = binomial())
summary(model1)
2.2 Ergebnisübersicht des Testmodells
Die Ergebnisse des Testmodells sehen wie folgt aus:
Call:
glm(formula = Krankheit ~ Geschlecht + BMI, family = binomial(),
data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4109 -0.6445 -0.3351 0.1878 2.3677
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.5638 3.9560 -2.923 0.00347 **
Geschlecht -1.4113 0.8375 -1.685 0.09194 .
BMI 0.4863 0.1700 2.861 0.00422 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 59.945 on 50 degrees of freedom
Residual deviance: 39.818 on 48 degrees of freedom
AIC: 45.818
Number of Fisher Scoring iterations: 5
3 Interpretation der binär logistischen Regression in R
3.1 Omnibus-Test
Bevor auf die Koeffizienten geschaut wird, sollte man sich versichern, dass das Modell besser als das sog. "Nullmodell" ist. Salopp formuliert könnte man sagen, dass geprüft wird, ob das Testmodell das Vorhandensein der Krankheit besser vorhersagt als ein Münzwurf.
Hierfür wird die sog. Deviance bzw. die Differenz der Null deviance und der Residual deviance herangezogen und der Omnibus-Test durchgeführt.
Null deviance: 59.945 on 50 degrees of freedom
Residual deviance: 39.818 on 48 degrees of freedom
Diese Deviance-Differenz wird in Verbindung mit der Freiheitsgrade-Differenz über die Chi-Quadrat-Verteilung ermittelt. Das ist der sog. Omnibus-Test.
Hierzu bildet man die Differenz aus der Null deviance (aus dem Nullmodell) und der Residual deviance (aus dem Testmodell):
# Omnibus-Test
modelchi <- model1$null.deviance - model1$deviance
Das Gleiche macht man für die Freiheitsgrade:
chidf <- model1$df.null - model1$df.residual
Und schließlich wird der p-Wert über 1-pchisq()-Funktion ermittelt. Dieser dient zur Prüfung der Nullhypothese, dass sich Nullmodell und Testmodell voneinander unterscheiden errechnet.
chisqp <- 1-pchisq(modelchi, chidf)
chisqp
[1] 4.26029e-05
Im Beispiel ist dieser p-Wert 4.26e-05 (=0,0000426) und die Nullhypothese wird verworfen. Das Testmodell ist demnach besser, also das Nullmodell und mit der Interpretation kann fortgefahren werden.
3.2 Testmodell interpretieren - Koeffizienten
Ist das Modell bzw. der Omnibus-Test signifikant, kann mit der Interpretation der Koeffizienten begonnen werden. Hierzu noch mal der Ausschnitt aus den Ergebnissen:
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.5638 3.9560 -2.923 0.00347 **
Geschlecht -1.4113 0.8375 -1.685 0.09194 .
BMI 0.4863 0.1700 2.861 0.00422 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
1: p-Werte betrachten
Zunächst ist von euch das Alphaniveau, zu dem getestet wird, festzulegen. Typisch (aber nicht unumstritten) sind 5% (=0.05) oder 1% (=0.01). Die p-Werte in der letzten Spalte sind dem Alpha jeweils gegenüberzustellen.
- Der Intercept wird nicht interpretiert.
- Der p-Wert des Geschlechts ist p = 0.09194. Die Nullhypothese für diesen Koeffizienten kann nicht verworfen werden.
- Der p-Wert des BMI ist p = 0.00422. Die Nullhypothese für diesen Koeffizienten wird verworfen. Ein Effekt ist in der Stichprobe beobachtbar.
2: Koeffizienten betrachten
Anschließend kann auf den Koeffizienten in der ersten Spalte geschaut werden. Da nur der BMI einen niedrigen p-Wert aufweist, kann auch nur dieser interpretiert werden. Er ist positiv, was bedeutet, das mit einer Zunahme des BMI die Chance auf das Vorhandensein der Krankheit steigt.
- Für ordinal-, intervall- und verhältnisskalierte Variablen ist das beim BMI beschriebene Vorgehen anzuwenden.
- Kategoriale Variablen sollten im Vorfeld in Dummys überführt werden und die bis auf die gewünschte Referenzkategorie mit ins Modell aufgenommen werden. Die Interpretation ist analog zu den im nächsten Punkt beschriebenen Dummy-Variablen.
- Eine dichotome Variable (=Dummy) wird im Falle eines hinreichend kleinen p-Wertes (im Beispiel nicht gegeben) wie folgt interpretiert: der Koeffizient bezieht sich immer auf den Vergleich zwischen Testkategorie und Referenzkategorie. Die Referenzkategorie besitzt immer die niedrigere Ausprägung, hier: Männer.
Das bedeutet, wie verändert sich die y-Variable bei Vorhandensein der Testkategorie im Vergleich zur Referenzkategorie.
Im Beispiel würde dies bei einem hinreichend kleinen p-Wert bedeuten:
Testkategorie (Geschlecht = 1; Frau) vergleichen mit Referenzkategorie (Geschlecht = 0; Mann) hat einen negativen Koeffizienten. Demnach sinkt die Chance auf das Vorhandensein einer Krankheit für Frauen, im Vergleich mit Männern.
Aber ACHTUNG: der p-Wert ist zu groß und die Interpretation erfolgt hier nur zu Demonstrationszwecken.
3.3 Testmodell interpretieren - Odd's ratios
Den Odd's ratios (Chancenverhältnisse) kommt im Rahmen der logistischen Regression und der Ergebnisinterpretation besondere Aufmerksamkeit zu. Sie müssen berechnet werden, was mit der exp()-Funktion einfach umzusetzen ist. Zur besseren Übersicht habe ich mit der round()-Funktion eine Rundung auf 3 Dezimalstellen vorgenommen:
# Ungerundet
exp(coef(model1))
# Gerundet
round(exp(coef(model1)),3)
(Intercept) Geschlecht BMI
0.000 0.244 1.626
Im Falle der Notwendigkeit von Konfidenzintervallen für das Odd's ratio wird folgender Code verwendet:
# OR mit KI
round(exp(cbind(OR = coef(model1), confint(model1))),3)
Das führt zu folgendem Ergebnis:
OR 2.5 % 97.5 %
(Intercept) 0.000 0.000 0.007
Geschlecht 0.244 0.040 1.171
BMI 1.626 1.222 2.413
Die Odd's ratios sind nun in der ersten Spalte (OR) zu sehen, das Konfindenzintervall (KI) dahinter.
Hinweise zum OR:
- Ein OR > 1 zeigt eine höhere Chance an, dass die abhängige Variable = 1 ist.
- Ein OR < 1 zeigt eine niederigere Chance an, dass die abhängige Variable = 1 ist.
- Ist das OR = 1 ist die Chance weder höher noch niedriger, das die abhängige Variable = 1 ist.
- Das OR-Konfidenzintervall und der p-Wert kommen zu identischen Entscheidungen. Eine 1 im 95%-KI ist gleichbedeutend mit einem p-Wert > 0.05
Interpretation des OR:
- Da nur der BMI einen hinreichend kleinen p-Wert hat (siehe oben), ist nur das OR dieser Variable interpretierbar: Es ist 1.626. Steigt der BMI um 1 Einheit, steigt die relative Wahrscheinlichkeit, dass die Krankheit vorliegt um 62,6% (1.626 - 1 = 0.626).
- Bei Dummy-Variablen ist die Interpretation analog, kann jedoch auch noch etwas vereinfacht werden. Auch hier erneut der Hinweis: Im Beispiel ist der p-Wert nicht hinreichend klein und die Interpretation des OR erfolgt nur zu Demonstrationszwecken:
- Das OR vom Geschlecht ist mit 0.244 < 1. Demnach nimmt die Chance ab, dass die Krankheit vorliegt, wenn man die Testkategorie (Geschlecht = 1; Frau) mit der Referenzkategorie (Geschlecht = 0 ; Mann) vergleicht.
- Schematisch ausgedrückt: Frauen habe eine 0.244-mal höhere Chance zu erkranken als Männer.
- Zur verständlicheren Interpretation wird das Reziproke gebildet: 1/0.244 = 4.098
- Die Interpretation kann dann umgedreht werden: Die Chance, dass Männer die Krankheit haben ist 4.098-mal höher als bei Frauen - wenn der p-Wert hinreichend klein wäre..
4 Modellgüte berechnen und einordnen
Die Modellgüte kann anhand sog. Pseudo-Bestimmtheitsmaße abgelesen werden. Sie sind so konstruiert, dass sie analog zum R² aus der linearen Regression gelesen werden können.
Mit R2cs wird das sog. Cox and Snell-R² berechnet. Es kann allerdings maximal 0.75 erreichen, weswegen es stets niedriger als andere R² ist.
Mit R2n wird das Nagelkerke-R² berechnet. Es korrigiert das Cox and Snell-R² und kann wie das normale R² Werte zwischen 0 und 1 annehmen, weshalb ich es vorziehe das R² nach Nagelkerke zu interpretieren.
# Gütemaße
n <- length(model1$residuals)
R2cs <- 1-exp((model1$deviance-model1$null.deviance)/n)
R2n <- R2cs/(1-exp(-(model1$null.deviance/n)))
R2n
Im Ergebnis erhält man ein R² nach Nagelkerke in Höhe von 0.4716974.
Pauschal kann eine Einordnung des R² - ähnlich zur linearen Regression - nicht vorgenommen werden. Höher ist natürlich besser, allerdings sollte hier ein Vergleich mit Modellen vorgenommen werden, die die gleiche abhängige Variable untersuchen.
5 Schnelle Berechnung mit dem rms-Paket
Möchte man schnell die Ergebnisse erhalten, kann man auch das rms-Paket verwenden.
Die Ergebnisse und die Interpretation ist analog zu oben.
library(rms)
rmsmodel <- lrm(Krankheit~Geschlecht+BMI, data=df, x=TRUE, y=TRUE)
rmsmodel
#OR
round(exp(coef(rmsmodel)),3)
Der Output ist folgender:
Logistic Regression Model
lrm(formula = Krankheit ~ Geschlecht + BMI, data = df, x = TRUE,
y = TRUE)
Model Likelihood Discrimination Rank Discrim.
Ratio Test Indexes Indexes
Obs 51 LR chi2 20.13 R2 0.472 C 0.839
0 37 d.f. 2 R2(2,51) 0.299 Dxy 0.678
1 14 Pr(> chi2) <0.0001 R2(2,30.5)0.448 gamma 0.679
max |deriv| 6e-08 Brier 0.117 tau-a 0.275
Coef S.E. Wald Z Pr(>|Z|)
Intercept -11.5638 3.9561 -2.92 0.0035
Geschlecht -1.4113 0.8375 -1.69 0.0919
BMI 0.4863 0.1700 2.86 0.0042
>
> #OR
> round(exp(coef(rmsmodel)),3)
Intercept Geschlecht BMI
0.000 0.244 1.626
6 Videotutorial
