Für eine optimale technische Stabilität empfehlen wir, dieses Online-Lehrbuch auf einem Notebook oder Desktop-Computer zu nutzen. Die interaktiven Komponente, insbesondere die R-Konsole (WebR), sind rechenintensiv und funktionieren auf mobilen Geräten nur eingeschränkt.
Viele der in diesem Companion behandelten Methoden basiert auf dem Konzept der Schätzung kausaler Effekte mit Regression. Als Regressionsansatz bezeichnet man eine Methode, welche die Beziehungen zwischen Variablen durch einen funktionalen Zusammenhang beschreibt und die Parameter der gewählten funktionalen Form anhand von beobachteten Daten schätzt. Lineare Regression nimmt eine lineare funktionale Form der Beziehung zwischen einer abhängigen Variable (Outcome-Variable) und erklärenden Variablen (Regressoren) an. Nicht-lineare Regressionsmethoden modellieren die Beziehung etwa durch Polynome höherer Ordnung, exponentielle Funktionen oder andere komplexere Formen. Die Wahl der funktionalen Form hängt von der Natur des datenerzeugenden Prozesses (DGP) und somit stets von der spezifischen Beziehung ab, die untersucht wird.
Regressionsansätze gehören zu den am häufigsten verwendeten Methoden für Kausalanalysen, da mit Regression in vielen Forschungsdesigns kausale Effekte identifiziert werden können, indem Backdoors geschlossen werden: Regression kann die durch die Behandlungsvariable verursachte Variation in der Outcome-Variable isolieren, indem für gemeinsame Einflussfaktoren von Behandlungs- und Outcome-Variable kontrolliert wird. In diesem Kapitel erläutern wir Erweiterungen der Spezifikation und Schätzung von Regressionsansätzen, die für spätere Kapitel relevant sind. Neben einer Motivation der Schätzung kausaler Effekte mit multipler Regression betrachten wir Modelle für verschiedene Kategorien von Outcome-Variablen und diskutieren deren Implementierung mit R.
Das Frisch-Waugh-Lovell-Theorem (FWL) besagt, dass die geschätzten Koeffizienten für eine interessierende Teilmenge der Regressoren in einer multiplen linearen Regression numerisch identisch mit den Koeffizientenschätzungen aus folgenden Schritten sind:
Rechne die Effekte der übrigen Regressoren auf (a) die Outcome-Variable und (b) die interessierende Teilmenge der Regressoren mit Regression heraus.
Regressiere anschließend die Residuen von Schritt (a) auf die Residuen aus Schritt (b).
In einem multiplen Modell mit zwei Regressoren \(X_1,\ X_2\), \[\begin{align}
Y = \beta_0 + \beta_1 X + \beta_2 X_2 + \epsilon \label{eq:fwlfullreg}
\end{align}\] kann der Effekt von \(X_1\) auf \(Y\) also mit der Regression \[\begin{align*}
\widehat{u}_{Y,X_2} = \beta_1 \widehat{u}_{X_1,X_2} + e \label{eq:fwl2reg}
\end{align*}\] geschätzt werden, wobei \(\widehat{u}_{Y,X_2}\) und \(\widehat{u}_{X_1,X_2}\) die Residuen der Regression von \(Y\) auf \(X_2\) und von \(X_1\) auf \(X_2\) sind.
FWL ermöglicht daher eine Vereinfachung der Schätzung komplexer Modelle durch die Zerlegung der Schätzung in Teilschritte.
Für das Verständnis der Schätzung kausaler Effekte mit linearer Regression ist FWL hilfreich, denn es zeigt, wie sowohl die Variation in der Outcome-Variable (\(\widehat{u}_{Y,X_2}\)) als auch die Variation in der Behandlungsvariable (\(\widehat{u}_{X_1,X_2}\)), die jeweils nicht durch Kovariablen (\(X_2\)) verursacht wird, mit multipler Regression isoliert werden kann, sodass Backdoors geschlossen werden.
Wir illustrieren dieses Konzept anhand einer multiplen Regression für das Gewicht (body_mass) von Pinguinen aus dem Datensatz palmerpenguins::penguins,
\[\begin{align}
\textup{body\_mass} = \beta_0 + \beta_1\cdot\textup{bill\_length} + \beta_2\cdot \textup{flipper\_length} + \epsilon,\label{eq:billdepthmodel}
\end{align}\] unter der Annahme, dass \(\beta_1\) der interessierende Effekt ist: Die erwartete Änderung des Gewichts eines Pinguins (in Gramm) für eine Änderung der Schnabel-Länge um 1mm.
Vor der Schätzung von Modell \(\eqref{eq:billdepthmodel}\) lesen wir den Datensatz ein und erstellen eine bereinigte Variante penguins_cleaned, analog zur Vorgehensweise in Kapitel 2.1.2.
Wir schätzen nun Modell \(\eqref{eq:billdepthmodel}\) mit lm() und erhalten eine Zusammenfassung der geschätzten Koeffizienten mit broom::tidy().
Das Ergebnis der Schätzung ist \(\widehat{\beta}_1\approx3.80\). Der nächste Code-Block berechnet die Residuen aus den Regressionen \[\begin{align*}
\textup{body\_mass} =&\, \alpha_0 + \alpha_1 \textup{flipper\_length} + u_{\textup{body\_mass},\,\textup{flipper\_length}},\\
\textup{bill\_length} =&\, \delta_0 + \delta_1 \textup{flipper\_length} + u_{\textup{bill\_length},\,\textup{flipper\_length}},
\end{align*}\]
und speichert diese in body_mass_res und bill_length_res.
Für den zweiten Schritt regressieren wir body_mass_res auf bill_length_res.
Der geschätzte Koeffizient aus der Regression der Residuen stimmt mit dem geschätzten Koeffizienten von bill_length aus der großen Regression \(\eqref{eq:billdepthmodel}\) überein.
Wir können den Effekt der Kontrolle für flipper_length visualisieren. Wir plotten hierzu:
Die originalen Datenpunkte für bill_length und body_mass1 gemeinsam mit der geschätzten Regressionslinie für das Modell \[ \textup{body\_mass} = \beta_0 + \beta_1\textup{bill\_length} + u \] (keine Kontrolle für flipper_length!)2.
Die um flipper_length bereinigten Datenpunkte und die zugehörige geschätzte Regressionslinie.
1 Für eine bessere Lesbarkeit der Grafik zentrieren wir beide Variablen um den jeweiligen Stichprobenmittelwert.
2 Der R-Befehl für diese Regression ist lm(I(body_mass - mean(body_mass)) ~ I(bill_length - mean(bill_length)) - 1, data = penguins_cleaned).
Wir erweitern p um die ursprünglichen Datenpunkte und die zugehörige Regressionslinie.
Der grafische Vergleich beider Vorgehensweisen zeigt den Effekt der Kontrolle für flipper_length: Die geschätzte (schwarze) Regressionslinie für die bereinigten Daten hat eine deutlich geringere Steigung als die anhand der ursprünglichen Daten geschätzte (lilane) Linie. Der Effekt von bill_length auf body_mass wird mit der einfachen Regression lm(body_mass ~ bill_length) vermutlich überschätzt, weil es andere Faktoren (wie flipper_length gibt, die mit bill_length und body_mass korrelieren. Kontrollieren für flipper_length in der multiplen Regression lm(body_mass ~ bill_length + flipper_length) schließt die Backdoor durch flipper_length. Die Konsequenz ist eine deutlich geringere Steigung der lilanen Regressionslinie.
4.2 Binäre Outcome-Variable
Eine binäre Variable, auch als dichotome Variable oder Indikator-Variable bezeichnet, ist eine Variable, die nur zwei Ausprägungen annehmen kann. Diese beiden Ausprägungen werden typischerweise durch die Werte 0 und 1 repräsentiert und dienen dazu, zwei verschiedene Zustände oder Kategorien zu unterscheiden. Formal kann eine binäre Variable \(B\) wie folgt definiert werden:
\[\begin{align}
B = \begin{cases}
1, & \text{Eigenschaft trifft zu,} \\
0, & \text{Eigenschaft trifft nicht zu.}
\end{cases}
\end{align}\]
Ein in späteren Kapiteln dieses Companions verwendeter binärer Regressor ist der Indikator für die Zuordnung von Beobachtungen zur Behandlungs- oder Kontrollgruppe (1 = Behandlungsgruppe, 0 = Kontrollgruppe).
Für viele ökonomische Forschungsfragen ist es hilfreich, eine binäre Outcome-Variable mit Regression zu modellieren. Hierzu gibt es verschiedene Ansätze, die wir nachfolgend zusammenfassen und ihre Anwendung mit R zeigen.
4.2.1 Das lineare Wahrscheinlichkeitsmodell
Das lineare Regressionsmodell
\[Y = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2} + \dots + \beta_k X_{k} + u\] mit einer binären abhängigen Variablen \(Y_i\in\{0,1\}\) wird als lineares Wahrscheinlichkeitsmodell bezeichnet. Wie üblich modellieren wir den Erwartungswert der abhängigen Variable gegeben der Regressoren \(X_1,\dots,X_k\) als lineare Funktion,
\[E(Y\vert X_1,X_2,\dots,X_k) = P(Y=1\vert X_1, X_2,\dots, X_3).\] Da \(Y\) eine binäre Variable ist, gilt hier
Das lineare Wahrscheinlichkeitsmodell beschreibt also die Wahrscheinlichkeit, dass \(Y=1\) als lineare Funktion der Regressoren: \(\beta_j\) misst die Änderung in der Wahrscheinlichkeit für das Ereignis \(Y_i=1\), unter der Bedingung, dass die anderen \(k-1\) Regressoren konstant gehalten werden. Genau wie bei multipler Regression mit einer kontinuierlichen abhängigen Variablen können die \(\beta_j\) mit der KQ-Methode geschätzt werden.
Aufgrund der Beschränktheit der \(Y_i\) auf \(\{0,1\}\) ist \(u_i\) heteroskedastisch. Folglich sollten Inferenzstatistiken mit robusten Standardfehlern berechnet werden. Weiterhin ist zu beachten, dass \(R^2\) in den meisten Anwendungen von linearen Wahrscheinlichkeitsmodellen keine hilfreiche Interpretation hat, da das geschätzte Modell die Daten nicht perfekt erklären kann, wenn die abhängige Variable binär, aber die Regressoren kontinuierlich verteilt sind.
Das lineare Wahrscheinlichkeitsmodell hat einen wesentlichen Nachteil: Das Modell nimmt an, dass die bedingte Wahrscheinlichkeitsfunktion linear ist und \(P(Y=1\vert X_1,\dots,X_k)\) über das für Wahrscheinlichkeiten definierte Intervall \([0,1]\) hinausgehen kann. Ein angepasstes Modell hat dann für Regressorwerte, die zu Vorhersagen von \(Y\) jenseits von \([0,1]\) führen keine sinnvolle Interpretation.
Dieser Umstand verlangt nach Regressionsansätzen, die \(P(Y=1)\) durch eine auf \([0,1]\) beschränkte (nicht-lineare) Funktion der Regressoren modellieren. Häufig verwendete Methoden sind Probit- und Logit-Regression.
Ein lineares Wahrscheinlichkeitsmodell kann mit lm() geschätzt werden, wobei die abhängige Variable den Typ numeric oder integer haben muss.
4.2.2 Probit-Regression
Bei der Probit-Regression wird die Standardnormalverteilungsfunktion \(\Phi(\cdot)\) verwendet, um die Regressionsfunktion einer binären abhängigen Variable zu modellieren. Wir nehmen an, dass \[\begin{align}
E(Y\vert X) = P(Y=1\vert X) = \Phi(\beta_0 + \beta_1 X), \label{eq:probitmodel}
\end{align}\] sodass der mit der Verknüpfungsfunktion (Link-Funktion) \[\textup{probit}(\cdot) = \Phi^{-1}(\cdot)\] transformierte Erwartungswert \(P(Y=1\vert X)\) dem linearen Prädiktor \(z=\beta_0 + \beta_1 X\) entspricht, \[\begin{align*}
\Phi^{-1}\big[P(Y=1\vert X)\big] = \beta_0 + \beta_1 X.
\end{align*}\]
\(z=\beta_0 + \beta_1 X\) in \(\eqref{eq:probitmodel}\) modelliert hier also ein Quantil der Standardnormalverteilung, \[\begin{align*}
\Phi(z) = P(Z \leq z) \ , \ Z \sim \mathcal{N}(0,1).
\end{align*}\] Der Koeffizient \(\beta_1\) in \(\eqref{eq:probitmodel}\) misst die Änderung in \(z\), die mit einer Änderung von \(X\) um eine Einheit verbunden ist. Obwohl der Effekt einer Änderung in \(X\) auf \(z\) linear ist, ist der Zusammenhang zwischen \(z\) und \(\textup{E}(Y\vert X) = P(Y=1\vert X)\)nicht linear, denn \(\Phi(\cdot)\) ist eine nicht-lineare Funktion von \(X\).
Aufgrund der Nicht-Linearität der Link-Funktion hat der Koeffizient von \(X\) keine einfache Interpretation hinsichtlich des Effekts auf \(P(Y=1\vert X)\). Die Änderung in der Wahrscheinlichkeit, dass \(Y=1\) ist, durch eine Änderung in \(X\) (partieller Effekt von \(X\)) kann berechnet werden als:
\[\begin{align*}
\frac{\partial\textup{E}(Y\vert X)}{\partial X} = \frac{\partial\textup{P}(Y=1\vert X)}{\partial X} = \frac{\partial\Phi(\beta_0 + \beta_1 X)}{\partial X} = \phi(\beta_0 + \beta_1 X) \beta_1,
\end{align*}\] wobei \(\phi(\cdot)\) die Dichtefunktion der Standardnormalverteilung ist. In empirischen Anwendungen wird der partielle Effekt häufig als Differenz in geschätzten Wahrscheinlichkeiten angegeben:
Berechne die geschätzte Wahrscheinlichkeit, dass \(Y=1\) für einen Bezugswert \(X\).
Berechne die geschätzte Wahrscheinlichkeit, dass \(Y=1\) für \(X + \Delta X\).
Berechne die Differenz zwischen der geschätzten Wahrscheinlichkeiten.
Wie im linearen Wahrscheinlichkeitsmodell kann das Modell \(\eqref{eq:probitmodel}\) auf eine Probit-Regression mit \(k+1\) Regressoren \(\boldsymbol{X} := (1, X_1, \dots, X_k)\) verallgemeinert werden, um das Risiko einer Verzerrung durch ausgelassene Variablen zu mindern. Die Schritte 1 bis 3 für die Berechnung des partiellen Effekts einer Änderung in \(X_j\) erfolgen dann unter der Annahme, dass die übrigen \(k-1\) Regressoren konstant gehalten werden, wobei der partielle Effekt von den jeweiligen Regressorwerten abhängt. Im nächsten Abschnitt erläutern wir die Schätzung von Probit-Modellen mit \(k+1\) Regressoren.
Hierbei ist \(\Phi(\mathbf{X}_i^\top \boldsymbol{\beta})\) die Wahrscheinlichkeit, dass \(Y_i = 1\) ist und \(1 - \Phi(\mathbf{X}_i^\top \boldsymbol{\beta})\) ist die Wahrscheinlichkeit, dass \(Y_i = 0\) ist.
Um den Maximum-Likelihood-Schätzer \(\widehat{\boldsymbol{\beta}}\) zu finden, muss die Log-Likelihood-Funktion \(\mathcal{L}(\boldsymbol{\beta})\) maximiert werden. In der Praxis erfolgt dies häufig durch numerische Optimierung, da die Log-Likelihood-Funktion der Probit-Regression im Allgemeinen keine geschlossene Form hat, sodass eine analytische Lösung nicht möglich ist.
Für eine Anwendung in R und einen Vergleich mit dem linearen Wahrscheinlichkeitsmodell erzeugen wir einen Beispieldatensatz simdata für einen datenerzeugenden Prozess (DGP)3 mit einem normalverteilten Regressor \(X\sim N(5,2^2)\) und \[\begin{align*}
P(Y=1\vert X) = \Phi(z), \quad z = -4 + 0.7 X.
\end{align*}\]
3 Siehe Kapitel 5 für Erläuterungen von Simulationsmethoden in R.
Das lineare Wahrscheinlichkeitsmodell schätzen wir wie gewohnt mit lm() und berechnen robuste Standardfehler mit lmtest::coeftest().
Beachte, dass lediglich die Vorzeichen der geschätzten Koeffizienten mit denen der wahren Werten übereinstimmmen. Da das lineare Modell fehlspezifiziert ist, sind die KQ-Schätzer der Koeffizienten hier inkonsistent.
Ein Probit-Modell kann mit stats::glm() geschätzt werden. Hierbei ist formula die Formel für den linearen Prädiktor \(z\) und family eine Link-Funktion für den Zusammenhang von \(z\) und \(P(Y\vert X)\). Mit family = binomial(link = "probit") wählen wir die Link-Funktion probit.
Die geschätzten Koeffizienten in der Probit-Regression liegen, wie erwartet, nahe bei Parametern des DGPs.
Um beide Schätzungen gemeinsam zu plotten, erzeugen wir mit predict() Vorhersagen für eine Menge von Werten X im Intervall \([0,11]\). Beachte, dass bei Vorhersagen für das Probit-Modell die gewünschte Transformation der vorherzusagenden Variable über type gewählt werden muss.4 Der Standardfall ist type = "link", d.h. wir erhalten Vorhersagen für den linearen Prädiktor: \(\widehat{z} = \widehat\beta_0 + \widehat{\beta}_1X\). Mit type = "response" werden diese Werte mit der Link-Funktion, hier \(\Phi(\cdot)\) zu vorhergesagten Wahrscheinlichkeiten \(\widehat{P}(Y=1\vert X = x)\) transformiert.
4 Dies gilt für jedes Modell mit einer anderen Link-Funktion als \(f(x) = x\) ist (lineare Regression).
Der nachfolgende Code-Chunk erstellt einen Punkte-Plot mit Jitter-Effekt (position_jitter), wobei die Beobachtungen zufällig leicht vertikal verschoben werden, um Überlappungen zu vermeiden. Zusätzlich zeichnen wir die geschätzten Regressionslinien für mod_lp (LPM) und mod_probit (Probit) sowie die tatsächliche Wahrscheinlichkeitsfunktion \(P(Y=1\vert X)\) ein.
Die Grafik zeigt, dass beide Modelle den positiven Zusammenhang für \(X\) und \(P(Y=1\vert X)\) erfassen. Das lineare Wahrscheinlichkeitsmodell approximiert die tatsächliche nicht-lineare Wahrscheinlichkeitsfunktion nur schlecht und liefert insbesondere in den Rändern der Verteilung von \(X\) (kleine und große Werte) unzulässige Vorhersagen außerhalb des Intervalls \([0,1]\). Die Probit-Spezifikation hingegen erfasst den nicht-linearen Verlauf der tatsächlichen Wahrscheinlichkeitsfunktion gut.
4.2.3 Logistische Regression
Bei logistischer Regression wird die logistische Funktion \(\Lambda(\cdot)\)\[\begin{align}
\Lambda(z) = \frac{1}{1 + \exp(-z)},
\end{align}\] als Link-Funktion genutzt, um die Wahrscheinlichkeitsfunktion von \(Y\) gegeben \(X\) zu modellieren. Ähnlich wie im Probit-Modell nehmen wir hier an, dass \[\begin{align}
E(Y\vert X) = P(Y=1\vert X) = \Lambda(\beta_0 + \beta_1 X). \label{eq:logitmodel}
\end{align}\] In diesem Modell ist die Link-Funktion \(\Lambda^{-1}(\cdot)\). Für \(z=\beta_0 + \beta_1 X\) ist \(\Lambda^{-1}(t)\) der sogenannte logit: Dass logarithmierte Verhältnis von \(p := P(Y=1\vert X)\) zu \(1 - p = P(Y=0\vert X)\), \[\begin{align*}
\textup{logit(p)} = \log\bigg(\frac{p}{1-p}\bigg) = \beta_0 + \beta_1 X.
\end{align*}\] Der Koeffizient \(\beta_1\) misst also die Veränderung des Logits pro Einheit Änderung im Regressor \(X\).
Ähnlich wie bei Probit-Regression ist der Einfluss von \(X\) auf den Logit linear, jedoch besteht auch hier eine nicht-lineare Beziehung zwischen dem linearen Prädiktor und der Wahrscheinlichkeit \(P(Y=1\vert X)\), denn \(\Lambda(\cdot)\) ist eine nicht-lineare Funktion mit dem Wertebereich \([0,1]\).
Die nachstehende interaktive Grafik illustriert, wie die Wahrscheinlichkeitsfunktion der latenten Variable \(z=\beta_0 + \beta_1 X\) von den Parametern \(\beta_0\) und \(\beta_1\) jeweils für Logit- und Probit-Regression beeinflusst wird.
Aufgrund der Nicht-Linearität von \(\Lambda(z)\) kann der Koeffizient \(\beta_1\) wie im Probit-Modell nicht direkt als Effekt auf die Wahrscheinlichkeit \(P(Y=1\vert X)\) interpretiert werden.
Um den partiellen Effekt einer Änderung in \(X\) am Punkt \(X\) auf \(P(Y=1\vert X)\) zu ermitteln, berechnen wir die Ableitung des bedingten Erwartungswerts:
\[\begin{align*}
\frac{\partial\textup{E}(Y\vert X)}{\partial X} = \frac{\partial\textup{P}(Y=1\vert X)}{\partial X} = \frac{\partial\Lambda(\beta_0 + \beta_1 X)}{\partial X} = \lambda(\beta_0 + \beta_1 X) \beta_1,
\end{align*}\] wobei \(\lambda(\cdot)\), ähnlich wie die Dichtefunktion der Normalverteilung im Probit-Modell, die Dichtefunktion der logistischen Verteilung darstellt. Diese ist gegeben durch \[\begin{align*}
\lambda(z) = \Lambda(z) \cdot (1 - \Lambda(z)).
\end{align*}\]
Die Interpretation des angepassten Modells ist aufgrund der Modellierung des Logits intuitiver als für Probit-Regression. Angenommen eine Bank möchte die Wahrscheinlichkeit modellieren, dass ein Kunde einen Kredit nicht zurückzahlt: \(P(Y=\textup{Zahlungsausfall}\vert X)\), wobei der Regressor \(X\) das Verhältnis von aktuellem Schuldenstand und Monatseinkommen ist. Die Schätzung einer logistischen Regression ergibt, dass ein Anstieg von \(X =x\) um 0.1 den Logit für die Wahrscheinlichkeit eines Kreditausfalls \(p(x)\) um 0.4 erhöht: \[\begin{align*}
\textup{logit} \big[P(Y=\textup{Zahlungsausfall}\vert X = x + 0.1)\big] = \textup{logit}[p(x)] + 0.4.
\end{align*}\]
Das Modell besagt dann, dass die Wahrscheinlichkeit der Zahlungsunfähigkeit etwa um den Faktor 1.5 ansteigt, denn \[\begin{align*}
\exp\big(\textup{logit}[p(x)] + 0.4\big) = \frac{p(x)}{1-p(x)} \cdot \exp(0.4)
\end{align*}\] mit \[\begin{align*}
\exp(0.4) \approx 1.50.
\end{align*}\]
4.2.4 Schätzung
Ähnlich wie bei Probit-Regressionen können Logit-Modelle mit Maximum-Likelihood geschätzt werden. Die Likelihood-Funktion für Logit-Regression lautet
mit \[\frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}\] der Wahrscheinlichkeit, dass \(Y_i = 1\). Für das Ereignis \(Y_i = 0\) ist die Wahrscheinlichkeit entsprechend \[1 - \frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}.\]
und erlaubt die Schätzung von \(\boldsymbol{\beta}\) durch Maximierung mit numerischen Methoden.
Für Anwendungen mit R nutzen wir stats::glm() und wählen die logistische Link-Funktion mit family = binomial(link = "logit"). Für die simulierten Daten simdata aus Kapitel 4.2.2:
Wir erweitern nun pred um die anhand von predict() für mod_logit geschätzten Wahrscheinlichkeiten von \(Y=1\vert X\) für die Regressorwerte in X.
Anhand der Vorhersagen in pred können wir die im Objekt p gespeicherte Grafik um die geschätzte Regressionsfunktion für das Logit-Modell erweitern.
Die Logit-Regression kann den wahren Zusammenhang gut abbilden und ist praktisch nicht von der Probit-Schätzung unterscheidbar.
4.2.5 Modellgüte
Ähnlich wie bei linearer Regression kann die Anpassung von generalisierten linearen Modellen an die Daten mit Anpassungsmaßen verglichen werden. In generalisierten linearen Modellen für binäre Outcome-Variablen kann hierzu die Deviance verwendet werden. Für eine Schätzung anhand der Beobachtungen \(i=1,\dots,n\) berechnet man die Deviance \(D\) als \[\begin{align*}
D = -2 \sum_{i=1}^{n} \left[ y_i \cdot \log\left(\widehat{P}_i\right) + (1 - y_i) \cdot \log\left(1 - \widehat{P}_i\right) \right].
\end{align*}\]\(D\) quantifiziert die Abweichung eines geschätzten Modells von einem perfekt passenden Modell.5 Eine geringere Deviance deutet darauf hin, dass das Modell die Daten besser erklärt. Eine R-Funktion zur Anwendung für Objekte des Typs glm ist deviance().
5 Im Allgemeinen vergleicht Deviance die Log-Likelihood des geschätzten Modells mit der Log-Likelihood des perfekten Modells.
Die Deviance für das Probit-Modell ist tatsächlich etwas geringer als für das Logit-Modell.
Ähnlich wie \(R^2\) für lineare Regression misst die Deviance \(D\) lediglich die Anpassung des Modells an die beobachteten Daten. Daher eignet sich \(D\) nur bedingt als Maß zur Einschätzung der Tauglichkeit eines Modells für Out-of-sample-Vorhersagen. Abhängig vom Ziel der Modellierung kann es hilfreich sein, eine robuste Einschätzung der Vorhersage-Güte unter Berücksichtigung der Modellkomplexität vorzunehmen. Hierfür werden oft Informationskriterien verwendet. Das Akaike-Informationskriterium (AIC) \[\begin{align*}
\text{AIC} = 2k - 2 \log(\widehat{\mathcal{L}})
\end{align*}\] ist ein Maß zur Bewertung der Güte eines Modells unter Berücksichtigung der Modellkomplexität. Das AIC wägt die Anpassung des Modells an die Daten (gemessen durch die maximierte Likelihood-Funktion \(\widehat{\mathcal{L}}\)6) und die Komplexität (gemessen als Funktion der Regressoranzahl \(k\)) gegeneinander ab. Hierbei werden Modelle mit weniger Parametern bevorzugt, wenn sie die Daten ähnlich gut erklären, sodass eine Überanpassung und damit eine schlechte Vorhersagefähigkeit vermieden werden. In R können wir das AIC für Modell-Objekte mit AIC() berechnen.
6\(\widehat{\mathcal{L}}\) meint die Likelihood-Funktion für die geschätzten Koeffizienten.
Auch das AIC zeigt an, dass die Schätzung des Probit-Modells in mod_probit besser geeignet für die Modiellierung von simdata ist als das Logit-Modell mod_logit.
4.2.6 Beispiel: Klassifikation von Palmer-Piniguinen
Mit den in diesem Kapitel betrachteten Methoden können wir Modelle zur Klassifikation von Pinguinen in palmerpenguins::penguins hinsichtlich ihres Geschlechts konstruieren. Hierfür nutzen wir den bereinigten Datensatz penguins_cleaned und transformieren zunächst die abhängige Variable sex in ein numerisches Format.
Wir schätzen nachfolgend ein lineares Wahrscheinlichkeitsmodell penguins_lp, ein Probit-Modell penguins_probit sowie ein Logit-Modell penguins_logit die jeweils sex anhand der simplen Spezikation des linearen Prädiktors \[z = \beta_0 + \beta_1 \textup{bill\_depth}\] erklären.
Mit modelsummary::modelsummary() können wir die wichtigsten Ergebnisse der Regressionen in einer publikationsfähigen Tabelle darstellen. Hierzu übergeben wir dem Argument models eine (benannte) Liste der geschätzten Modelle. Bei Angabe von vcov = "HC1" werden robuste Standardfehler ausgegeben.
Die geschätzten negativen und signifikanten Koeffizienten von \(\textup{bill\_depth}\) erfassen den Zusammenhang, dass Pinguine mit größeren Schnäbeln tendenziell männlich sind: Je größer \(\textup{bill\_depth}\), desto geringer die geschätzte Wahrscheinlichkeit, dass ein Pinguin weiblich ist. Ein Vergleich anhand des AIC zeigt, dass das Probit-Modell penguins_probit am besten geeignet ist.
Analog zur Vorgehensweise in Kapitel 4.2.3 können wir die geschätzten Modelle vergleichen, in dem wir geschätzte Wahrscheinlichkeiten \(P(Y=\textup{weiblich}\vert\textup{bill\_depth})\) für eine Menge repräsentativer Werte des Regressors bill_depth berechnen und gemeinsam mit den Beobachtungen plotten.
Die geschätzten Regressionsfunktionen unterscheiden sich für den Bereich beobachteter Werte von \(bill\_depth\) nur wenig. Die Grafik zeigt, dass die Schätzungen mit Probit- und Logit-Ansatz nur wenig Nicht-Linearität aufweisen, sodass eine (leichter interpretierbare) Modellierung mit dem linearen Wahrscheinlichkeitsmodell penguins_lp hier zulässig scheint.
Die Klassifizierung der Pinguine im Datensatz hinsichtlich ihres Geschlechts anhand vorhergesagter Wahrscheinlichkeiten \(\widehat{P}_i\) kann durch Abgleich mit einem Grenzwert erfolgen. Ein Beispiel ist \[\begin{align}
\widehat{Y}_i = \begin{cases}
\textup{weiblich}, & \widehat{P}_i \geq .5,\\
\textup{männlich}, & \textup{sonst.}
\end{cases}\label{eq:pengclass}
\end{align}\]
Wir können diesen Ansatz in R implementieren, indem wir zunächst die für den Datensatz angepassten Wahrscheinlichkeiten fitted aus den Modell-Objekten auslesen.
Für die Klassifikation der Pinguine anhand der jeweiligen geschätzten Wahrscheinlichkeiten wenden wir die Regel \(\eqref{eq:pengclass}\) mit across() spaltenweise für die Modelle an und erzeugen neue Spalten mit vorhersagen des Geschlechts. Mit .names = "{.col}_pred" erhalten die neuen Spalten das Suffix _pred.
Mit summarise() berechnen wir nun den Anteil korrekter Vorhersagen.
In diesem Beispiel unterscheiden sich die Vorhersagen der drei Modelle für den Grenzwert 0.5 nicht.
4.3 Modellierung von Zählvariablen
Eine weitere Klasse von \(Y\) für die sich der Regressionsansatz von einfacher linearer Regression unterscheidet, sind Zählvariablen: Variablen, die diskrete, nicht-negative Werte annehmen. Wenn die abhängige Variable \(Y\) Ereignisse in einem bestimmten Zeitraum misst (wie z.B. die Anzahl der Verkehrsunfälle in Essen innerhalb eines Monats) und weitere Bedingungen erfüllt sind, folgt \(Y\) einer Poisson-Verteilung und kann mit Poisson-Regression modelliert werden. Wir erläutern nachfolgend die wesentlichen Komponenten von Poisson-Regression und diskutieren ein R-Beispiel mit fiktiven Daten.
4.3.1 Poisson-Verteilung
Die Zufallsvariable \(Y\) folgt einer Poisson-Verteilung mit Parameter \(\lambda\), wenn
die Ereignisse unabhängig voneinander auftreten: Das Auftreten eines Ereignisses beeinflusst nicht die Wahrscheinlichkeit, dass ein weiteres Ereignis auftritt.
die Ereignisse einzeln auftreten: Die Wahrscheinlichkeit, dass in einem sehr kleinen Intervall (im Sinne von \([a,b]\) für \(a \to b\)) mehr als ein Ereignis auftritt, ist vernachlässigbar.
die Rate der Ereignisse konstant ist: Die durchschnittliche Anzahl der Ereignisse pro Zeiteinheit oder pro Raumeinheit \(\lambda\) ist konstant.
Die Wahrscheinlichkeitsfunktion von \(Y\) ist
\[\begin{align}
P(Y = y) = \frac{\lambda^y e^{-\lambda}}{y!} \quad \text{für} \quad y = 0, 1, 2, \ldots,\label{eq:poissonpmf}
\end{align}\] wobei \(\lambda\) sowohl der Erwartungswert als auch die Varianz der Verteilung ist, \[\textup{E}(Y) = \textup{Var}(Y) = \lambda.\]
Die nächste Grafik zeigt die diskrete (Wahrscheinlichkeitsmasse-)Funktion \(\eqref{eq:poissonpmf}\) für eine Poisson-verteilte Zufallsvariable \(Y\) mit \(\lambda = 5\) für den Wertebereich von 0 bis 15.
4.3.2 Poisson-Regression
Wie für die bisher betrachteten Regressionsansätzen modelliert Poisson-Regression den Erwartungswert (und damit gleichzeitig die Varianz) \(\lambda\) der abhängigen Variable \(Y\) als eine Funktion der Regressoren \(\mathbf{X} = (X_1, X_2, \ldots, X_k)\). Unter Annahme einer Poisson-Verteilung kann Poisson-Regression als generalisiertes lineares Modell verstanden werden, wobei eine logarithmische Verknüpfungsfunktion für den linearen Prädiktor \(\mathbf{X}_i^\top \boldsymbol{\beta}\) verwendet wird:
Den ML-Schätzer \(\widehat{\boldsymbol{\beta}}\) erhalten wir durch Maximierung der Log-Likelihoodfunktion \(\mathcal{L}(\boldsymbol{\beta})\). Die R-Implementierung von Poisson-Regression in stats::glm() wird mittels family = poisson(link = "log") aufgerufen.
4.3.4 Interpretation der Koeffizienten
Die Koeffizienten \(\boldsymbol{\beta}\) in der Poisson-Regression haben eine logarithmisch-lineare Beziehung zur Zählvariable. Für einen bestimmten Koeffizienten \(\beta_j\) ist die Interpretation wie folgt:
Eine Änderung der unabhängigen Variable \(X_j\) um eine Einheit führt zu einer Änderung des Logarithmus des Erwartungswertes von \(Y\) um \(\beta_j\). Der Erwartungswert \(\lambda\) ändert sich also multiplikativ um den Faktor \(\exp(\beta_j)\).
Als Beispiel betrachten wir den linearen Prädiktor
Angenommen \(X\) ist die Anzahl durchgeführter Werbekampagnen und die abhängige Zählvariable \(Y\) misst die Anzahl der Verkäufe des beworbenen Produkts pro Monat.
Wenn \(\beta_1 = 0.5\), bedeutet dies, dass jede zusätzliche Werbekampagne (\(\Delta X = 1\)) die erwartete Anzahl der Verkäufe pro Monat um einen Faktor von \(\exp(0.5) \approx 1.65\) erhöht:
Das heißt, die Rate der Verkäufe steigt um 65% für jede zusätzliche Werbekampagne.
Zur Veranschaulichung der Schätzung einer Poisson-Regression für dieses Beispiel erzeugen wir Poisson-verteilte Daten für \[\begin{align*}
\log(\lambda) = 2 + 0.4 \cdot X
\end{align*}\] und ziehen \(X\) gleichverteilt aus \(\{1,2,\dots,8\}\).
Die Züge aus der abhängigen Variablen visualieren wir mit einem Häufigkeits-Histogramm.
Wir schätzen das Modell und extrahieren die robuste Zusammenfassung mit broom::tidy().
Die ML-Schätzung des Poisson-Modells liefert für die verwendete Stichprobe von 500 Beobachtungen Koeffizienten-Schätzungen nahe der wahren Parameter. Zur Veranschaulichung des geschätzten Modells überlagern wir die simulierten Datenpunkte aus dat mit der anhand von mod_poisson geschätzten Anzahl der Verkäufe je Anzahl an Werbekampagnen.
4.4 Zusammenfassung
In diesem Kapitel haben wir erweiterte Konzepte der Regressionsanalyse diskutiert. Anhand des Frisch-Waugh-Lovell-Theorems wurde die Wirkungsweise der Kontrolle von Kovariablen mit multipler Regression auf den interessierenden Koeffizienten (Effekt) einer Variable veranschaulicht. Weiterhin haben wir gängige Typen generalisierter linearer Modelle für binäre und Poisson-verteilte Outcome-Variablen eingeführt. Anhand relevanter Beispiele wurde gezeigt, wie diese generalisierten Regressionsmodelle in R spezifiziert, geschätzt und der resultierende Output mit Paketen wie broom, ggplot2 und modelsummary interpretiert werden kann.
Quellcode
---format: live-htmlengine: knitrwebr: packages: [ 'broom', 'cowplot', 'lmtest', 'dplyr', 'ggplot2', 'gt', 'modelsummary', 'palmerpenguins', 'purrr', 'tidyr', 'readr', 'sandwich' ]---{{< include ./_extensions/r-wasm/live/_knitr.qmd >}}# Regression {#sec-regression}Viele der in diesem Companion behandelten Methoden basiert auf dem Konzept der Schätzung kausaler Effekte mit *Regression*. Als Regressionsansatz bezeichnet man eine Methode, welche die Beziehungen zwischen Variablen durch einen funktionalen Zusammenhang beschreibt und die Parameter der gewählten funktionalen Form anhand von beobachteten Daten schätzt. *Lineare Regression* nimmt eine lineare funktionale Form der Beziehung zwischen einer abhängigen Variable (Outcome-Variable) und erklärenden Variablen (Regressoren) an. *Nicht-lineare Regressionsmethoden* modellieren die Beziehung etwa durch Polynome höherer Ordnung, exponentielle Funktionen oder andere komplexere Formen. Die Wahl der funktionalen Form hängt von der Natur des datenerzeugenden Prozesses (DGP) und somit stets von der spezifischen Beziehung ab, die untersucht wird. Regressionsansätze gehören zu den am häufigsten verwendeten Methoden für Kausalanalysen, da mit Regression in vielen Forschungsdesigns kausale Effekte identifiziert werden können, indem Backdoors geschlossen werden: Regression kann die durch die Behandlungsvariable verursachte Variation in der Outcome-Variable isolieren, indem für gemeinsame Einflussfaktoren von Behandlungs- *und* Outcome-Variable kontrolliert wird. In diesem Kapitel erläutern wir Erweiterungen der Spezifikation und Schätzung von Regressionsansätzen, die für spätere Kapitel relevant sind. Neben einer Motivation der Schätzung kausaler Effekte mit multipler Regression betrachten wir Modelle für verschiedene Kategorien von Outcome-Variablen und diskutieren deren Implementierung mit R.```{r}library(tidyverse)library(cowplot)```## Regression schließt Backdoors: Frisch-Waugh-Lovell-TheoremDas Frisch-Waugh-Lovell-Theorem (FWL) besagt, dass die geschätzten Koeffizienten für eine interessierende Teilmenge der Regressoren in einer multiplen linearen Regression numerisch identisch mit den Koeffizientenschätzungen aus folgenden Schritten sind: 1. Rechne die Effekte der übrigen Regressoren auf (a) die Outcome-Variable und (b) die interessierende Teilmenge der Regressoren mit Regression heraus. 2. Regressiere anschließend die Residuen von Schritt (a) auf die Residuen aus Schritt (b). In einem multiplen Modell mit zwei Regressoren $X_1,\ X_2$,\begin{align} Y = \beta_0 + \beta_1 X + \beta_2 X_2 + \epsilon \label{eq:fwlfullreg}\end{align}kann der Effekt von $X_1$ auf $Y$ also mit der Regression\begin{align*} \widehat{u}_{Y,X_2} = \beta_1 \widehat{u}_{X_1,X_2} + e \label{eq:fwl2reg}\end{align*}geschätzt werden, wobei $\widehat{u}_{Y,X_2}$ und $\widehat{u}_{X_1,X_2}$ die Residuen der Regression von $Y$ auf $X_2$ und von $X_1$ auf $X_2$ sind.FWL ermöglicht daher eine Vereinfachung der Schätzung komplexer Modelle durch die Zerlegung der Schätzung in Teilschritte. Für das Verständnis der Schätzung kausaler Effekte mit linearer Regression ist FWL hilfreich, denn es zeigt, wie sowohl die Variation in der Outcome-Variable ($\widehat{u}_{Y,X_2}$) als auch die Variation in der Behandlungsvariable ($\widehat{u}_{X_1,X_2}$), die jeweils *nicht* durch Kovariablen ($X_2$) verursacht wird, mit multipler Regression isoliert werden kann, sodass Backdoors geschlossen werden.Wir illustrieren dieses Konzept anhand einer multiplen Regression für das Gewicht (`body_mass`) von Pinguinen aus dem Datensatz `palmerpenguins::penguins`,\begin{align} \textup{body\_mass} = \beta_0 + \beta_1\cdot\textup{bill\_length} + \beta_2\cdot \textup{flipper\_length} + \epsilon,\label{eq:billdepthmodel}\end{align}unter der Annahme, dass $\beta_1$ der interessierende Effekt ist: Die erwartete Änderung des Gewichts eines Pinguins (in Gramm) für eine Änderung der Schnabel-Länge um 1mm.Vor der Schätzung von Modell \eqref{eq:billdepthmodel} lesen wir den Datensatz ein und erstellen eine bereinigte Variante `penguins_cleaned`, analog zur Vorgehensweise in @sec-pp.```{webr}library(palmerpenguins)data(penguins)# Datensatz bereinigenpenguins_cleaned <- penguins %>% rename( bill_depth = bill_depth_mm, bill_length = bill_length_mm, flipper_length = flipper_length_mm, body_mass = body_mass_g ) %>% drop_na() %>% filter( body_mass < quantile(body_mass, .95) )# Überblickslice_head(penguins_cleaned, n = 10)```Wir schätzen nun Modell \eqref{eq:billdepthmodel} mit `lm()` und erhalten eine Zusammenfassung der geschätzten Koeffizienten mit `broom::tidy()`.```{webr}library(broom)# "Großes" Modell schätzenlm( formula = body_mass ~ bill_length + flipper_length, data = penguins_cleaned) %>% tidy()```Das Ergebnis der Schätzung ist $\widehat{\beta}_1\approx3.80$. Der nächste Code-Block berechnet die Residuen aus den Regressionen\begin{align*} \textup{body\_mass} =&\, \alpha_0 + \alpha_1 \textup{flipper\_length} + u_{\textup{body\_mass},\,\textup{flipper\_length}},\\ \textup{bill\_length} =&\, \delta_0 + \delta_1 \textup{flipper\_length} + u_{\textup{bill\_length},\,\textup{flipper\_length}},\end{align*}und speichert diese in `body_mass_res` und `bill_length_res`.```{webr}# FWL-Schritt 1 (a)peng_mod_fwl1a <- lm( formula = body_mass ~ flipper_length, data = penguins_cleaned)body_mass_res <- residuals(peng_mod_fwl1a)head(body_mass_res, n = 10)``````{webr}# FWL-Schritt 1 (b)peng_mod_fwl1b <- lm( formula = bill_length ~ flipper_length, data = penguins_cleaned)bill_length_res <- residuals(peng_mod_fwl1b)head(bill_length_res, n = 10)```Für den zweiten Schritt regressieren wir `body_mass_res` auf `bill_length_res`.```{webr}# FWL-Schritt 2lm(formula = body_mass_res ~ bill_length_res - 1) %>% tidy()```Der geschätzte Koeffizient aus der Regression der Residuen stimmt mit dem geschätzten Koeffizienten von `bill_length` aus der großen Regression \eqref{eq:billdepthmodel} überein.Wir können den Effekt der Kontrolle für `flipper_length` visualisieren. Wir plotten hierzu:- Die originalen Datenpunkte für `bill_length` und `body_mass`^[Für eine bessere Lesbarkeit der Grafik zentrieren wir beide Variablen um den jeweiligen Stichprobenmittelwert.] gemeinsam mit der geschätzten Regressionslinie für das Modell $$ \textup{body\_mass} = \beta_0 + \beta_1\textup{bill\_length} + u $$ (keine Kontrolle für `flipper_length`!)^[Der R-Befehl für diese Regression ist `lm(I(body_mass - mean(body_mass)) ~ I(bill_length - mean(bill_length)) - 1, data = penguins_cleaned)`.].- Die um `flipper_length` bereinigten Datenpunkte und die zugehörige geschätzte Regressionslinie.```{webr}# Residuen in tibble sammelnp <- tibble( body_mass_res, bill_length_res) %>% ggplot() + # Bereinigte Datenpunkte plotten geom_point( mapping = aes( x = bill_length_res, y = body_mass_res, color = "Um flipper_length bereinigte Daten" ) ) + # Regressionslinie für bereinigte Datenpunkte geom_smooth( mapping = aes(x = bill_length_res, y = body_mass_res), method = "lm", formula = "y ~ x - 1", se = F, col = "purple" ) + scale_color_manual( name = "", values = c( "Um flipper_length bereinigte Daten" = "purple", "Ursprüngliche Datenpunkte" = "black" ) ) + labs( title = "Penguins: Anwendung von FWL", x = "bill_length", y = "body_mass" ) + theme_cowplot() + theme(legend.position = "top")p```Wir erweitern `p` um die ursprünglichen Datenpunkte und die zugehörige Regressionslinie.```{webr}# Ursprüngliche Datenpunkte hinzufügenp + geom_point( data = penguins_cleaned, mapping = aes( x = bill_length - mean(bill_length), y = body_mass - mean(body_mass), color = "Ursprüngliche Datenpunkte" ), alpha = .5, ) + # Regressionslinie für ursprüngliche Datenpunkte geom_smooth( data = penguins_cleaned, mapping = aes( x = bill_length - mean(bill_length), y = body_mass - mean(body_mass) ), method = "lm", formula = "y ~ x - 1", se = F, col = alpha("black", .5) )```Der grafische Vergleich beider Vorgehensweisen zeigt den Effekt der Kontrolle für `flipper_length`: Die geschätzte (schwarze) Regressionslinie für die bereinigten Daten hat eine deutlich geringere Steigung als die anhand der ursprünglichen Daten geschätzte (lilane) Linie. Der Effekt von `bill_length` auf `body_mass` wird mit der einfachen Regression `lm(body_mass ~ bill_length)` vermutlich *überschätzt*, weil es andere Faktoren (wie `flipper_length` gibt, die mit `bill_length` und `body_mass` korrelieren. Kontrollieren für `flipper_length` in der multiplen Regression `lm(body_mass ~ bill_length + flipper_length)` schließt die Backdoor durch `flipper_length`. Die Konsequenz ist eine deutlich geringere Steigung der lilanen Regressionslinie.## Binäre Outcome-Variable {#sec-bov}Eine binäre Variable, auch als dichotome Variable oder Indikator-Variable bezeichnet, ist eine Variable, die nur zwei Ausprägungen annehmen kann. Diese beiden Ausprägungen werden typischerweise durch die Werte 0 und 1 repräsentiert und dienen dazu, zwei verschiedene Zustände oder Kategorien zu unterscheiden. Formal kann eine binäre Variable $B$ wie folgt definiert werden:\begin{align} B = \begin{cases} 1, & \text{Eigenschaft trifft zu,} \\ 0, & \text{Eigenschaft trifft nicht zu.} \end{cases}\end{align}Ein in späteren Kapiteln dieses Companions verwendeter *binärer Regressor* ist der Indikator für die Zuordnung von Beobachtungen zur Behandlungs- oder Kontrollgruppe (1 = Behandlungsgruppe, 0 = Kontrollgruppe). Für viele ökonomische Forschungsfragen ist es hilfreich, eine *binäre Outcome-Variable* mit Regression zu modellieren. Hierzu gibt es verschiedene Ansätze, die wir nachfolgend zusammenfassen und ihre Anwendung mit R zeigen.### Das lineare Wahrscheinlichkeitsmodell {#sec-lpm}Das lineare Regressionsmodell$$Y = \beta_0 + \beta_1 X_{1} + \beta_2 X_{2} + \dots + \beta_k X_{k} + u$$mit einer binären abhängigen Variablen $Y_i\in\{0,1\}$ wird als *lineares Wahrscheinlichkeitsmodell* bezeichnet. Wie üblich modellieren wir den Erwartungswert der abhängigen Variable gegeben der Regressoren $X_1,\dots,X_k$ als lineare Funktion,$$E(Y\vert X_1,X_2,\dots,X_k) = P(Y=1\vert X_1, X_2,\dots, X_3).$$ Da $Y$ eine binäre Variable ist, gilt hier$$ P(Y = 1 \vert X_1, X_2, \dots, X_k) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k.$$Das lineare Wahrscheinlichkeitsmodell beschreibt also die *Wahrscheinlichkeit*, dass $Y=1$ als lineare Funktion der Regressoren: $\beta_j$ misst die Änderung in der Wahrscheinlichkeit für das Ereignis $Y_i=1$, unter der Bedingung, dass die anderen $k-1$ Regressoren konstant gehalten werden. Genau wie bei multipler Regression mit einer kontinuierlichen abhängigen Variablen können die $\beta_j$ mit der KQ-Methode geschätzt werden.Aufgrund der Beschränktheit der $Y_i$ auf $\{0,1\}$ ist $u_i$ heteroskedastisch. Folglich sollten Inferenzstatistiken mit robusten Standardfehlern berechnet werden. Weiterhin ist zu beachten, dass $R^2$ in den meisten Anwendungen von linearen Wahrscheinlichkeitsmodellen keine hilfreiche Interpretation hat, da das geschätzte Modell die Daten nicht perfekt erklären kann, wenn die abhängige Variable binär, aber die Regressoren kontinuierlich verteilt sind.Das lineare Wahrscheinlichkeitsmodell hat einen wesentlichen Nachteil: Das Modell nimmt an, dass die bedingte Wahrscheinlichkeitsfunktion linear ist und $P(Y=1\vert X_1,\dots,X_k)$ über das für Wahrscheinlichkeiten definierte Intervall $[0,1]$ hinausgehen kann. Ein angepasstes Modell hat dann für Regressorwerte, die zu Vorhersagen von $Y$ jenseits von $[0,1]$ führen keine sinnvolle Interpretation.Dieser Umstand verlangt nach Regressionsansätzen, die $P(Y=1)$ durch eine auf $[0,1]$ beschränkte (nicht-lineare) Funktion der Regressoren modellieren. Häufig verwendete Methoden sind Probit- und Logit-Regression.Ein lineares Wahrscheinlichkeitsmodell kann mit `lm()` geschätzt werden, wobei die abhängige Variable den Typ `numeric` oder `integer` haben muss.### Probit-Regression {#sec-probitreg}Bei der Probit-Regression wird die Standardnormalverteilungsfunktion $\Phi(\cdot)$ verwendet, um die Regressionsfunktion einer binären abhängigen Variable zu modellieren. Wir nehmen an, dass\begin{align} E(Y\vert X) = P(Y=1\vert X) = \Phi(\beta_0 + \beta_1 X), \label{eq:probitmodel}\end{align}sodass der mit der Verknüpfungsfunktion (Link-Funktion) $$\textup{probit}(\cdot) = \Phi^{-1}(\cdot)$$ transformierte Erwartungswert $P(Y=1\vert X)$ dem linearen Prädiktor $z=\beta_0 + \beta_1 X$ entspricht,\begin{align*} \Phi^{-1}\big[P(Y=1\vert X)\big] = \beta_0 + \beta_1 X.\end{align*}$z=\beta_0 + \beta_1 X$ in \eqref{eq:probitmodel} modelliert hier also ein *Quantil* der Standardnormalverteilung,\begin{align*}\Phi(z) = P(Z \leq z) \ , \ Z \sim \mathcal{N}(0,1).\end{align*}Der Koeffizient $\beta_1$ in \eqref{eq:probitmodel} misst die Änderung in $z$, die mit einer Änderung von $X$ um eine Einheit verbunden ist. Obwohl der Effekt einer Änderung in $X$ auf $z$ linear ist, ist der Zusammenhang zwischen $z$ und $\textup{E}(Y\vert X) = P(Y=1\vert X)$ *nicht linear*, denn $\Phi(\cdot)$ ist eine nicht-lineare Funktion von $X$.```{webr}# N(0,1)-Quantilsfunktion / Probit-Linkggplot() + geom_function(fun = qnorm) + scale_x_continuous( name = "z", limits = c(0, 1) ) + scale_y_continuous(name = "Quantil") + labs(title = "Probit-Funktion") + theme_cowplot()```Aufgrund der Nicht-Linearität der Link-Funktion hat der Koeffizient von $X$ keine einfache Interpretation hinsichtlich des Effekts auf $P(Y=1\vert X)$. Die Änderung in der Wahrscheinlichkeit, dass $Y=1$ ist, durch eine Änderung in $X$ (partieller Effekt von $X$) kann berechnet werden als:\begin{align*} \frac{\partial\textup{E}(Y\vert X)}{\partial X} = \frac{\partial\textup{P}(Y=1\vert X)}{\partial X} = \frac{\partial\Phi(\beta_0 + \beta_1 X)}{\partial X} = \phi(\beta_0 + \beta_1 X) \beta_1,\end{align*}wobei $\phi(\cdot)$ die Dichtefunktion der Standardnormalverteilung ist. In empirischen Anwendungen wird der partielle Effekt häufig als Differenz in geschätzten Wahrscheinlichkeiten angegeben:1. Berechne die geschätzte Wahrscheinlichkeit, dass $Y=1$ für einen Bezugswert $X$.2. Berechne die geschätzte Wahrscheinlichkeit, dass $Y=1$ für $X + \Delta X$.3. Berechne die Differenz zwischen der geschätzten Wahrscheinlichkeiten.Wie im linearen Wahrscheinlichkeitsmodell kann das Modell \eqref{eq:probitmodel} auf eine Probit-Regression mit $k+1$ Regressoren $\boldsymbol{X} := (1, X_1, \dots, X_k)$ verallgemeinert werden, um das Risiko einer Verzerrung durch ausgelassene Variablen zu mindern. Die Schritte 1 bis 3 für die Berechnung des partiellen Effekts einer Änderung in $X_j$ erfolgen dann unter der Annahme, dass die übrigen $k-1$ Regressoren konstant gehalten werden, wobei der partielle Effekt von den jeweiligen Regressorwerten abhängt. Im nächsten Abschnitt erläutern wir die Schätzung von Probit-Modellen mit $k+1$ Regressoren.#### SchätzungDie Likelihood-Funktion für Probit-Regression ist\begin{align*} L(\boldsymbol{\beta}) = \prod_{i=1}^n \Phi(\mathbf{X}_i^\top \boldsymbol{\beta})^{y_i} \left[1 - \Phi(\mathbf{X}_i^\top \boldsymbol{\beta})\right]^{1-y_i}\end{align*}Hierbei ist $\Phi(\mathbf{X}_i^\top \boldsymbol{\beta})$ die Wahrscheinlichkeit, dass $Y_i = 1$ ist und $1 - \Phi(\mathbf{X}_i^\top \boldsymbol{\beta})$ ist die Wahrscheinlichkeit, dass $Y_i = 0$ ist.Die Log-Likelihood-Funktion ergibt sich als\begin{align*} \mathcal{L}(\boldsymbol{\beta}) = \sum_{i=1}^n \left[ y_i \log \Phi(\mathbf{X}_i^\top \boldsymbol{\beta}) + (1 - y_i) \log \left(1 - \Phi(\mathbf{X}_i^\top \boldsymbol{\beta})\right) \right]\end{align*}Um den Maximum-Likelihood-Schätzer $\widehat{\boldsymbol{\beta}}$ zu finden, muss die Log-Likelihood-Funktion $\mathcal{L}(\boldsymbol{\beta})$ maximiert werden. In der Praxis erfolgt dies häufig durch numerische Optimierung, da die Log-Likelihood-Funktion der Probit-Regression im Allgemeinen keine geschlossene Form hat, sodass eine analytische Lösung nicht möglich ist.Für eine Anwendung in R und einen Vergleich mit dem linearen Wahrscheinlichkeitsmodell erzeugen wir einen Beispieldatensatz `simdata` für einen datenerzeugenden Prozess (DGP)^[Siehe @sec-sim für Erläuterungen von Simulationsmethoden in R.] mit einem normalverteilten Regressor $X\sim N(5,2^2)$ und\begin{align*} P(Y=1\vert X) = \Phi(z), \quad z = -4 + 0.7 X.\end{align*}```{webr}# Daten simulierenset.seed(1234)# Stichprobengrößen <- 500 simdata <- tibble( X = rnorm(n = n, mean = 5, sd = 2), # Regressor P = pnorm(-4 + 0.7 * X),)# Binäre Outcome-Variable hinzufügensimdata <- simdata %>% mutate( Y = as.integer(runif(n) < P) )# Überblickslice_head(simdata, n = 10)```Das lineare Wahrscheinlichkeitsmodell schätzen wir wie gewohnt mit `lm()` und berechnen robuste Standardfehler mit `lmtest::coeftest()`.```{webr}# lineares Wahrscheinlichkeitsmodell schätzenmod_lp <- lm( formula = Y ~ X, data = simdata)# Zusammenfassungmod_lp %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()```Beachte, dass lediglich die Vorzeichen der geschätzten Koeffizienten mit denen der wahren Werten übereinstimmmen. Da das lineare Modell fehlspezifiziert ist, sind die KQ-Schätzer der Koeffizienten hier inkonsistent.Ein Probit-Modell kann mit `stats::glm()` geschätzt werden. Hierbei ist `formula` die Formel für den linearen Prädiktor $z$ und `family` eine Link-Funktion für den Zusammenhang von $z$ und $P(Y\vert X)$. Mit `family = binomial(link = "probit")` wählen wir die Link-Funktion `probit`.```{webr}# Probit-Modell schätzenmod_probit <- glm( formula = Y ~ X, data = simdata, family = binomial(link = "probit"))# Zusammenfassungmod_probit %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()```Die geschätzten Koeffizienten in der Probit-Regression liegen, wie erwartet, nahe bei Parametern des DGPs.Um beide Schätzungen gemeinsam zu plotten, erzeugen wir mit `predict()` Vorhersagen für eine Menge von Werten `X` im Intervall $[0,11]$. Beachte, dass bei Vorhersagen für das Probit-Modell die gewünschte Transformation der vorherzusagenden Variable über `type` gewählt werden muss.^[Dies gilt für jedes Modell mit einer anderen Link-Funktion als $f(x) = x$ ist (lineare Regression).] Der Standardfall ist `type = "link"`, d.h. wir erhalten Vorhersagen für den linearen Prädiktor: $\widehat{z} = \widehat\beta_0 + \widehat{\beta}_1X$. Mit `type = "response"` werden diese Werte mit der Link-Funktion, hier $\Phi(\cdot)$ zu vorhergesagten Wahrscheinlichkeiten $\widehat{P}(Y=1\vert X = x)$ transformiert.```{webr}# geschätzte Wahrscheinlichkeitsfunktion# für lineares ModellX <- seq(0, 11, 0.01)pred <- tibble( X = X, LP = predict( object = mod_lp, newdata = tibble(X) ), probit = predict( object = mod_probit, newdata = tibble(X), type = "response" ))# Überblick verschaffenslice_head(pred, n = 10)```Der nachfolgende Code-Chunk erstellt einen Punkte-Plot mit Jitter-Effekt (`position_jitter`), wobei die Beobachtungen zufällig leicht vertikal verschoben werden, um Überlappungen zu vermeiden. Zusätzlich zeichnen wir die geschätzten Regressionslinien für `mod_lp` (LPM) und `mod_probit` (Probit) sowie die tatsächliche Wahrscheinlichkeitsfunktion $P(Y=1\vert X)$ ein.```{webr}# Daten und geschätzte Modelle plottenp <- simdata %>% ggplot(mapping = aes(x = X, y = Y)) + geom_point( position = position_jitter( height = .025, seed = 1234 ), alpha = .25, color = "gray" ) + geom_line( data = pred, mapping = aes(y = LP, color = "LPM"), lwd = .75 ) + geom_line( data = pred, mapping = aes(y = probit, color = "Probit"), lwd = .75 ) + geom_function( fun = \(x) pnorm(-4 + .7 * x), mapping = aes(color = "Wahrheit") ) + scale_color_manual( name = "", values = c( "Wahrheit" = "black", "LPM" = "steelblue", "Probit" = "orange" ) ) + labs(title = "Gesch. Modelle für binäre abhängige Variable") + theme_cowplot() + theme(legend.position = "top")p```Die Grafik zeigt, dass beide Modelle den positiven Zusammenhang für $X$ und $P(Y=1\vert X)$ erfassen. Das lineare Wahrscheinlichkeitsmodell approximiert die tatsächliche nicht-lineare Wahrscheinlichkeitsfunktion nur schlecht und liefert insbesondere in den Rändern der Verteilung von $X$ (kleine und große Werte) unzulässige Vorhersagen außerhalb des Intervalls $[0,1]$. Die Probit-Spezifikation hingegen erfasst den nicht-linearen Verlauf der tatsächlichen Wahrscheinlichkeitsfunktion gut.### Logistische Regression {#sec-logreg}Bei logistischer Regression wird die logistische Funktion $\Lambda(\cdot)$ \begin{align}\Lambda(z) = \frac{1}{1 + \exp(-z)},\end{align}als Link-Funktion genutzt, um die Wahrscheinlichkeitsfunktion von $Y$ gegeben $X$ zu modellieren. Ähnlich wie im Probit-Modell nehmen wir hier an, dass\begin{align}E(Y\vert X) = P(Y=1\vert X) = \Lambda(\beta_0 + \beta_1 X). \label{eq:logitmodel}\end{align}In diesem Modell ist die Link-Funktion $\Lambda^{-1}(\cdot)$. Für $z=\beta_0 + \beta_1 X$ ist $\Lambda^{-1}(t)$ der sogenannte *logit*: Dass logarithmierte Verhältnis von $p := P(Y=1\vert X)$ zu $1 - p = P(Y=0\vert X)$,\begin{align*} \textup{logit(p)} = \log\bigg(\frac{p}{1-p}\bigg) = \beta_0 + \beta_1 X.\end{align*}Der Koeffizient $\beta_1$ misst also die Veränderung des Logits pro Einheit Änderung im Regressor $X$. Ähnlich wie bei Probit-Regression ist der Einfluss von $X$ auf den Logit linear, jedoch besteht auch hier eine nicht-lineare Beziehung zwischen dem linearen Prädiktor und der Wahrscheinlichkeit $P(Y=1\vert X)$, denn $\Lambda(\cdot)$ ist eine nicht-lineare Funktion mit dem Wertebereich $[0,1]$.```{webr}# Plot: Logistische Funktionggplot() + stat_function( fun = \(x) 1 / (1 + exp(-x)) ) + scale_x_continuous( name = "z", limits = c(-4, 4) ) + scale_y_continuous(name = "Lambda(z)") + labs(title = "Logistische Funktion") + theme_cowplot()``````{webr}# Plot: Logitggplot() + stat_function( fun = \(x) log(x / (1 - x)) ) + scale_x_continuous( name = "p", limits = c(0, 1) ) + scale_y_continuous(name = "Logit(p)") + labs(title = "Logit-Funktion") + theme_cowplot()```Die nachstehende interaktive Grafik illustriert, wie die Wahrscheinlichkeitsfunktion der latenten Variable $z=\beta_0 + \beta_1 X$ von den Parametern $\beta_0$ und $\beta_1$ jeweils für Logit- und Probit-Regression beeinflusst wird.<iframe class="obs-soft-box-shadow" width="100%" height="680" frameborder="0" src="https://observablehq.com/embed/@mca91/latent-variable-cdf?cells=plot%2Cviewof+beta0%2Cviewof+beta%2CMathJax"></iframe>Aufgrund der Nicht-Linearität von $\Lambda(z)$ kann der Koeffizient $\beta_1$ wie im Probit-Modell nicht direkt als Effekt auf die Wahrscheinlichkeit $P(Y=1\vert X)$ interpretiert werden.Um den partiellen Effekt einer Änderung in $X$ am Punkt $X$ auf $P(Y=1\vert X)$ zu ermitteln, berechnen wir die Ableitung des bedingten Erwartungswerts:\begin{align*}\frac{\partial\textup{E}(Y\vert X)}{\partial X} = \frac{\partial\textup{P}(Y=1\vert X)}{\partial X} = \frac{\partial\Lambda(\beta_0 + \beta_1 X)}{\partial X} = \lambda(\beta_0 + \beta_1 X) \beta_1,\end{align*}wobei $\lambda(\cdot)$, ähnlich wie die Dichtefunktion der Normalverteilung im Probit-Modell, die Dichtefunktion der logistischen Verteilung darstellt. Diese ist gegeben durch\begin{align*}\lambda(z) = \Lambda(z) \cdot (1 - \Lambda(z)).\end{align*}Die Interpretation des angepassten Modells ist aufgrund der Modellierung des Logits intuitiver als für Probit-Regression. Angenommen eine Bank möchte die Wahrscheinlichkeit modellieren, dass ein Kunde einen Kredit nicht zurückzahlt: $P(Y=\textup{Zahlungsausfall}\vert X)$, wobei der Regressor $X$ das Verhältnis von aktuellem Schuldenstand und Monatseinkommen ist. Die Schätzung einer logistischen Regression ergibt, dass ein Anstieg von $X =x$ um 0.1 den Logit für die Wahrscheinlichkeit eines Kreditausfalls $p(x)$ um 0.4 erhöht:\begin{align*} \textup{logit} \big[P(Y=\textup{Zahlungsausfall}\vert X = x + 0.1)\big] = \textup{logit}[p(x)] + 0.4. \end{align*}Das Modell besagt dann, dass die Wahrscheinlichkeit der Zahlungsunfähigkeit etwa um den *Faktor* 1.5 ansteigt, denn\begin{align*} \exp\big(\textup{logit}[p(x)] + 0.4\big) = \frac{p(x)}{1-p(x)} \cdot \exp(0.4)\end{align*}mit\begin{align*} \exp(0.4) \approx 1.50.\end{align*}### SchätzungÄhnlich wie bei Probit-Regressionen können Logit-Modelle mit Maximum-Likelihood geschätzt werden. Die Likelihood-Funktion für Logit-Regression lautet\begin{align} L(\boldsymbol{\beta}) &= \prod_{i=1}^n \left(\frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}\right)^{y_i} \left(1 - \frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}\right)^{1-y_i},\end{align}mit $$\frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}$$ der Wahrscheinlichkeit, dass $Y_i = 1$. Für das Ereignis $Y_i = 0$ ist die Wahrscheinlichkeit entsprechend $$1 - \frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}.$$Die Log-Likelihood-Funktion lautet\begin{align*} \mathcal{L}(\boldsymbol{\beta}) &= \sum_{i=1}^n \left[ y_i \log \left(\frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}\right) + (1 - y_i) \log \left(1 - \frac{1}{1 + \exp(-\mathbf{X}_i^\top \boldsymbol{\beta})}\right) \right]\end{align*}und erlaubt die Schätzung von $\boldsymbol{\beta}$ durch Maximierung mit numerischen Methoden.Für Anwendungen mit R nutzen wir `stats::glm()` und wählen die logistische Link-Funktion mit `family = binomial(link = "logit")`. Für die simulierten Daten `simdata` aus @sec-probitreg:```{webr}# Logit-Modell schätzenmod_logit <- glm( formula = Y ~ X, data = simdata, family = binomial(link = "logit"))# Robuste Zusammenfassungmod_logit %>% coeftest(vcov = vcovHC, type = "HC1") tidy()```Wir erweitern nun `pred` um die anhand von `predict()` für `mod_logit` geschätzten Wahrscheinlichkeiten von $Y=1\vert X$ für die Regressorwerte in `X`.```{webr}# gesch. WSK-Funktion gem. Logit-Modell hinzufügenpred <- pred %>% mutate( Logit = predict(mod_logit, tibble(X), type = "response") )```Anhand der Vorhersagen in `pred` können wir die im Objekt `p` gespeicherte Grafik um die geschätzte Regressionsfunktion für das Logit-Modell erweitern.```{webr}# Grafik mit Logit-fit erweiternp + geom_line( data = pred, mapping = aes(y = Logit, color = "Logit"), lwd = .75 ) + scale_color_manual( name = "", values = c( "Wahrheit" = "black", "LPM" = "steelblue", "Probit" = "orange", "Logit" = "darkred" ) )```Die Logit-Regression kann den wahren Zusammenhang gut abbilden und ist praktisch nicht von der Probit-Schätzung unterscheidbar. ### ModellgüteÄhnlich wie bei linearer Regression kann die Anpassung von generalisierten linearen Modellen an die Daten mit Anpassungsmaßen verglichen werden. In generalisierten linearen Modellen für binäre Outcome-Variablen kann hierzu die *Deviance* verwendet werden. Für eine Schätzung anhand der Beobachtungen $i=1,\dots,n$ berechnet man die Deviance $D$ als\begin{align*} D = -2 \sum_{i=1}^{n} \left[ y_i \cdot \log\left(\widehat{P}_i\right) + (1 - y_i) \cdot \log\left(1 - \widehat{P}_i\right) \right].\end{align*}$D$ quantifiziert die Abweichung eines geschätzten Modells von einem perfekt passenden Modell.^[Im Allgemeinen vergleicht Deviance die Log-Likelihood des geschätzten Modells mit der Log-Likelihood des perfekten Modells.] Eine geringere Deviance deutet darauf hin, dass das Modell die Daten besser erklärt. Eine R-Funktion zur Anwendung für Objekte des Typs `glm` ist `deviance()`.```{webr}# Deviance berechnenlist( "Probit" = mod_probit, "Logit" = mod_logit) %>% map_dbl(.f = deviance)```Die Deviance für das Probit-Modell ist tatsächlich etwas geringer als für das Logit-Modell. Ähnlich wie $R^2$ für lineare Regression misst die Deviance $D$ lediglich die Anpassung des Modells an die beobachteten Daten. Daher eignet sich $D$ nur bedingt als Maß zur Einschätzung der Tauglichkeit eines Modells für Out-of-sample-Vorhersagen. Abhängig vom Ziel der Modellierung kann es hilfreich sein, eine robuste Einschätzung der Vorhersage-Güte unter Berücksichtigung der Modellkomplexität vorzunehmen. Hierfür werden oft Informationskriterien verwendet. Das Akaike-Informationskriterium (AIC) \begin{align*} \text{AIC} = 2k - 2 \log(\widehat{\mathcal{L}})\end{align*}ist ein Maß zur Bewertung der Güte eines Modells unter Berücksichtigung der Modellkomplexität. Das AIC wägt die Anpassung des Modells an die Daten (gemessen durch die maximierte Likelihood-Funktion $\widehat{\mathcal{L}}$^[$\widehat{\mathcal{L}}$ meint die Likelihood-Funktion für die geschätzten Koeffizienten.]) und die Komplexität (gemessen als Funktion der Regressoranzahl $k$) gegeneinander ab. Hierbei werden Modelle mit weniger Parametern bevorzugt, wenn sie die Daten ähnlich gut erklären, sodass eine Überanpassung und damit eine schlechte Vorhersagefähigkeit vermieden werden. In R können wir das AIC für Modell-Objekte mit `AIC()` berechnen. ```{webr}# AIC berechnenlist( Probit = mod_probit, Logit = mod_logit) %>% map_dbl(.f = AIC)```Auch das AIC zeigt an, dass die Schätzung des Probit-Modells in `mod_probit` besser geeignet für die Modiellierung von `simdata` ist als das Logit-Modell `mod_logit`.### Beispiel: Klassifikation von Palmer-PiniguinenMit den in diesem Kapitel betrachteten Methoden können wir Modelle zur Klassifikation von Pinguinen in `palmerpenguins::penguins` hinsichtlich ihres Geschlechts konstruieren. Hierfür nutzen wir den bereinigten Datensatz `penguins_cleaned` und transformieren zunächst die abhängige Variable `sex` in ein numerisches Format.```{webr}# 'sex' transformierenpenguins_cleaned <- penguins_cleaned %>% mutate( # sex in 0,1-codierte Variable umwandeln sex = if_else(sex == "female", 1, 0) )# Überprüfenpenguins_cleaned %>% select(sex) %>% slice_head(n = 10)```Wir schätzen nachfolgend ein lineares Wahrscheinlichkeitsmodell `penguins_lp`, ein Probit-Modell `penguins_probit` sowie ein Logit-Modell `penguins_logit` die jeweils `sex` anhand der simplen Spezikation des linearen Prädiktors $$z = \beta_0 + \beta_1 \textup{bill\_depth}$$ erklären.```{webr}# Lineares Modellpenguins_lp <- lm( formula = sex ~ bill_depth, data = penguins_cleaned) # Robuste Zusammenfassungpenguins_lp %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()``````{webr}# Probit-Modellpenguins_probit <- glm( formula = sex ~ bill_depth, data = penguins_cleaned, family = binomial("probit"))# Robuste Zusammenfassungpenguins_probit %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()``````{webr}# Logit-Modellpenguins_logit <- glm( formula = sex ~ bill_depth, data = penguins_cleaned, family = binomial("logit"))# Robuste Zusammenfassungpenguins_logit %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()```Mit `modelsummary::modelsummary()` können wir die wichtigsten Ergebnisse der Regressionen in einer publikationsfähigen Tabelle darstellen. Hierzu übergeben wir dem Argument `models` eine (benannte) Liste der geschätzten Modelle. Bei Angabe von `vcov = "HC1"` werden robuste Standardfehler ausgegeben.```{webr}library(modelsummary)# Geschätzte Modellemodelsummary( models = list( "LPM" = penguins_lp, "Probit" = penguins_probit, "Logit" = penguins_logit ), output = "gt", vcov = "HC1", stars = TRUE, ) %>% tab_header( title = "Modelle für Geschlecht v. Pinguinen" )```Die geschätzten negativen und signifikanten Koeffizienten von $\textup{bill\_depth}$ erfassen den Zusammenhang, dass Pinguine mit größeren Schnäbeln tendenziell männlich sind: Je größer $\textup{bill\_depth}$, desto geringer die geschätzte Wahrscheinlichkeit, dass ein Pinguin weiblich ist. Ein Vergleich anhand des AIC zeigt, dass das Probit-Modell `penguins_probit` am besten geeignet ist.Analog zur Vorgehensweise in @sec-logreg können wir die geschätzten Modelle vergleichen, in dem wir geschätzte Wahrscheinlichkeiten $P(Y=\textup{weiblich}\vert\textup{bill\_depth})$ für eine Menge repräsentativer Werte des Regressors `bill_depth` berechnen und gemeinsam mit den Beobachtungen plotten.```{webr}# Sequenz für bill_depth erstellenbill_depth_new <- seq( from = min(penguins_cleaned$bill_depth), to = max(penguins_cleaned$bill_depth), by = .1)# Vorhersagen WSK mit 'bill_depth_new'preds <- tibble( bill_depth = bill_depth_new, lp = predict( object = penguins_lp, newdata = tibble(bill_depth = bill_depth_new) ), probit = predict( object = penguins_probit, newdata = tibble(bill_depth = bill_depth_new), type = "response" ), logit = predict( object = penguins_logit, newdata = tibble(bill_depth = bill_depth_new), type = "response" ))# Prüfen:slice_head(preds, n = 10)``````{webr}# Geschätzte bedingte WSK plottenpenguins_cleaned %>% ggplot( mapping = aes(x = bill_depth, y = sex) ) + geom_point( position = position_jitter( height = .025, seed = 1234 ), alpha = .25, color = "blue" ) + geom_line( data = preds, mapping = aes(y = lp, color = "LPM") ) + geom_line( data = preds, mapping = aes(y = probit, color = "Probit") ) + geom_line( data = preds, mapping = aes(y = logit, color = "Logit") ) + scale_color_manual( name = "Modelle", values = c( "LPM" = "steelblue", "Probit" = "orange", "Logit" = "darkred" ) ) + labs(title = "Modellierung v. Pinguin-Geschlecht mit Regression") + theme_cowplot() + theme(legend.position = "top")```Die geschätzten Regressionsfunktionen unterscheiden sich für den Bereich beobachteter Werte von $bill\_depth$ nur wenig. Die Grafik zeigt, dass die Schätzungen mit Probit- und Logit-Ansatz nur wenig Nicht-Linearität aufweisen, sodass eine (leichter interpretierbare) Modellierung mit dem linearen Wahrscheinlichkeitsmodell `penguins_lp` hier zulässig scheint.Die Klassifizierung der Pinguine im Datensatz hinsichtlich ihres Geschlechts anhand vorhergesagter Wahrscheinlichkeiten $\widehat{P}_i$ kann durch Abgleich mit einem Grenzwert erfolgen. Ein Beispiel ist\begin{align} \widehat{Y}_i = \begin{cases} \textup{weiblich}, & \widehat{P}_i \geq .5,\\ \textup{männlich}, & \textup{sonst.} \end{cases}\label{eq:pengclass}\end{align}Wir können diesen Ansatz in R implementieren, indem wir zunächst die für den Datensatz angepassten Wahrscheinlichkeiten `fitted` aus den Modell-Objekten auslesen.```{webr}( WSK <- penguins_cleaned %>% select(sex) %>% mutate( list( linear = penguins_lp, probit = penguins_probit, logit = penguins_logit ) %>% map_dfc(.f = ~ .$fitted) ))```Für die Klassifikation der Pinguine anhand der jeweiligen geschätzten Wahrscheinlichkeiten wenden wir die Regel \eqref{eq:pengclass} mit `across()` spaltenweise für die Modelle an und erzeugen neue Spalten mit vorhersagen des Geschlechts. Mit `.names = "{.col}_pred"` erhalten die neuen Spalten das Suffix `_pred`.```{webr}# Klassifikation f. Geschlecht durchführenvorh <- WSK %>% mutate( across( c(linear:logit), .fns = ~ (. > 0.5) == sex, .names = "{.col}_pred") )```Mit `summarise()` berechnen wir nun den Anteil korrekter Vorhersagen.```{webr}# Genauigkeit der Vorhersage berechnenvorh %>% summarise( across(linear_pred:logit_pred, mean) )```In diesem Beispiel unterscheiden sich die Vorhersagen der drei Modelle für den Grenzwert 0.5 nicht.## Modellierung von Zählvariablen {#sec-poissonreg}Eine weitere Klasse von $Y$ für die sich der Regressionsansatz von einfacher linearer Regression unterscheidet, sind *Zählvariablen*: Variablen, die diskrete, nicht-negative Werte annehmen. Wenn die abhängige Variable $Y$ Ereignisse in einem bestimmten Zeitraum misst (wie z.B. die Anzahl der Verkehrsunfälle in Essen innerhalb eines Monats) und weitere Bedingungen erfüllt sind, folgt $Y$ einer Poisson-Verteilung und kann mit Poisson-Regression modelliert werden. Wir erläutern nachfolgend die wesentlichen Komponenten von Poisson-Regression und diskutieren ein R-Beispiel mit fiktiven Daten.### Poisson-VerteilungDie Zufallsvariable $Y$ folgt einer Poisson-Verteilung mit Parameter $\lambda$, wenn1. **die Ereignisse unabhängig voneinander auftreten:** Das Auftreten eines Ereignisses beeinflusst nicht die Wahrscheinlichkeit, dass ein weiteres Ereignis auftritt.2. **die Ereignisse einzeln auftreten:** Die Wahrscheinlichkeit, dass in einem sehr kleinen Intervall (im Sinne von $[a,b]$ für $a \to b$) mehr als ein Ereignis auftritt, ist vernachlässigbar.3. **die Rate der Ereignisse konstant ist:** Die durchschnittliche Anzahl der Ereignisse pro Zeiteinheit oder pro Raumeinheit $\lambda$ ist konstant.Die Wahrscheinlichkeitsfunktion von $Y$ ist\begin{align}P(Y = y) = \frac{\lambda^y e^{-\lambda}}{y!} \quad \text{für} \quad y = 0, 1, 2, \ldots,\label{eq:poissonpmf}\end{align}wobei $\lambda$ sowohl der Erwartungswert als auch die Varianz der Verteilung ist, $$\textup{E}(Y) = \textup{Var}(Y) = \lambda.$$Die nächste Grafik zeigt die diskrete (Wahrscheinlichkeitsmasse-)Funktion \eqref{eq:poissonpmf} für eine Poisson-verteilte Zufallsvariable $Y$ mit $\lambda = 5$ für den Wertebereich von 0 bis 15.```{webr}#| label: fig-poissondensitytibble( x = 0:15, dens = dpois(x, lambda = 5)) %>% ggplot( mapping = aes(x = x, y = dens) ) + geom_bar( stat = "identity", fill = "steelblue", color = "white" ) + scale_x_continuous(breaks = 0:15) + labs( title = "Poisson-Dichtefunktion für lambda = 5", y = "Dichte" ) + theme_cowplot()```### Poisson-RegressionWie für die bisher betrachteten Regressionsansätzen modelliert Poisson-Regression den Erwartungswert (und damit gleichzeitig die Varianz) $\lambda$ der abhängigen Variable $Y$ als eine Funktion der Regressoren $\mathbf{X} = (X_1, X_2, \ldots, X_k)$. Unter Annahme einer Poisson-Verteilung kann Poisson-Regression als generalisiertes lineares Modell verstanden werden, wobei eine logarithmische Verknüpfungsfunktion für den linearen Prädiktor $\mathbf{X}_i^\top \boldsymbol{\beta}$ verwendet wird:\begin{align*}\log(\lambda_i) &= \mathbf{X}_i^\top \boldsymbol{\beta},\end{align*}sodass\begin{align*}\lambda_i &= \exp(\mathbf{X}_i^\top \boldsymbol{\beta}),\end{align*}wobei- $\lambda_i = \textup{E}(Y_i\vert \mathbf{X}_i) = \textup{Var}(Y_i\vert \mathbf{X}_i)$ für Beobachtung $i$,- $\mathbf{X}_i$ der Vektor der unabhängigen Variablen für Beobachtung $i$ ist, und- $\boldsymbol{\beta}$ der Vektor der Regressionskoeffizienten ist.### SchätzungDie Likelihood-Funktion für $n$ Beobachtungen ist \begin{align}L(\boldsymbol{\beta}) = \prod_{i=1}^n \frac{\lambda_i^{y_i} e^{-\lambda_i}}{y_i!}.\end{align}Die Log-Likelihood-Funktion ist daher\begin{align}\mathcal{L}(\boldsymbol{\beta}) = \sum_{i=1}^n \left( y_i \log(\lambda_i) - \lambda_i - \log(y_i!). \right)\end{align}Da $\lambda_i = \exp(\mathbf{X}_i^\top \boldsymbol{\beta})$, wird die Log-Likelihood-Funktion zu\begin{align}\mathcal{L}(\boldsymbol{\beta}) = \sum_{i=1}^n \left( y_i (\mathbf{X}_i^\top \boldsymbol{\beta}) - \exp(\mathbf{X}_i^\top \boldsymbol{\beta}) - \log(y_i!) \right)\end{align}Den ML-Schätzer $\widehat{\boldsymbol{\beta}}$ erhalten wir durch Maximierung der Log-Likelihoodfunktion $\mathcal{L}(\boldsymbol{\beta})$. Die R-Implementierung von Poisson-Regression in `stats::glm()` wird mittels `family = poisson(link = "log")`aufgerufen.### Interpretation der KoeffizientenDie Koeffizienten $\boldsymbol{\beta}$ in der Poisson-Regression haben eine logarithmisch-lineare Beziehung zur Zählvariable. Für einen bestimmten Koeffizienten $\beta_j$ ist die Interpretation wie folgt:Eine Änderung der unabhängigen Variable $X_j$ um eine Einheit führt zu einer Änderung des *Logarithmus* des Erwartungswertes von $Y$ um $\beta_j$. Der Erwartungswert $\lambda$ ändert sich also *multiplikativ* um den Faktor $\exp(\beta_j)$.Als Beispiel betrachten wir den linearen Prädiktor\begin{align*} \log(\lambda) = \beta_0 + \beta_1 X.\end{align*}Für eine Änderung von $X$ um $\Delta X$ erhalten wir\begin{align*} \lambda =&\, \exp(\beta_0 + \beta_1 X) \\ \lambda + \delta\lambda =&\, \exp(\beta_0 + \beta_1 (X + \Delta X)) \\ =&\, \exp(\beta_0 + \beta_1 X) \cdot \exp(\beta\Delta X) \\ =&\, \lambda \cdot \exp(\beta\Delta X)\end{align*}Angenommen $X$ ist die Anzahl durchgeführter Werbekampagnen und die abhängige Zählvariable $Y$ misst die Anzahl der Verkäufe des beworbenen Produkts pro Monat.Wenn $\beta_1 = 0.5$, bedeutet dies, dass jede zusätzliche Werbekampagne ($\Delta X = 1$) die erwartete Anzahl der Verkäufe pro Monat um einen Faktor von $\exp(0.5) \approx 1.65$ erhöht:Das heißt, die *Rate* der Verkäufe steigt um 65\% für jede zusätzliche Werbekampagne.Zur Veranschaulichung der Schätzung einer Poisson-Regression für dieses Beispiel erzeugen wir Poisson-verteilte Daten für\begin{align*} \log(\lambda) = 2 + 0.4 \cdot X\end{align*}und ziehen $X$ gleichverteilt aus $\{1,2,\dots,8\}$.```{webr}# Setze den Zufallszahlengenerator für Reproduzierbarkeitset.seed(1234)# Anzahl der Beobachtungenn <- 500# Züge für unabhängige Variable X (Anzahl der Werbekampagnen)X <- sample(1:8, replace = T, size = n)# Erwartungswertelambda <- exp(2 + 0.4 * X)# Abhängige Variable Y (Anzahl der Verkäufe) als Poisson-verteilte Zufallsvariable, gegeben lambdaY <- map_dbl( .x = 1:n, .f = ~ rpois( n = 1, lambda = lambda[.x] ))# Daten in tibble sammelndat <- tibble( Kampagnen = X, Verkaeufe = Y )# Überblickslice_head(dat, n = 10)```Die Züge aus der abhängigen Variablen visualieren wir mit einem Häufigkeits-Histogramm.```{webr}# Histogramm der Abh. Variableggplot( data = dat, mapping = aes(x = Verkaeufe) ) + geom_histogram(binwidth = 2) + labs( title = "Poisson-verteilte Daten" ) + theme_cowplot()```Wir schätzen das Modell und extrahieren die robuste Zusammenfassung mit `broom::tidy()`.```{webr}# Poisson-Regression schätzenmod_poisson <- glm( formula = Verkaeufe ~ Kampagnen, family = poisson(link = "log"), data = dat)# Zusammenfassung des gesch. Modellsmod_poisson %>% coeftest(vcov = vcovHC, type = "HC1") %>% tidy()```Die ML-Schätzung des Poisson-Modells liefert für die verwendete Stichprobe von 500 Beobachtungen Koeffizienten-Schätzungen nahe der wahren Parameter. Zur Veranschaulichung des geschätzten Modells überlagern wir die simulierten Datenpunkte aus `dat` mit der anhand von `mod_poisson` geschätzten Anzahl der Verkäufe je Anzahl an Werbekampagnen. ```{webr}# Vektor für KampagnenKampagnen_new <- tibble(Kampagnen = 1:8)# Vorhersagen mit Modellpred <- tibble( Kampagnen = 1:8, predicted = predict( mod_poisson, newdata = Kampagnen_new, type = "response" ))# Simulierte Daten und Schätzungenggplot() + geom_point( data = dat, mapping = aes( x = Kampagnen, y = Verkaeufe ), color = "gray", alpha = 0.25, position = position_jitter(width = .1) ) + geom_point( data = pred, mapping = aes( x = Kampagnen, y = predicted ), color = "red" ) + scale_x_continuous(breaks = 1:8) + labs( title = "Datensatz und Poisson-Schätzungen der Verkäufe", x = "Anzahl der Werbekampagnen", y = "Anzahl der Verkäufe" ) + theme_cowplot() + theme(legend.position = "top")```## ZusammenfassungIn diesem Kapitel haben wir erweiterte Konzepte der Regressionsanalyse diskutiert. Anhand des Frisch-Waugh-Lovell-Theorems wurde die Wirkungsweise der Kontrolle von Kovariablen mit multipler Regression auf den interessierenden Koeffizienten (Effekt) einer Variable veranschaulicht. Weiterhin haben wir gängige Typen generalisierter linearer Modelle für binäre und Poisson-verteilte Outcome-Variablen eingeführt. Anhand relevanter Beispiele wurde gezeigt, wie diese generalisierten Regressionsmodelle in R spezifiziert, geschätzt und der resultierende Output mit Paketen wie `broom`, `ggplot2` und `modelsummary` interpretiert werden kann.