Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.

Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an .

Versionshistory:

  • v1.0: erste online-gestellte Version (03.3.2024)

9 Logistische Regression

In den Teilen 7-8 haben wir uns mit Varianten der (multiplen) Regression beschäftigt. Während wir für die Prädiktoren bereits verschiedene Skalenniveaus zugelassen haben, sind wir immer von einem (mindestens) intervallskalierten Kriterium ausgegangen. In diesem Teil 9 behandeln wir nun die logistische Regression, bei der ein binäres Kriterium vorhergesagt wird; also eines, welches nur die Werte 0 und 1 annimmt. Dies trifft auf viele Variablen zu, die potenziell erhoben werden und bei denen “ja” vs. “nein” gemessen wird (z.B. Krankheit ja vs. nein, Gabelstaplerschein ja vs. nein, …). Auch auf kognitionspsychologische Experimenten trifft dies zu, wenn z.B. fehlerhafte vs. korrekte Durchgänge unterschieden werden. Weitergehende Regressionsmodelle werden wir gegen Ende dieses Teils als Ausblick nur knapp behandeln.

# Pakete die in diesem Teil benutzt werden:
library(car)

9.1 Grundlagen

9.1.1 Warum logistische Regression? Das “Verallgemeinerte Lineare Modell”

Wir gehen nun immer davon aus, dass die binäre abhängige Variable nur die Werte 0 und 1 annehmen kann. In einer (fiktiven) Studie wurde nun einerseits das Alter von Versuchspersonen erhoben und andererseits, ob die Versuchspersonen verheiratet sind oder nicht. Zur Vorhersage des letzteren Status soll nun das Alter benutzt werden: “Alter” ist also der Prädiktor und “Verheiratet” ist das Kriterium, welches aber nur die Werte 0 (“nein”) und 1 (“ja”) annehmen kann. Die Ergebnisse sind in der folgenden Abbildung eingezeichnet. Mit diesen Daten können wir natürlich eine normale einfache lineare Regression berechnen und wir erhalten dann auch Koeffizienten und können eine Regressionsgerade in die Abbildung einzeichnen. Dies ist die blaue Linie in der folgenden Abbildung:

modell <- lm(verheiratet ~ alter, 
             data = daten)          # Regression berechnen
coef(modell)                        # Koeffizienten ausgeben
## (Intercept)       alter 
## -1.63795510  0.07144378

Hierbei wird ein Problem allerdings sehr deutlich: Benutzt man zur Vorhersage ein normales Regressionsmodell, so können Werte vorhergesagt werden, die außerhalb des zulässigen Wertebereichs (hier also 0 und 1) liegen.

Nun ist das Ziel der logistischen Regression aber auch, eine Funktion zu finden, die die Daten möglichst gut beschreibt: Anstatt einer Geraden wie bei der einfachen linearen Regression wird hier die sog. logistische Funktion benutzt (die die Verteilungsfunktion der logistischen Verteilung ist). Die logistische Funktion hat den Vorteil, dass sie nur Werte zwischen 0 und 1 annimmt, also Werte im zulässigen Wertebereich. In der folgenden Abbildung sind nochmals die Datenpunkte von gerade eingezeichnet, allerdings ist die blaue Linie die logistische Funktion, mit denjenigen Parametern, die sich aus der logistischen Regression ergeben:

Zusätzlich zu diesem Problem der Nicht-Linearität gibt es noch zwei weitere Punkte, weswegen bei einem binären Kriterium eine lineare Regression nicht adäquat ist. Insbesondere um Nullhypothesen-Signifikanztests über die Parameter durchführen zu können, haben wir nämlich eine Normalverteiltheit der Residuen und Homogenität der bedingten Varianzen der Residuen vorausgesetzt. Beide Punkte sind bei binärem Kriterium nicht erfüllt.

Die Normalverteiltheit von Residuen setzt eine kontinuierliche Variable voraus. Binäre Variablen können aber nur zwei Werte annehmen und und damit können auch die bedingten Residuen nur zwei Werte annehmen und somit kann die Verteilungsvoraussetzung hier nicht erfüllt werden. Auch Homogenität der bedingten Varianzen kann nicht gelten. Dies wird intuitiv klar, wenn man bedenkt, wie die Varianz einer binären Variablen berechnet wird (vgl. Teil 4, Statistik I): \[S_X^2=P\cdot(1-P)\] Dies bedeutet, dass die Varianz eine Funktion von \(P\) ist. Sie erreicht ihr Maximum \(0.25\) bei \(P=0.5\) und wird ansonsten kleiner:

Im Rahmen der logistischen Regression werden diese Probleme dadurch gelöst, dass die “Zielfunktion” der logistischen Funktion durch eine sog. Link-Funktion linearisiert wird und das Ergebnis wiederum durch eine Linearkombination von Prädiktoren modelliert werden kann. Insofern kann das Verallgemeinerte Lineare Modell mit drei Komponenten beschrieben werden als: \[E(\mathbf{Y})=g(\mu)=a+b_1X_1+\ldots+b_qX_q\]

  1. Eine systematische Komponente (ganz rechts) als Linearkombination von Prädiktoren, ganz ähnlich wie bei der (multiplen) linearen Regression.
  2. Eine zufällige Komponente (ganz links), die das Kriterium und dessen Verteilung spezifiziert. Auch dies ist analog zum linken Teil der linearen Regression, wo eben eine Normalverteilung des Kriteriums angenommen wird. Im Fall eines binären Kriteriums wird hier typischerweise eine Binomialverteilung angenommen; bei Häufigkeiten als Kriterium wird üblicherweise eine Poisson-Verteilung angenommen.
  3. Die Link-Funktion \(g(\mu)\) spezifiziert schließlich die Beziehung zwischen der zufälligen und der systematischen Komponente. Dadurch wird außerdem erreicht, dass die zufällige Komponente nicht mehr normalverteilt sein muss. Nimmt man für die zufällige Komponente eine stetige Variable mit Normalverteilung an und setzt als Link-Funktion die sog. Identity-Funktion \(g(\mu)=\mu\), dann resultiert nichts anderes als das Allgemeine Lineare Modell (ALM). Das ALM ist also ein Spezialfall des Verallgemeinerten Linearen Modells. Für eine normalverteilte Variable ist die Identity-Funktion die sog. kanonische Link-Funktion, was in etwa als “Standard-Link-Funktion” aufgefasst werden kann (es gibt dafür auch sehr technische Begründungen, die wir uns hier sparen). Im Fall einer Binomial-Verteilung ist die Inverse der logistischen Funktion (die sog. Logit-Transformation) die kanonische Link-Funktion und für eine Poisson-Verteilung ist es eine logarithmische Funktion.

Der Weg zum logistischen Regressionsmodell für den Fall eines binären Kriteriums wird nun in den nächsten drei Abschnitten beschrieben.

9.1.2 Die logistische Funktion

Die zu modellierende logistische Funktion hängt mit der logistischen Verteilung zusammen, die ihrerseits eine Dichtefunktion stetiger Zufallsvariablen ist. Man sagt, eine stetige Zufallsvariable sei logistisch verteilt mit den Parametern \(a\) und \(b\) (wobei \(b>0\) sei), wenn ihre Dichte gegeben ist durch die folgende Funktion: \[ f(x)=\frac{e^{-\frac{x-a}{b}}}{b\cdot \left(1+ e^{-\frac{x-a}{b}} \right)^2} \] Man bezeichnet \(a\) auch als Lageparameter und \(b\) als Skalierungsparameter (\(e\) ist die Euler’sche Zahl \(e\approx 2.71828\)). Der linke Teil der nächsten Abbildung zeigt exemplarisch drei logistische Verteilungen mit verschiedenen Parametern. Wie man sieht, ist eine gewisse Ähnlichkeit zur Normalverteilung nicht zu verleugnen. Im rechten Teil sind die dazugehörigen Verteilungsfunktionen dargestellt, die ihrerseits als sog. logistische Funktionen bezeichnet werden und beschrieben werden können als \[ F(x)=\frac{e^{\frac{x-a}{b}}}{1+ e^{\frac{x-a}{b}} } = \frac{1}{1+ e^{-\frac{x-a}{b}} } \]

Nachdem wir nun wissen, welche Zielfunktion wir modellieren wollen, gehen wir den nächsten Schritt und betrachten sog. Odds und die Logit-Transformation und deren Auswirkungen auf Werte, die bei einer logistischen Regression auftreten können.

9.1.3 Odds und Logits

Eines der Hauptprobleme im Umgang mit einem binären Kriterium war der begrenzte Wertebereich. Wir haben oben die Werte 0 und 1 verwendet und gesehen, dass eine normale lineare Regression unplausible Wertebereiche produziert. Statt einer Geraden hatten wir dann gesagt, es soll eine logistische Funktion resultieren. Um nun mit einer Linearkombination eine logistische Funktion modellieren zu können, müssen die Werte des Kriteriums so transformiert werden, dass aus der logistischen Funktion eine Gerade wird. Dazu wenden wir zwei Transformationen an, die dafür sorgen, dass der Wertebereich von \(-\infty\) bis \(+\infty\) geht.

Als Odds (oder auf Deutsch manchmal “(Wett-)Chance”) bezeichnet man das Verhältnis einer Wahrscheinlichkeit \(P\) zu ihrer Gegenwahrscheinlichkeit \(1-P\), also: \[O=\frac{P}{1-P}\] Wenn also ein Ereignis \(A\) eine Wahrscheinlichkeit von \(p(A)=0.75\) hat, dann ist die Chance für das Eintreten von \(A\) also \(O_A=\frac{0.75}{0.25}\) bzw. 3 zu 1. Umgekehrt kann aus Odds auch die Wahrscheinlichkeit berechnet werden: \[P=\frac{O}{1+O}\] Odds tragen also Informationen über Wahrscheinlichkeiten in sich, haben aber einen weiteren Vorteil: Sie haben zwar 0 als untere Grenze, nach oben geht der Wertebereich aber bis \(+\infty\). Dies wird im linken Teil der folgenden Abbildung deutlich:

Die Odds werden nun als nächstes in sog. Logits transformiert, indem der natürliche Logarithmus auf sie angewendet wird. Logits sind also definiert als: \[\text{Logit} = \ln (O) = \ln \left( \frac{P}{1-P}\right)\] Dadurch wird bewirkt, dass der Wertebereich auch seine untere Grenze verliert und nun also von \(-\infty\) bis \(+\infty\) geht (siehe rechter Teil der letzten Abbildung).

Nun schauen wir noch einmal die logistische Funktion für den Spezialfall \(a=0\) und \(b=1\) an: \[ F(X)=\frac{e^{\frac{x-a}{b}}}{1+ e^{\frac{x-a}{b}} } = \frac{e^x}{1+e^x} \] Zwischen logistischer Funktion und Logit-Transformation besteht nämlich eine sehr hilfreiche Beziehung: Die Logit-Transformation ist die Umkehrfunktion der logistischen Funktion. In anderen Worten: \[ \text{Logit}\left( \frac{e^x}{1+e^x} \right) = x \]

Das können wir sehr einfach einmal in R überprüfen:

x <- 0.6  # als Beispiel
( Fx <- (exp(x)) / (1+exp(x)) )    # Wert von x auf der logistischen Funktion
## [1] 0.6456563
( logit.Fx <- log( Fx / (1-Fx)) )  # Logit berechnen... kommt 0.6 heraus
## [1] 0.6
# Alternativ: 
( logit(Fx) )                      # Funktion logit() aus dem Paket car
## [1] 0.6

Das bedeutet aber auch: Wenden wir auf die Werte einer logistischen Funktion die Logit-Transformation an, so resultiert eine Gerade – und genau diese Gerade können wir mit einer Linearkombination von Prädiktoren modellieren. Dass tatsächlich eine Gerade resultiert wird im folgenden R-Code noch einmal visualisiert. Statt der Formel für die logistische Regression greifen wir direkt auf die Verteilungsfunktion \(F(x)\) der logistischen Verteilung logis in R zurück, genau wie wir es von anderen Verteilungen schon kennen. In der linken Abbildung ist für einen Wertebereich von -10 bis 10 die logistische Funktion abgebildet; die rechte Abbildung stellt die resultierenden Werte Logit-transformiert – wie erwartet als Gerade – dar:

par(mfrow = c(1, 2))
x <- seq(-10, 10, 0.01)     # Sequenz von x-Werten schaffen

# linke Abbildung: 
x.logis <- plogis(x)        # Verteilungsfunktion der logistischen Verteilung
plot(x, 
     x.logis, 
     type = "l",
     ylab = expression(italic(F)~plain("(")~italic(x)~")"))

#rechte Abbildung
logit.x.logis <- logit(x.logis)  # Logit der logistischen Verteilung
plot(x, logit.x.logis, type = "l", 
     ylab = expression(plain("Logit (")~italic(F)~plain("(")~italic(x)~"))"))

9.1.4 Die Regressionsgleichung(en) der logistischen Regression mit einem Prädiktor und Maximum-Likelihood Schätzung

Mit diesem Vorwissen können wir nun Gleichungen für die logistische Regression aufstellen. Hierbei sind drei verschiedene Varianten möglich, die allesamt ihre Vor- und Nachteile haben und ein wenig ist es auch fachspezifisch, wie die Regressionsgleichung bevorzugt formuliert wird: bezogen auf (1) Logits, (2) Odds oder (3) Wahrscheinlichkeiten. Wir berücksichtigen hier erst einmal nur den Fall der Regression mit einem Prädiktor. \(P\) meint dabei jeweils die Wahrscheinlichkeit derjenigen Kategorie anzugehören, die mit 1 kodiert wird, wenn ein bestimmter Wert auf dem Prädiktor \(X\) vorliegt.

  1. Beziehen wir die Regressionsgleichung auf Logits können wir nun einfach – fast wie eine normale lineare Regression – schreiben: \[\text{Logit}(P)=\ln \left(\frac{P}{1-P} \right) = a+bX\] Die Linearkombination des Prädiktors sagt also keine Wahrscheinlichkeiten \(P\) vorher, sondern nur deren Logit-Wert. Wir sehen aber, dass durch die Transformation im Prinzip das Grundgerüst der bereits bekannten Regressionen resultiert.
  2. Einfacher zu interpretieren sind daher i.d.R. eine der zwei weiteren Varianten. Beziehen wir die Regressionsgleichung auf die Odds, indem wir beide Seiten exponieren, dann resultiert: \[O = \frac{P}{1-P} = e^{a+bX}\]
  3. Schließlich können wir noch einen Schritt weitergehen und die Regressionsgleichung direkt auf Wahrscheinlichkeiten beziehen: \[P=\frac{O}{1+O} = \frac{e^{a+bX}}{1+e^{a+bX}} = \frac{1}{1+e^{-(a+bX)}} \]

Wie kommen wir nun zu den gesuchten Werten für \(a\) und \(b\) im hier beschriebenen Fall? Für den Spezialfall der Normalverteilung und der Identity-Link-Funktion können wir – genau wie wir es bisher getan haben – die Methode der kleinsten Quadrate verwenden. Das ist ansonsten und auch hier leider nicht mehr möglich, sondern es kommen sog. iterative Maximum-Likelihood Schätzungen zum Einsatz (siehe Teil 6 von Statistik II im Allgemeinen und diese Ergänzung für eine Ausformulierung am Beispiel der logistischen Regression).

9.2 Durchführung mit R und Interpretation

9.2.1 Berechnung der Regression

Die praktische Berechnung einer logistischen Regression mit R ist sehr einfach. Im Wesentlichen verwenden wir statt lm() nun die Funktion glm() und geben noch die Verteilungsannahmen explizit an. Das passiert für unseren Fall mit dem Argument family = binomial(link = "logit") (wobei wir link = "logit") auch weglassen können, weil es die Voreinstellung ist:

glm.ergebnis <- glm(verheiratet ~ alter, 
                    data = daten,
                    family = binomial(link = "logit"))

Auch viele der Ausgabefunktionen sind nun genau wie bei lm(). Zunächst können wir auch Werte für die Koeeffizienten, also für \(a\) und \(b\), ausgeben:

coef(glm.ergebnis)
## (Intercept)       alter 
## -21.3794357   0.7140123

Eine genauere Zusammenfassung erhalten wir mit summary() oder S():

# summary(glm.ergebnis)               # R Standard
S(glm.ergebnis, brief = TRUE)         # aus car
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -21.3794     7.7154  -2.771  0.00559 **
## alter         0.7140     0.2564   2.785  0.00535 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.452  on 39  degrees of freedom
## Residual deviance: 25.105  on 38  degrees of freedom
## 
## logLik     df    AIC    BIC 
## -12.55      2  29.10  32.48 
## 
## Number of Fisher Scoring iterations: 6

Während einiges hier bekannt vorkommen mag, ist doch die Interpretation der Koeffizienten anders als gewohnt (siehe nächster Abschnitt). Daneben gibt es auch Informationen, die wir bisher nicht kennen (z.B. Null deviance und Residual deviance), die wir im übernächsten Abschnitt ansprechen werden. Zusätzlich wird unter Number of Fisher Scoring iterations angegeben, wie viele Iterationen benötigt wurden, um das Modell zu schätzen. Eine hohe Zahl (manchmal wird gesagt \(>25\)) deutet darauf hin, dass das Modell nicht gut mit den Daten zusammenpasst.

9.2.2 Interpretation der Ergebnisse

Im Prinzip haben die Koeffizienten ganz ähnliche Interpretationen wie in der (einfachen) linearen Regression. Das Problem entsteht allerdings, wenn man bedenkt, welche Größe hier vorhergesagt wird: nämlich Logits. Im Prinzip ist aber \(a=-21.38\) der Wert in Logits, der vorhergesagt wird, wenn alle Prädiktoren (hier also nur einer) den Wert 0 annehmen. Der Koeffizient \(b=0.714\) gibt an, um welchen Wert sich der Logit erhöht, wenn der Prädiktor um den Wert 1 erhöht wird. Soweit so gut, aber klar ist auch, dass dies eine sehr gewöhnungsbedürftige Interpretation ist, da das Denken in logarithmierten Odds nicht sehr angenehm ist.

Daher betrachten wir das Ergebnis nun anders, nämlich wenn wir die Regressionsgleichung auf Odds beziehen, also beide Seiten der Regressionsgleichung auf Logits bezogen exponieren und direkt eine kleine Umformung durchführen: \[O=e^{a+bX}=e^a \cdot e^{bX}\] Setzen wir nun zunächst \(X=0\), dann gibt \(e^a\) (statt \(a\)!) die Odds eines Erfolges (also, dass das Kriterium den Wert 1 annimmt) an, wenn alle Prädiktoren auf Null gesetzt sind. Den entsprechenden Wert können wir in R ganz einfach mit der Funktion exp() bestimmen:

exp(coef(glm.ergebnis))
##  (Intercept)        alter 
## 5.188347e-10 2.042169e+00

Wie man sehen kann, ist der Wert sehr nahe an Null. Wenn \(x=0\) eingesetzt wird, sind die Odds, also die Chance, verheiratet zu sein, quasi Null.

Die Interpretation von \(e^b\) ist leider noch etwas komplizierter, wird aber aus der Umformung von gerade eben auch deutlich: Wenn \(x=0\) ist, dann ist \(O=e^a\). Wenn nun die Variable \(X\) um einen Schritt erhöht wird und \(x=1\) ist, dann muss also \(e^a\) mit \(e^b\) multipliziert werden und das Ergebnis ist dann wieder ein Odds, also die Chance verheiratet zu sein wenn \(x=1\) ist.

Wir machen uns das noch einmal am Beispiel \(x=0\) und \(x=1\) deutlich. Nachdem was wir gerade gesagt haben, gilt also: \[ \frac{P(Y=1|x=0)}{1-P(Y=1|x=0)}=e^a \quad\text{und}\quad \frac{P(Y=1|x=1)}{1-P(Y=1|x=1)}=e^a\cdot e^b \] Nun setzen wir \(e^a\) in den rechten Teil ein und erhalten: \[\begin{equation*} \begin{aligned} &\frac{P(Y=1|x=1)}{1-P(Y=1|x=1)}=\frac{P(Y=1|x=0)}{1-P(Y=1|x=0)}\cdot e^b\\ \Leftrightarrow & \frac{\frac{P(Y=1|x=1)}{1-P(Y=1|x=1)}}{\frac{P(Y=1|x=0)}{1-P(Y=1|x=0)}}=e^b \end{aligned} \end{equation*}\]

Zusammengefasst ist also \(e^b\) (nicht \(b\)!) der Quotient zweier Odds. Dies wiederum wird als Odds Ratio bezeichnet und man könnte auch sagen: \[ e^b=\frac{O_{\text{nach Erhöhung um 1}}}{O_{\text{vor Erhöhung um 1}}} \] Der Koeffizient \(e^b\) gibt also die Veränderung der Odds an, wenn der Wert von \(X\) um 1 erhöht wird.

Schließlich können wir noch die Regressionsgleichung bezogen auf Wahrscheinlichkeiten anschauen und diese auch gleich minimal umformen: \[P=\frac{e^{a+bX}}{1+e^{a+bX}}=\frac{e^a\cdot e^{bX}}{1+(e^a\cdot e^{bX})}\] Setzen wir wieder \(x=0\) bleibt also \[P=\frac{e^a}{1+e^a}\] und damit können wir die Wahrscheinlichkeit eines Erfolges (also, dass das Kriterium den Wert 1 annimmt) berechnen, wenn alle Prädiktoren den Wert 0 annehmen. Mit der vorherigen Formel können wir auch die respektiven Wahrscheinlichkeiten für andere Werte von \(X\) berechnen.

Die folgende Abbildung stellt Logits, Odds und Wahrscheinlichkeiten für unser gerade berechnetes Modell grafisch dar, indem Logits, Odds und Wahrscheinlichkeiten für verschiedene Werte von \(X\) über die Formeln berechnet werden. Hier wird noch einmal klar, dass für Logits eine lineare Funktion vorhergesagt wird, während für Wahrscheinlichkeiten die logistische Funktion resultiert:

a <- coef(glm.ergebnis)[1]  # a extrahieren
b <- coef(glm.ergebnis)[2]  # b extrahieren
par(mfrow = c(1, 3)) # 3 Plots nebeneinander

X <- seq(0:50)    # Vektor mit X-Werten

## Logits
logits <- a + b * X
plot(X,
     logits, 
     type = "l",
     main = "Logits")

## Odds
odds <- exp(a + b * X)
plot(X, 
     odds, 
     type = "l",
     main = "Odds")

## Wahrscheinlichkeit P
P <- exp(a + b * X) / ( 1 + exp(a + b * X))
plot(X, 
     P,
     type = "l",
     main = "P")

9.2.3 Modellgüte und Inferenzstatistik

In der linearen Regression (als ein Beispiel des Allgemeinen Lineares Modells) hatten wir als Maß der Modellgüte \(R^2\) verwendet und das Gesamtmodell mit einem \(F\)-Test gegen die Nullhypothese \(H_0:R^2=0\) getestet. Einzelne Koeffizienten wurden wiederum mit \(t\)-Tests auf Signifikanz getestet und verschiedene Modelle konnten wir auch mit einem \(F\)-Test vergleichen. Im Verallgemeinerten Linearen Modell, und damit auch für die logistische Regression, ist die Vorgehensweise etwas anders. Wir unterscheiden auch wieder Tests über das ganze Modell und Tests einzelner Parameter.

9.2.3.1 Modellgüte und Modellvergleich: der Likelihood-Ratio-Test

\(R^2\) hatten wir auch beschrieben als Korrelation der beobachteten und der vorhergesagten Werte. Die Maximum-Likelihood Methode zur Schätzung der Parameter wirft als “Nebenprodukt” die (maximierte) Log-Likelihood \(\ell\) ab, also der logarithmierte Maximalwert der Likelihood (der typischerweise erreicht wird, wenn der Schätz-Algorithmus konvergiert). Aus der Log-Likelihood \(\ell\) wiederum kann sehr einfach die sog. Deviance-Statistik errechnet werden als (vgl. Teil 6 von Statistik II): \[\text{Deviance}=-2\cdot\ell\]

In gewisser Weise ist die Deviance-Statistik analog zur Quadratsumme der Residuen, da auch sie damit zu tun hat, wie viel Information in den Daten unerklärt bleibt, wenn die “besten” Koeffizienten gefunden wurden. Je höher der Wert ist, desto schlechter kann das Modell die Daten erklären.

Diese Statistik wiederum hat approximativ, also annähernd, eine \(\chi^2\)-Verteilung. Differenzen der Deviance-Statistik für verschiedene aufeinander aufbauende Modelle können dann an einer \(\chi^2\)-Verteilung mit \(m\) Freiheitsgraden getestet werden, wobei \(m\) die Differenz in der Anzahl der Parameter ist, die beide Modelle haben. Dies wird oft als Likelihood-Ratio-Test (LRT) bezeichnet und ist quasi analog zum Vergleich verschiedener linearer Modelle mit dem \(F\)-Test (vgl. Teil 7 von Statistik II). Die Ausgabe von summary(glm.ergebnis) beinhaltet zwei Zeilen mit Werten für Deviance-Statistiken:

  1. Die Zeile Null deviance: 55.452 ist die Deviance-Statistik für das sog. Nullmodell, welches nur den Parameter \(a\) und keine weiteren Prädiktoren beinhaltet.
  2. Die Zeile Residual deviance: 25.105 ist die Deviance-Statistik für das Modell mit alter als zusätz-lichem Prädiktor.

Zunächst sollte der Wert für die Residual deviance kleiner sein als der für die Null deviance. Anders ausgedrückt: Je größer die Differenz beider Deviance-Statistiken ist, umso besser kann das Modell mit Prädiktor die Daten erklären. Die Differenz kann nun mit einem LRT auf Signifikanz (i.S.v.: die Differenz ist größer als Null) getestet werden. Dazu berechnen wir zusätzlich noch das Nullmodell, indem wir nur den Parameter \(a\) zulassen. Dies wird getan, indem rechts von der Tilde ~ explizit eine 1 steht:

glm.ergebnis.null <- glm(verheiratet ~ 1, 
                         data = daten, 
                         family = "binomial")

Den eigentlichen LRT führen wir mit der Funktion anova() durch, indem wir die zu vergleichenden Modelle hintereinander übergeben und zusätzlich den zu verwendenden \(\chi^2\)-Test spezifizieren:

anova(glm.ergebnis.null, glm.ergebnis, 
      test = "Chisq")

Die Differenz der Deviance-Statistiken beträgt also 30.35 und getestet an einer \(\chi^2\)-Verteilung mit 1 Freiheitsgrad ergibt sich \(p<.001\). Die Hinzunahme des Prädiktors alter verbessert das Modell also signifikant (bei \(\alpha=0.05\)).

9.2.3.2 Wald-Test und Konfidenzintervall einzelner Parameter

Die Koeffizienten in der (multiplen) linearen Regression wurden mit \(t\)-Tests gegen die Nullhypothese \(H_0:\beta=0\) getestet. Im hier betrachteten Fall greifen wir auf den Wald-Test (benannt nach Abraham Wald) zurück: Der Quotient eines quadrierten Koeffizienten und dessen quadrierten Standardfehlers ist approximativ \(\chi^2\)-verteilt mit 1 Freiheitsgrad: \[ \frac{b^2}{\hat{\sigma}^2_b}\approx \chi^2_1 \] Dieser Test hat den Vorteil, dass er leicht auf \(J\) Parameter gleichzeitig ausgedehnt werden kann: \[ \sum_{j=1}^J\frac{b^2_j}{\hat{\sigma}^2_{b_j}}\approx \chi^2_J \] Für den Fall eines einzelnen Parameters kann eine Ableitung aus diesen Fällen benutzt werden, bei der auf die Normalverteilung zurückgegriffen werden kann: \[ \frac{b}{\hat{\sigma}_b}\approx z \sim N(0,1) \] Dies ist auch der Test, den summary() im Fall der logistischen Regression ausgibt. Konfidenzintervalle für die Koeffizienten können leicht mit der Funktion confint() angefordert werden, die auf das glm()-Objekt angewendet wird:

confint(glm.ergebnis)
## Es wird auf das Profilieren gewartet ...
##                   2.5 %    97.5 %
## (Intercept) -41.8631447 -10.14622
## alter         0.3404704   1.39368

Ebenso können wir die exponierten Koeffizienten, die ja eine leichtere Interpretation als Odds Ratio erlauben, und die dazugehörigen Konfidenzintervalle anfordern:

exp(coef(glm.ergebnis))     # exponierte Koeffizienten...
##  (Intercept)        alter 
## 5.188347e-10 2.042169e+00
exp(confint(glm.ergebnis))  # ... und deren Konfidenzintervalle
## Es wird auf das Profilieren gewartet ...
##                    2.5 %       97.5 %
## (Intercept) 6.592760e-19 3.922398e-05
## alter       1.405609e+00 4.029652e+00

Diese Werte werden üblicherweise in Ergebnisteilen mitberichtet. Wenn wir statt summary() die Funktion S() benutzen, werden diese Werte direkt mit ausgeben (ganz am Ende):

S(glm.ergebnis)
## Call: glm(formula = verheiratet ~ alter, family = binomial(link = "logit"),
##           data = daten)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept) -21.3794     7.7154  -2.771  0.00559 **
## alter         0.7140     0.2564   2.785  0.00535 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 55.452  on 39  degrees of freedom
## Residual deviance: 25.105  on 38  degrees of freedom
## 
## logLik     df    AIC    BIC 
## -12.55      2  29.10  32.48 
## 
## Number of Fisher Scoring iterations: 6
## 
## Exponentiated Coefficients and Confidence Bounds
##                 Estimate        2.5 %       97.5 %
## (Intercept) 5.188347e-10 6.592760e-19 3.922398e-05
## alter       2.042169e+00 1.405609e+00 4.029652e+00

9.3 Ausblick: Weiterführende Regressionsmodelle

Wir haben hier nur einen einzigen Fall des Verallgemeinerten Linearen Modells betrachtet, nämlich die logistische Regression mit einem Prädiktor. Abschließend erwähnen wir hier weitere Varianten, ohne im Detail auf diese einzugehen.

Zunächst einmal kann eine logistische Regression mit mehreren Prädiktoren durchgeführt werden, ganz ähnlich wie es bei der multiplen linearen Regression der Fall ist. Die Syntax der Funktion glm() bleibt dabei die gleiche, nur werden rechts der Tilde ~ mehrere Prädiktoren angegeben, wie hier z.B.:

glm.ergebnis.multiple <- glm(kriterium ~ praed1 + praed2,
                             data = daten, 
                             family = binomial())

Auch die meisten Aspekte der Interpretation ähneln dann denen der multiplen linearen Regression: Jeder Koeffizient gibt nun eine Veränderung (des Logits, der Odds oder der Wahrscheinlichkeit) an, wenn die anderen Prädiktoren konstant gehalten werden. Auch die Durchführung von LRTs kann leicht verallgemeinert werden, um verschiedene Modelle miteinander zu vergleichen.

Die logistische Regression kann auch erweitert werden auf den Fall nominalskalierter Variablen mit mehr als zwei Ausprägungen als Kriterium. Dann spricht man von einer Multinomialen Regression (von der die logistische Regression ein Spezialfall mit zwei Ausprägungen ist). Die Annahme hierbei ist, dass das Kriterium als multinomial-verteilt angenommen wird (die Binomial-Verteilung ist dann wieder ein Spezialfall). Dann werden die \(J\) Kategorien durch \(J-1\) Variablen dargestellt (ganz ähnlich wie bei der Dummykodierung) und es werden \(J-1\) logistische Regressionen simultan gelöst, d.h. die Likelihood wird über alle Gleichungen simultan maximiert, um die Koeffizienten zu bestimmen.

Häufigkeiten als Kriterium stellen eine weitere Besonderheit dar, da sie eine untere Grenze haben (nämlich 0). Im Wesentlichen wird hier angenommen, dass Häufigkeiten einer Poisson-Verteilung folgen und als kanonische Link-Funktion wird der Logarithmus verwendet. Daher redet man dann auch von einer Poisson-Regression

Schließlich sind Hierarchische Lineare Modelle (HLM) (oder Multilevel Modelle oder Mixed (Effect) Modelle) noch eine wichtige Erweiterung der Regressionsrechnung. Ausgangslage hierbei ist, dass Beobachtungen (also z.B. Versuchspersonen) in übergeordneten Ebenen strukturiert sind. Ein typisches Beispiel (siehe nächste Abbildung) hierfür sind individuelle Schüler:innen (Level 1), die verschiedenen Klassen angehören (Level 2), die sogar wiederum in verschiedenen Schulen klassifiziert sind (Level 3). Aber auch einzelne Reaktionszeiten als Level 1-Beobachtungen gehören zu bestimmten Level 2-Versuchspersonen und können mit derartigen Modellen ausgewertet werden. Hierfür wird das Regressionsmodell noch einmal erweitert und es wird i.W. die Funktion lmer() zur Auswertung verwendet. Eine Einführung findet sich z.B: unter https://rpubs.com/favstats/multilevel.