[[5]{.chapter-number}  [Simulation]{.chapter-title}]{#sec-sim .quarto-section-identifier}

5 Simulation

Hinweis

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.

Bei der Analyse quantitativer Methoden treten häufig Fragestellungen auf, die sich rein mathematisch nur schwer oder gar nicht lösen lassen. So ist es sogar in einfachsten ökonometrischen Modellen oft unmöglich, die exakte Verteilung von Koeffizientenschätzern unter “realistischen” Bedingungen (etwa eine endliche Stichprobengröße und schwache Annahmen über die Verteilung der Fehlerterme) formal zu bestimmen. Asymptotische Theorie kann dann helfen, indem eine Approximation der Stichproben-Eigenschaften eines Schätzers aus Eigenschaften der asymptotischen Verteilung abgeleitet werden. Häufig sind solche Approximationen jedoch nur bei großen Stichproben verlässlich. Wichtige Stichprobeneigenschaften von Schätzern kausaler Effekte in empirischen Anwendungen, wie die Varianz oder das zu erwartende Ausmaß der Verzerrung bei Verletzungen von Modellannahmen hingegen können mit asymptotischer Theorie oftmals gar nicht approximiert werden.

Um diesen Schwierigkeiten zu begegnen, können wir Simulationstechniken (Monte-Carlo-Methoden) verwenden. In Monte-Carlo-Simulationen wird die Möglichkeit der Erzeugung von Pseudozufallszahlen mit Computern und ihre Rechenleistung genutzt, um eine Vielzahl von Datensätzen gemäß eines interessierenden datenerzeugenden Prozesses (DGP) zu generieren. Anhand dieser simulierten Datensätze können wir die Stichprobenverteilung von Schätzern approximieren und so die statistischen Eigenschaften für verschiedene Parameterwerte des DGP analysieren. Insbesondere erlauben Monte-Carlo-Studien den Vergleich verschiedener Schätzfunktionen hinsichtlich ihrer Robustheit gegenüber Verletzungen zugrundeliegender Annahmen über den DGP.

In diesem Kapitel erläutern wir Grundlagen in der Anwendung von Simulations-Methoden mit R und diskutieren Beispiele für die Schätzung kausaler Effekte mit Regressionsmodellen. Wir erläutern weiterhin das Konzept und die Anwendung des Bootstraps (Efron 1979), eine simulations-basierte Methode, die für viele Schätzer gültige Inferenz in empirischen Anwendungen erlaubt, wenn Modellannahmen derart verletzt sind, dass asymptotische kritische Werte (sogar für sehr große Stichproben) ungültig sind.

5.1 Simulation der Stichprobenverteilung eines Schätzers

Zur Illustration der Simulation von Stichprobenverteilungen eines Schätzers betrachten wir zunächst das kanonische Beispiel der Schätzung des Erwartungswerts $\mu_X$ einer normalverteilten Zufallsvariable $X$ anhand einer $n$-elementigen Stichprobe, \[\begin{align} X_i \sim\,u.i.v. N(\mu_X,\ \sigma_X^2), \quad i=1,\dots,n.\label{eq:distx} \end{align}\] Die Vorschrift $\eqref{eq:distx}$ ist der datenerzeugende Prozess (DGP)¹. Das arithmetische Mittel der $X_i$, \[\begin{align*} \overline{X} := \sum_{i=1}^n X_i, \end{align*}\] ist ein konsistenter Schätzer für $\mu_X$ und hat (aufgrund der Normalverteilung der $X_i$) die formal herleitbare Eigenschaft, dass die Stichprobenverteilung normal ist², \[\begin{align*} \overline{X} \sim N(\mu_X, \sigma_X^2 / n), \qquad \overline{X}\xrightarrow{p}\mu_X. \end{align*}\]

¹ Englische Abkürzung für data generating process.

² D. h. hier ist $\overline{X}\sim N(50, 2.88)$.

Wir können die Eigenschaften von $\overline{X}$ mit wenig Aufwand in R prüfen indem wir die Stichprobenentnahme durch wiederholtes Ziehen einer (großen) Anzahl von $N$ Stichproben aus der Populationsverteilung von $X$ nach der Vorschrift $\eqref{eq:distx}$ simulieren und jeweils $\overline{X}$ für diese simulierten Stichproben berechnen. Wir simulieren so Züge aus der Stichprobenverteilung von $\overline{X}$, deren Eigenschaften – in Abhängigkeit von $\mu_X$, $\sigma_X^2$ und $n$ – wir anhand der Realisierungen von $\overline{X}$ ableiten können. In den nachfolgenden Code-Beispielen wählen wir \[\begin{align*} \mu_X = 50,\quad \sigma_X^2=12^2,\quad n = 50. \end{align*}\]

Eine einfache Zufallsstichprobe von $X$ erzeugen wir mit rnorm() und berechnen das arithmetische Mittel mit mean().

Dieses einfache Simulationsexperiment mit $N=1$ ergibt zwar eine Schätzung des Erwartungswerts von $X$ nahe des wahren Werts $\mu_X = 50$, liefert jedoch wenig Information über die Qualität des Schätzers, d. h. die Eigenschaften der Stichprobenverteilung von $\overline{X}$. Diese Realisation könnte zufällig nahe bei $\mu_X = 50$ liegen, obwohl die Stichprobenverteilung $\overline{X}$ einen von $50$ verschiedenen Erwartungswert hat und/oder eine große Varianz aufweist.

Um die Stichprobenverteilung zu simulieren, wiederholen wir das Simulationsexperiment $N=1000$ mal: Wir erzeugen $N=1000$ Realisierungen aus der Stichprobenverteilung von $\overline{X}$ für die gewählten Parameter. Hierfür verwenden wir Iteration in R: Eine for()-Schleife wiederholt das Simulationsexperiment über 1:N Iterationen und speichert das Ergebnis des j-ten Experiments in einem zuvor definierten numerischen Vektor sim_est.³

³ Iteration mit Funktionen aus dem Paket purrr (siehe Tab Iteration mit purrr) ist eine Alternative zum “klassischen” Ansatz mit einer for()-Schleife.

Iteration mit for()-Schleife
Iteration mit purrr

Die zusammenfassenden Statistiken sind Schätzungen der Verteilungseigenschaften von $\overline{X}$. Da wir die Anzahl der Replikation $N$ selber wählen, können wir die wahren Parameter theoretisch “beliebig” genau approximieren: Lediglich die Rechenleistung des Computers limitiert die Güte dieser Approximation. In vielen Anwendungen liefert eine vierstellige Zahl an Simulationsdurchläufen gute Approximationen. Tatsächlich liegen die simulierten Parameter der Stichprobenverteilung von $\overline{X}$ nahe bei den wahren Werten.

Für einen grafischen Vergleich der simulierten Stichprobenverteilung mit der theoretischen $N(50, 12^2/50)$-Verteilung plotten wir die theoretische Dichte sowie ein Dichte-Histogramm der simulierten Züge von $\overline{X}$ in sim_est.

Die Abbildung zeigt, dass die simulierte Stichprobenverteilung (graues Dichte-Histogramm) die Dichte der theoretischen Verteilung (gestrichelte Linie) gut approximiert.

5.1.1 Simulation mit nicht-normalen Daten

Simulationen sind insbesondere dann hilfreich, wenn die Stichproben-Eigenschaften eines Schätzers unbekannt sind. Im nachfolgenden Beispiel wiederholen wir die Simulation für eine Log-normal-verteilte Zufallsvariable⁴ $X$ mit Erwartungswert $\mu_X = 50$ und plotten die geschätzten Dichtefunktionen der Stichprobenverteilungen von $\overline{X}$ für einen Vergleich.

⁴ Ist $X$ log-normal-verteilt mit Erwartungswert $\mu_X$ und Varianz $\sigma_X^2$, dann ist $\log X$ normal-verteilt mit $\mu_{\log X} = \log \mu_X - \sigma_X^2/2$ und $\sigma^2_{\log X} = \mu_X^2 [\exp(\sigma_X^2)-1]$.

Die funktionale Form der Stichprobenverteilung von $\overline{X}$ für eine log-normalverteilte Population ist unbekannt. Unsere Simulation erlaubt jedoch Rückschlüsse auf Eigenschaften dieser Verteilung: Die Verteilung ist rechtsschief und weist eine höhere Varianz auf als die Stichprobenverteilung von $\overline{X}$ für die normal-verteilte Population.

5.1.2 Regression: Verzerrung durch Confounder

Der DGP in Abbildung 5.1 zeigt eine typische Situation, in der wir den kausalen Effekt von $X$ auf $Y$ nicht ohne Weiteres mit Regression identifizieren können: Die Variablen $U$ und $e$ sind unbeobachtbar. $e$ ist ein unsystematischer Fehler, der lediglich $Y$ beeinflusst. Die Variable $U$ hingegen beeinflusst sowohl $X$ als auch $Y$ und ist damit eine Backdoor-Variable: $U$ ist ein Confounder, sodass der kausale Effekt von $X$ auf $Y$ in einer Regression, die nicht für den Effekt von $U$ kontrolliert, verzerrt geschätzt wird. Wie stark die Verzerrung ist, hängt von den Eigenschaften des DGP ab. Simulation kann hilfreich sein, um die Verzerrung eines Schätzers in Abhängigkeit der Parameter des DGP zu untersuchen.

Abbildung 5.1: DGP mit unbeobachtbarem Confounder U

Für die Simulation von Stichproben gemäß Abbildung 5.1 nutzen wir den linearen Prozess

\[\begin{align*} X =& \, a + \alpha_1 U,\\ Y =& \, \beta_1 X + \beta_2 U + e, \end{align*}\]

mit $a \sim \mathcal{U}[0,1]$, $U\sim N(0, 2^2)$ und $e\sim N(0, 10)$ und wählen die Parameter $\alpha_1 = 1$, $\beta_1 = 2$ und $\beta_2 = -1$.⁵ Der nächste Code-Block implementiert den DGP und zeigt die KQ-Schätzung des fehlspezifizierten Modells

⁵ $\mathcal{U}[a,b]$ meint die stetige Gleichverteilung auf dem Intervall $[a,b]$

\[\begin{align*} Y = \beta_1 X + \epsilon \end{align*}\]

anhand einer simulierten Stichprobe mit einem Umfang von 100 Beobachtungen.

Analog zu Kapitel 5.1 ist diese Simulation mit nur einer Iteration wenig aussagekräftig für die Verzerrung der KQ-Schätzers. Für eine ausführlichere Analyse implementieren wir zunächst beide Schritte der Simulation (Stichprobe generieren, Modell schätzen) in einer Funktion sim_function_conf(), der wir sämtliche Parameter sowie das zu schätzende Regressionsmodell als Fomel übergeben können. Diese Funktion gibt die aus dem lm-Objekt ausgelesene KQ-Schätzung von $\beta_1$ zurück. Wir testen die Funktion sim_function_conf() durch Verwendung des selben Seeds.

Der mit sim_function_conf() simulierte Zug aus der Stichprobenverteilung stimmt mit dem Ergebnis des vorherigen Zufallsexperiments überein. Wir können nun die Simulation durch iterative Berechnung von sim_function_conf() über die Indexmenge 1:N durchführen. Wie zuvor visualisieren wir die Simulationsergebnisse mit einem Plot der geschätzten Dichtefunktion von $\widehat{\beta}_1$.

Der grafische Vergleich der simulierten Stichprobenverteilung von $\widehat{\beta}_1$ mit dem wahren Koeffizienten $\beta_1 = 2$ (gestrichelte rote Linie) zeigt, dass $\widehat{\beta}_1$ aufgrund des ausgelassenen Confounders $U$ für die hier gewählte Spezifikation des DGP tatsächlich ein verzerrter Schätzer für den kausalen Effekt $\beta_1$ ist: Die geschätzte Dichtefunktion ist ein Indikator dafür, dass die Verteilungsfunktion von $\widehat{\beta}_1$ einen Großteil der Wahrscheinlichkeitsmasse im Wertebereich oberhalb von $\beta_1 = 2$ hat, sodass wir den kausalen Effekt im Mittel überschätzen.

Wir können die Simulationsergebnisse in sim_est nutzen, um die Verzerrung⁶ von $\widehat\beta_1$ zu schätzen.

⁶ Die Verzerrung (engl. bias) eines Schätzers $\widehat{\beta}_1$ für $\beta_1$ ist $\textup{bias}(\widehat{\beta}_1) = \textup{E}(\widehat{\beta}_1 - \beta_1)$.

Diese Berechnung anhand von Zügen aus der Stichprobenverteilung von $\widehat{\beta}_1$ zeigt eine deutliche positive Verzerrung des KQ-Schätzers von $\beta_1$ aufgrund des Confounders $U$ von etwa 0.94.

Anhand von sim_function_conf() können wir die Identifizierbarkeit des kausalen Effekts $\beta_1$ bei Kontrolle für $U$ im Modell

\[\begin{align*} Y = \beta_0 + \beta_1 X + \beta_2 U + \varepsilon \end{align*}\]

überprüfen, indem wir die Simulation unter Angabe dieses Modells über das Argument formula durchführen.

Die geschätzte Stichprobenverteilung von $\widehat{\beta}_1$ bei Kontrolle für den Confounder $U$ zeigt keine Anzeichen für eine Verzerrung des Schätzers für den kausalen Effekt von $X$ auf $Y$. Die Monte-Carlo-Schätzung der Verzerrung von $\widehat{\beta}_1$ unterstützt das Ergebnis der grafischen Analyse.

5.1.3 Regression: Verzerrung durch Collider

Der Graph in Abbildung 5.2 zeigt einen DGP, bei dem die Variable $C$ ein Collider auf dem Pfad $X\rightarrow C \leftarrow Y$ ist. Für die Identifikation des Effekts von $X$ auf $Y$ ist $C$ damit eine irrelavante Variable, für die nicht kontrolliert werden darf: Statistische Verfahren, die den Effekt von $X$ auf $Y$ unter Kontrolle des Effekts von $C$ auf $Y$ schätzen sind verzerrt, weil ein Backdoor-Pfad durch $C$ besteht.

Abbildung 5.2: DGP mit Collider-Variable C

Ähnlich wie in Kapitel 5.1.2 können wir die Konsequenzen der Kontrolle für eine Collider-Variable in einem linearen Regressionsmodell anhand einer Simulationsstudie untersuchen. Wir implementieren hierzu den DGP

\[\begin{align*} X \sim& \, N(0,1),\\ Y =& \, \beta_1 X + e,\\ C =& \, N(0,1) + .25X + .25Y, \end{align*}\] wobei $e \sim N(0,3^2)$ erneut ein unsystematischer Fehlerterm für $Y$ ist. Wie im DGP aus Kapitel 5.1.2 wählen wir $\beta_1 = 2$ für den wahren kausalen Effekt von $X$ auf $Y$.

In der nachfolgenden Simulation iterieren wir sim_function_col() für zwei verschiedene Modelle und speichern Ergebnisse in einer tibble: Die erste Spalte von sim_est_col enthält simulierte KQ-Schätzungen von $\widehat{\beta}_1$ in einem einfachen Regressionsmodell, dass nicht für den Collider $C$ kontrolliert (formula = "Y ~ X"). Die zweite Spalte hingegen sammelt Schätzungen für ein multiples Modell mit Kontrolle für $C$ (formula = "Y ~ X + C").

Für die Auswertung der Simulationsergebnisse hinsichtlich der Auswirkungen einer Kontrolle für den Collider $C$ plotten wir die geschätzen Dichtefunktionen für die Stichprobenverteilung von $\widehat{\beta}_1$ aus beiden Monte-Carlo-Experimenten gemeinsam und vergleich mit dem wahren Wert $\beta_1 = 2$.

Die Abbildung zeigt deutlich die (negative) Verzerrung des KQ-Schätzers für den interessierenden kausalen Effekt im Modell mit Kontrolle für den Collider $C$. Die Simulationsergebnisse für eine einfache Regression von $Y$ auf $X$ hingegen stimmen mit dem theretischen Resultat überein, dass der KQ-Schätzer für $\beta_1$ in diesem Modell erwartungstreu ist.

Durch Zusammenfassung der Simulationsergebnisse mit dplyr::summarize() können wir den Collider-Bias im Modell mit Kontrolle für $C$ schätzen und die Unverzerrtheit von $\widehat{\beta}_1$ in der einfachen Regression von $Y$ auf $X$ prüfen.

5.2 Regression: Durchschnittlicher Behandlungseffekt

Wir können Simulation einsetzen, um die Güte der Schätzung von Behandlungseffekten in linearen Regressionsmodellen zu untersuchen. Hierzu implementieren wir einen DGP gemäß Abbildung 5.3.

Abbildung 5.3: DGP mit heterogenen Behandlungseffekten

Für die Komponenten wählen wir: \[\begin{align} \begin{split} \text{B} &\sim \text{Bernoulli}(0.5) \\ X &\sim N(0, 2) \\ \beta_1 &\sim \mathcal{U}(0, 3) \\ e &\sim N(0, 2) \\ \\ Y &= 1 + \beta_1 \cdot \text{B} + 2 \cdot X + e \end{split}\label{eq:dgpatesim} \end{align}\]

In diesem linearen DGP bestimmt $B$ über die Zuteilung der Behandlung. $X$ ist eine Determinante der Outcome-Variable $Y$ ohne Backdoor-Pfad. Der Behandlungseffekt $\beta_1$ ist heterogen, da $\beta_1$ aus einer kontinuierlichen Gleichverteilung auf dem Intervall $[0, 3]$ gezogen wird. Der durchschnittliche Behandlungseffekt (ATE) entspricht hier dem Erwartungswert von $\beta_1$, was zu einem ATE von 1.5 führt, da

\[\begin{align*} \textup{ATE} =&\, \textup{E}(\beta_1)\\ =&\, \frac{1}{2}(a + b)\\ =&\, \frac{1}{2}(0 + 3) = 1.5. \end{align*}\]

Wir verifizieren nachfolgend mit Monte-Carlo-Simulation für Stichproben von DGP $\eqref{eq:dgpatesim}$ mit $n=150$, dass

der KQ-Schätzer für $\beta_1$ im Modell $Y=\beta_0 + \beta_1 B + \epsilon$ ein erwartungstreuer Schätzer für den ATE ist und
kontrollieren für $X$ hier die Präzision der Schätzung des ATE verbessert.

Im nachfolgenden Code-Block erzeugen wir die Simulationsergebnisse in jeder Iteration durch Berechnung der interessierenen Schätzer auf derselben simulierten Stichprobe und geben die Ergebnisse als tibble() mit einer Reihe und zwei benannten Einträgen zurück. Die Funktion map_dfr() fasst die Ergebnisse nach der Iteration über 1:N als tibble() zusammen. Wir nutzen hier einen Ansatz, bei dem die Schritte für die Simulation dem Argument .f über den Syntax ~ { ... } innerhalb einer anonymen Funktion ohne Argumente.⁷

⁷ Anonyme Funktionen in R (auch Lambda-Funktionen genannt) sind Funktionen, die ohne Namen definiert werden und in der Regel sofort angwendet werden. ~ ist eine Kurznotation für anonyme Funktionen innerhalb von purrr-Funktionen.

Für die grafische Auswertung transformieren wir die Simulationsergebnisse in sim_res_ATE mit pivot_longer() in ein langes Format.

Die grafische Analyse der Stichprobenverteilungen bestätigt, dass der KQ-Schätzer für $\beta_1$ sowohl im Modell mit Kontrolle für $X$ als auch im Modell ohne $X$ erwartungstreu für den ATE von 1.5 ist. Da $X$ gemäß DGP $\eqref{eq:dgpatesim}$ einen größeren Teil der Variation in $Y$ verursacht, verbessert die Kontrolle für $X$ hier die Präzision der Schätzung des ATE mit KQ.

Die Variabilität der Ansätze können wir durch Zusammenfassung der Spalten in sim_res_ATE mit summarize() und across() quantifizieren, wobei wir sd als zusammenfassende Statistik nutzen und .names = "sd_{.col}" ein entsprechendes Präfix für die Variablennamen einführt.

5.3 Verletzung von Modellannahmen: Heteroskedastizität

Unter der Annahme von Homoskedastizität (konstante Varianz) der Fehlerterme des Modells ist die KQ-Schätzung von Behandlungseffekten anhand linearer Regression effizient, d.h. der KQ-Schätzer hat die geringste Varianz unter allen unverzerrten Schätzern des interessierenden Koeffizienten.⁸ Hängt die Varianz der Fehlerterme von den Regressoren ab, liegt Heteroskedastizität vor: Die Fehlerterm-Varianz unterscheided sich zwischen den Beobachtungen, in Abhängigkeit der Ausprägungen der Regressoren. Unter Heteroskedastie ist der KQ-Schätzer nicht länger der effizienteste Schätzer: In Verfahren wie der gewichteten KQ-Methode (WKQ)⁹ haben Beobachtungen mit geringer Fehlerterm-Varianz einen größeren Einfluss als Beobachtungen mit großer Varianz, sodass der WKQ-Schätzer präziser als der KQ-Schätzer sein kann.

⁸ Die ist die Kern-Aussage des Gauss-Markov-Theorems.

⁹ Engl. weighted least squares.

Weiterhin sind bei heteroskedastischen Fehlertermen die (aus historischen Gründen) standardmäßig von summary() berechneten Standardfehler ungültig, und können die Unsicherheit von KQ unterschätzen. Dies ist problematisch, da unter Verwendung von ungültigen Standardfehlern berechnete Inferenzsstatistiken ebenfalls ungültig sind. Als Konsequenz erhöht sich die Gefahr falscher Schlussfolgerungen hinsichtlich des zu schätzenden Behandlungseffekts anhand von Teststatistiken und Konfidenzintervallen.

Wir zeigen nachfolgend, wie die Konsequenzen von Heteroskedastizität für den KQ-Schätzer mit Simulation untersucht werden können. Hierfür betrachten wir den DGP

\[\begin{align} \begin{split} X \sim&\, \chi^2_{50}\\ e \sim&\, N(0, \sigma_e^2(X))\\ Y =&\, \beta_1 X + e, \end{split}\label{eq:simhetdgp} \end{align}\]

wobei die Fehlerterm-Varianz $\sigma_e^2(X)$ eine Funktion von $X$ ist, \[\begin{align*} \sigma_e^2(X) = \textup{Var}(e\vert X) = \lvert X-50 \rvert^p. \end{align*}\] Der Parameter $p$ bestimmt die funktionale Form dieser bedingten Varianz.

Der nächste Code-Chunk plottet die Funktion $f(x) = \lvert X-50 \rvert^p$ für verschiedene $p$ über den Wertebereich von $X$. Hierfür berechnen wir $f(x)$ für ein mit expand_grid() erstelltes Raster von Werten für $X$ und $p$.

Für die Simulation heteroskedastischer Daten mit $\eqref{eq:simhetdgp}$ verwenden wir $p = 1.5$ und $\beta_1 = 0$ (kein Effekt). Die nachfolgende Grafik zeigt eine simulierte Stichprobe mit 150 Beobachtungen.

In den Simulationsdurchläufen schätzen wir jeweils $\beta_1$ im Modell $Y=\beta_1X+\epsilon$ mit KQ, berechnen den nur bei Homoskedastie gültigen Standardfehler und eine bei Heteroskedastie robuste Schätzung. Für lezteres verwenden wir sandwich::vcovHC() und mit type = "HC1".¹⁰ Wir fassen diese Schritte in der Funktion sim_function_het() zusammen.

¹⁰ Der Typ HC1, auch bekannt als Huber-White-Schätzer, ist ein häufig genutzter Standardfehler und wird in späteren Kapiteln dieses Companions häufig angewendet.

Die quantitative Auswertung der Ergebnisse deuten darauf hin, dass der nur bei Homoskedastie gültige Standardfehler die Standardabweichung von $\widehat{\beta}_1$ deutlich unterschätzt. Der anhand der Simulationsergebnisse geschätzte Erwartungswert des robusten Standardfehlers hingegen liegt nahe bei der Standardabweichung der simulierten Stichprobenverteilung von $\widehat{\beta}_1$. Die Tendenz des nur bei Homoskedastie gültigen Standardfehlers, kleinere Schätzungen der Unsicherheit von $\widehat{\beta}_1$ zu liefern, können wir anhand der geschätzten Dichtefunktionen für se und se_rob grafisch zeigen.

Als Konsequenz falscher Standardfehler kann die nominale $\alpha$-Fehlerrate von Signifikanztests verletzt sein: Der Hypothesentest lehnt dann eine korrekte Nullhypothese mit einer geringeren oder höheren Fehlerwahrscheinlichkeit (das nominale Signifikanzniveau $\alpha$) als gewünscht ab.

Inwiefern die in der Simulation gewählte Form von Heteroskedastizität die $\alpha$-Fehlerrate eines t-Tests beeinflusst, können wir anhand der Simulationsergebnisse schätzen: Wir berechnen die Ablehnrate für t-Tests der (wahren) Null-Hypothese $H_0:\beta_1=0$ zum 5%-Niveau jeweils für robuste und nicht-robuste Standardfehler und vergleichen. Wir vergleichen die so berechneten Test-Statistiken mit dem entsprechenden kritischen Wert der t-Verteilung mit $n-1$ Freiheitsgraden.

Die Simulationsergebnisse deuten darauf hin, dass ein Fehler erster Art für den Signifikanztest ohne korrigierte Standardfehler oberhalb des 5%-Niveaus liegt. Für die robuste Teststatistik hingegen liegt die Schätzung des $\alpha$-Fehlers nahe 5%.

5.4 Der Bootstrap

In den zuvor diskutierten Simulationsstudien haben wir iterativ Stichproben basierend auf einem bekannten DGP generiert, um die Eigenschaften der Stichprobenverteilung eines Schätzers (für bestimmte Parameterwerte) zu untersuchen. Der Bootstrap (Efron 1979) ist eine Methode zur Approximation der unbekannten Stichprobenverteilung eines Schätzers, bei der zufällige Stichproben mit Zurücklegen aus einer beobachteten Stichprobe gezogen werden.¹¹ Bootstrapping approximiert die Verteilung eines Schätzers aus beobachteten Datenpunkten und ohne Kenntnis des wahren DGP. Der Bootstrap ist daher eine Simulationsmethode, die in empirischen Anwendungen für die Bestimmung der Unsicherheit (Standardfehler) eines Schätzers genutzt werden kann. Bootstrapping ist hilfreich für Hypothesentests, wenn die beobachtete Stichprobe klein oder untypisch ist, was dazu führen kann, dass Standard-Inferenz für Behandlungseffekte nicht zuverlässig anwendbar oder ungültig ist. Insbesondere wenn Modellannahmen verletzt sind wird Bootstrapping häufig in empirischen Studien eingesetzt, um Konfidenzintervalle, Standardfehler oder p-Werte zu berechnen.

¹¹ Diese Vorgehensweise wird auch als Resampling bezeichnet.

Die nachfolgenden Key Facts fassen die wesentlichen Eigenschaften des Bootstraps und die grundsätzliche Vorgehensweise zusammen.

Key Facts zum Bootstrap

Der Bootstrap ist eine statistische Methode, die durch wiederholtes Ziehen von Stichproben mit Zurücklegen aus einer beobachteten Stichprobe die Stichprobenverteilung eines Schätzers approximieren kann. Bootstrapping ermöglicht die Berechnung gültiger Inferenzstatistiken, ohne starke Annahmen über den zugrundeliegenden unbekannten DGP zu treffen.

Vorgehensweise

Stichproben mit Zurücklegen ziehen: Ziehe wiederholt ($B$ mal) Stichproben der Größe $N$ (Größe der ursprünglichen Stichprobe) mit Zurücklegen aus der originalen Stichprobe.
Parameterschätzung für jede Bootstrap-Stichprobe: Wende die Schätzfunktion $\widehat{\beta}$ auf jede der $B$ Bootstrap-Stichproben an und erhalte $B$ Bootstrap-Schätzwerte $\widehat{\beta}^1,\dots,\widehat{\beta}^B$ des interessierenden Koeffizienten $\beta$.¹²
Verteilungsparameter Schätzen: Anhand der $B$ Bootstrap-Schätzungen aus Schritt 2 können die Eigenschaften der Verteilung des Schätzers $\widehat{\beta}$ (oder für Funktionen von $\widehat{\beta}$) geschätzt werden:
- Das arithmetische Mittel der $B$ Bootstrap-Schätzungen ist ein konsistenter Schätzer des Erwartungswerts von $\widehat{\beta}$
- Der Standardfehler des Schätzers $\widehat{\beta}$ kann durch die Standardabweichung der $B$ Bootstrap-Schätzungen approximiert werden
- Ein $1-\alpha\%$-Bootstrap-Konfidenzintervall für $\beta_1$ kann anhand des $\alpha/2\%$- und des $1-\alpha/2\%$-Perzentils der Bootstrap-Verteilung approximiert werden

Bootstrapping wird häufig genutzt, um die Unsicherheit eines Schätzers zu ermitteln. Dies ist insbesondere dann hilfreich, wenn keine analytische Darstellung für die Verteilung eines Schätzers existiert oder “starke” Verteilungsannahmen (z.B. Normalverteilung) unplausibel sind.
Ein Bootstrap-Schätzer für einen unbekannten Parameter¹³ ist im Allgemeinen nicht erwartungstreu für den interessierenden Effekt, selbst wenn ein erwartungstreuer Schätzer verwendet wird. Bootstrap-Schätzer sind jedoch in vielen Fällen unter schwachen Annahmen konsistent für den wahren Parameter.

¹³ Das arithmetische Mittel von $B$ Bootstrap-Schätzwerten ist ein einfacher Bootstrap-Schätzer.

¹² Häufig wird für $B$ eine ungerade Anzahl (z.B. 199) gewählt. Dies ermöglicht ein genaue Bestimmung des Medians der Bootstrap-Verteilung, was die Berechnung von Konfidenzintervallen erleichtert.

Der nächste Abschnitt diskutiert Bootstrapping der Stichprobenverteilung von $\overline{X}$ anhand des Beispiels aus Kapitel 5.1.

5.4.1 Interaktive Visualisierung des Bootstraps

Die nachfolgende interaktive Visualisierung illustriert die Vorgehensweise des oben beschriebenen Bootstrap-Algorithmus für die Schätzung eines $95\%$-Konfidenzintervalls für den Erwartungswert einer normalen Verteilung basierend auf einer Stichprobe, wenn das arithmetische Mittel als Schätzer verwendet wird. Diese Pupulationsverteilung ist $N(\mu = 50, \sigma^2 = 12^2)$. In jeder Iteration des Bootstraps wird eine Bootstrap-Stichprobe mit Zurücklegen aus der originalen Stichprobe gezogen. Für jede Bootstrap-Stichprobe berechnen wir das Arithmetische Mittel (grüne Punkte). Für hinreichend viele Bootstrap-Schätzungen erhalten wir eine repräsentative Bootstrap-Verteilung unseres Schätzers.¹⁴ Basierend auf dieser Verteilung berechnen wir eine Bootstrap-Schätzung des Erwartungswerts $\mu$ (schwarzer Punkt) und ein $95\%$-Konfidenzintervall (schwarze Linie).

¹⁴ Die Grüne Linie zeigt eine Kerndichte-Schätzung der Bootstrap-Verteilung. Die Bandweite dieser Schätzung kann in der Applikation variiert werden.

In der interaktiven Grafik können wir insbesondere folgende Eigenschaften des Bootstrap-Verfahrens überprüfen:

Die Populationsverteilung ist normal. Daher ist auch die Stichproben-Verteilung des arithmetischen Mittels $\overline{X}$ für jede Stichprobengröße normal. Die Illustration zeigt, dass die Bootstrap-Verteilung diese Normalverteilung des Schätzers (die typische Glockenform) tatsächlich gut approximieren kann, wenn $B$ hinreichend groß gewählt wird.
Der Bootstrap-Schätzer von $\mu$ ist nicht erwartungstreu für $\mu$, da wir mit Zurücklegen aus einer Stichprobe ziehen und so die zentrale Tendenz in den beobachteten Daten reproduzieren. Diese Verzerrung verringert sich mit der Größe der originalen Stichprobe.
Größere Stichproben verringern die Unsicherheit des interessierenden Schätzers $\overline{X}$: Die Varianz der Stichprobenverteilung von $\overline{X}$ nimmt mit der Stichprobengröße ab. Da der Bootstrap die Verteilung von $\overline{X}$ approximiert, spiegelt sich diese Eigenschaft auch in der Bootstrap-Verteilung.

5.4.2 Bootstrap mit R durchführen

Mit dem R-Paket boot kann Bootstrapping komfortabel für eine selbst definierte Schätzfunktion durchgeführt werden. Diese Funktion (boot_fun im nachfolgenden Code-Beispiel) muss den Schätzer unter Angabe einer Indexmenge für die beobachteten Daten berechnen und zurückgeben.¹⁵

¹⁵ boot::boot() zieht in jeder Bootstrap-Iteration mit Zurücklegen aus der Indexmenge der originalen Stichprobe, um die Bootstrap-Stichprobe festzulegen.

Der nächste Code-Chunk zeigt, wie boot::boot() angewendet werden kann, um den Bootstrap von Efron (1979) für $\overline{X}$ anhand einer Stichprobe mit $n=250$ aus der in Kapitel 5.1 betrachteten $N(\mu = 50, \sigma^2 = 12^2)$-Populationsverteilung für $B=999$ Iterationen zu berechnen.

Der Output in boot_est beinhaltet die Schätzung von $\mu$ mit $\overline{X}$ anhand der originalen Stichprobe (original), eine Schätzung der Verzerrung von $\overline{X}$ (bias) sowie den Bootstrap-Standardfehler (std. error). Wir können die Berechnung dieser Maße mit der originalen Stichprobe sowie den Bootstrap-Berechnungen von $\overline{X}$ (boot_est$t) nachvollziehen.

Die Bootstrap-Verteilung von $\overline{X}$ ist recht gut mit der tatsächlichen $N(\mu=50,\sigma^2 = 12^2/250)$-Stichprobenverteilung vergleichbar.

Mit boot::boot.ci() erhalten wir ein symmetrisches 95%-Bootstrap-Konfidenzintervall für $\mu$, \[\begin{align} 95\%\textup{-KI}_\textup{boot} = \big[2\cdot \overline{X}_o - q_{.975,\,\textup{boot}}, \ 2\cdot \overline{X}_o - q_{.025,\,\textup{boot}}\big].\label{eq:95bootci} \end{align}\]

Für die Berechnung per Hand mit der Formel $\eqref{eq:95bootci}$ bestimmen wir die Quantile der Bootstrap-Verteilung mit der Funktion quantile().

Anhand dieses Konfidenzintervalls können wir Hypothesen zum 5%-Niveau testen.

Wir können die (wahre) Nullhypothese $H_0: \mu = 50$ also zum 5%-Niveau nicht ablehnen.

5.5 Zusammenfassung

In diesem Kapitel haben wir Monte-Carlo-Simulationen als Hilfsmittel zur Analyse der Eigenschaften von Koeffizientenschätzern in ökonometrischen Modellen betrachtet. In Simulationen können die Parameter eines DGP vollständig kontrolliert und insbesondere variiert werden, um deren Auswirkungen auf das Verhalten eines Schätzers systematisch zu untersuchen: Simulation erlaubt es die Verteilung eines Schätzers durch iterative Berechnung auf simulierten Datensätzen näherungsweise zu bestimmen. Solche Approximationen sind hilfreich, da analytische Lösungen für die Verteilungseigenschaften von Schätzfunktionen oft nicht verfügbar sind. In methodischer Forschung erlauben Simulations-Studien es insbesondere, verschiedene anwendbare Schätzer gemäß ihrer Güte bei Verletzungen von Modellannahmen zu bewerten und so Empfehlungen hinsichtlich der Verlässlichkeit in empirischen Anwendungen zu treffen.

Der Bootstrap ist eine Simulationsmethode zur Abschätzung der Unsicherheit eines Schätzers im empirischen Anwenungen ohne starke Verteilungsannahmen zu treffen. Bootstrapping ist besonders nützlich, wenn es plausible ist, dass herkömmliche Inferenz unzuverlässig ist oder versagt, weil die Approximationen mit einer asymptotischen Verteilung aufgrund zu kleiner Stichprobengröße schlecht ist, oder gar die asymptotische Verteilung unbekannt ist. Der Bootstrap wird in empirischen Studien häufig zur Berechnung verlässlicher Konfidenzintervalle für kausale Effekte angewendet.

Die gezeigten Code-Beispiele illustrieren die Anwendung dieser Techniken in R mit Paketen für Iteration (purrr) und Resampling (boot).

--- format: live-html engine: knitr webr: packages: [ 'boot', 'cowplot', 'dplyr', 'ggplot2', 'purrr', 'sandwich', 'tidyr' ] --- {{< include ./_extensions/r-wasm/live/_knitr.qmd >}} # Simulation {#sec-sim} Bei der Analyse quantitativer Methoden treten häufig Fragestellungen auf, die sich rein mathematisch nur schwer oder gar nicht lösen lassen. So ist es sogar in einfachsten ökonometrischen Modellen oft unmöglich, die exakte Verteilung von Koeffizientenschätzern unter “realistischen” Bedingungen (etwa eine endliche Stichprobengröße und schwache Annahmen über die Verteilung der Fehlerterme) formal zu bestimmen. Asymptotische Theorie kann dann helfen, indem eine Approximation der Stichproben-Eigenschaften eines Schätzers aus Eigenschaften der asymptotischen Verteilung abgeleitet werden. Häufig sind solche Approximationen jedoch nur bei großen Stichproben verlässlich. Wichtige Stichprobeneigenschaften von Schätzern kausaler Effekte in empirischen Anwendungen, wie die Varianz oder das zu erwartende Ausmaß der Verzerrung bei Verletzungen von Modellannahmen hingegen können mit asymptotischer Theorie oftmals gar nicht approximiert werden. Um diesen Schwierigkeiten zu begegnen, können wir Simulationstechniken (Monte-Carlo-Methoden) verwenden. In Monte-Carlo-Simulationen wird die Möglichkeit der Erzeugung von Pseudozufallszahlen mit Computern und ihre Rechenleistung genutzt, um eine Vielzahl von Datensätzen gemäß eines interessierenden datenerzeugenden Prozesses (DGP) zu generieren. Anhand dieser simulierten Datensätze können wir die Stichprobenverteilung von Schätzern approximieren und so die statistischen Eigenschaften für verschiedene Parameterwerte des DGP analysieren. Insbesondere erlauben *Monte-Carlo-Studien* den Vergleich verschiedener Schätzfunktionen hinsichtlich ihrer Robustheit gegenüber Verletzungen zugrundeliegender Annahmen über den DGP. In diesem Kapitel erläutern wir Grundlagen in der Anwendung von Simulations-Methoden mit R und diskutieren Beispiele für die Schätzung kausaler Effekte mit Regressionsmodellen. Wir erläutern weiterhin das Konzept und die Anwendung des *Bootstraps* [@Efron1979], eine simulations-basierte Methode, die für viele Schätzer gültige Inferenz in empirischen Anwendungen erlaubt, wenn Modellannahmen derart verletzt sind, dass asymptotische kritische Werte (sogar für sehr große Stichproben) ungültig sind. ## Simulation der Stichprobenverteilung eines Schätzers {#sec-simspv} Zur Illustration der Simulation von Stichprobenverteilungen eines Schätzers betrachten wir zunächst das kanonische Beispiel der Schätzung des Erwartungswerts $\mu_X$ einer normalverteilten Zufallsvariable $X$ anhand einer $n$-elementigen Stichprobe, \begin{align} X_i \sim\,u.i.v. N(\mu_X,\ \sigma_X^2), \quad i=1,\dots,n.\label{eq:distx} \end{align} Die Vorschrift \eqref{eq:distx} ist der *datenerzeugende Prozess* (DGP)[^simulation-1]. Das arithmetische Mittel der $X_i$, \begin{align*} \overline{X} := \sum_{i=1}^n X_i, \end{align*} ist ein *konsistenter Schätzer* für $\mu_X$ und hat (aufgrund der Normalverteilung der $X_i$) die formal herleitbare Eigenschaft, dass die Stichprobenverteilung normal ist[^simulation-2], \begin{align*} \overline{X} \sim N(\mu_X, \sigma_X^2 / n), \qquad \overline{X}\xrightarrow{p}\mu_X. \end{align*} [^simulation-1]: Englische Abkürzung für *data generating process*. [^simulation-2]: D. h. hier ist $\overline{X}\sim N(50, 2.88)$. Wir können die Eigenschaften von $\overline{X}$ mit wenig Aufwand in R prüfen indem wir die Stichprobenentnahme durch wiederholtes Ziehen einer (großen) Anzahl von $N$ Stichproben aus der *Populationsverteilung* von $X$ nach der Vorschrift \eqref{eq:distx} simulieren und jeweils $\overline{X}$ für diese simulierten Stichproben berechnen. Wir simulieren so Züge aus der Stichprobenverteilung von $\overline{X}$, deren Eigenschaften -- in Abhängigkeit von $\mu_X$, $\sigma_X^2$ und $n$ -- wir anhand der Realisierungen von $\overline{X}$ ableiten können. In den nachfolgenden Code-Beispielen wählen wir \begin{align*} \mu_X = 50,\quad \sigma_X^2=12^2,\quad n = 50. \end{align*} Eine einfache Zufallsstichprobe von $X$ erzeugen wir mit `rnorm()` und berechnen das arithmetische Mittel mit `mean()`. ```{webr-r} #| autorun: true library(dplyr) # Seed für Reproduzierbarkeit set.seed(1234) # Parameter der Simulation n <- 50 mu_X <- 50 sigma_X <- 12 # Stichprobe ziehen, Mittel berechnen mean( rnorm( n = n, mean = mu_X, sd = sigma_X ) ) ``` Dieses einfache Simulationsexperiment mit $N=1$ ergibt zwar eine Schätzung des Erwartungswerts von $X$ nahe des wahren Werts $\mu_X = 50$, liefert jedoch wenig Information über die Qualität des Schätzers, d. h. die Eigenschaften der Stichprobenverteilung von $\overline{X}$. Diese *Realisation* könnte zufällig nahe bei $\mu_X = 50$ liegen, obwohl die Stichprobenverteilung $\overline{X}$ einen von $50$ verschiedenen Erwartungswert hat und/oder eine große Varianz aufweist. Um die Stichprobenverteilung zu simulieren, wiederholen wir das Simulationsexperiment $N=1000$ mal: Wir erzeugen $N=1000$ Realisierungen aus der Stichprobenverteilung von $\overline{X}$ für die gewählten Parameter. Hierfür verwenden wir *Iteration* in R: Eine `for()`-Schleife wiederholt das Simulationsexperiment über `1:N` Iterationen und speichert das Ergebnis des `j`-ten Experiments in einem zuvor definierten numerischen Vektor `sim_est`.[^simulation-3] [^simulation-3]: Iteration mit Funktionen aus dem Paket `purrr` (siehe Tab *Iteration mit purrr*) ist eine Alternative zum "klassischen" Ansatz mit einer `for()`-Schleife. ::: panel-tabset ## Iteration mit `for()-Schleife` ```{webr-r} #| autorun: true # Parameter für Iteration N <- 1000 sim_est <- numeric(N) set.seed(1234) # Simulation durchführen for (j in 1:N) { sim_est[j] <- mean( x = rnorm( n = n, mean = 50, sd = 12 ) ) } # Statistische Zusammenfassung tibble(sim_est) %>% summarise( mean_Xbar = mean(sim_est), sd_Xbar = sd(sim_est) ) ``` ## Iteration mit `purrr` ```{webr-r} #| autorun: true library(purrr) # Parameter N <- 1000 set.seed(1234) # Simulation mit purrr sim_est <- map_dbl( .x = 1:N, .f = ~ mean( x = rnorm( n = n, mean = 50, sd = 12 ) ) ) # Statistische Zusammenfassung tibble(sim_est) %>% summarise( mean_Xbar = mean(sim_est), sd_Xbar = sd(sim_est) ) ``` ::: Die zusammenfassenden Statistiken sind *Schätzungen* der Verteilungseigenschaften von $\overline{X}$. Da wir die Anzahl der Replikation $N$ selber wählen, können wir die wahren Parameter theoretisch "beliebig" genau approximieren: Lediglich die Rechenleistung des Computers limitiert die Güte dieser Approximation. In vielen Anwendungen liefert eine vierstellige Zahl an Simulationsdurchläufen gute Approximationen. Tatsächlich liegen die simulierten Parameter der Stichprobenverteilung von $\overline{X}$ nahe bei den wahren Werten. ```{webr-r} #| autorun: true # Simulierter Erwartungswert der Stichprobenverteilung mean(sim_est) # Theoretische Standardabweichung 12/sqrt(n) # Simulierte Standardabweichung der Stichprobenverteilung sd(sim_est) ``` Für einen grafischen Vergleich der simulierten Stichprobenverteilung mit der theoretischen $N(50, 12^2/50)$-Verteilung plotten wir die theoretische Dichte sowie ein Dichte-Histogramm der simulierten Züge von $\overline{X}$ in `sim_est`. ```{webr-r} #| autorun: true library(ggplot2) # Simulierte und theoretische Verteilungen des Schätzers plotten tibble(est = sim_est) %>% ggplot( mapping = aes(x = sim_est) ) + geom_histogram( mapping = aes(y = after_stat(density)), bins = 20, fill = "lightgray", color = "white" ) + stat_function( fun = dnorm, args = list( mean = 50, sd = 12/sqrt(n) ), color = "red", lty = 2 ) + labs(title = "Simulierte Verteilung des Stichprobenmittels") + theme_cowplot() ``` Die Abbildung zeigt, dass die simulierte Stichprobenverteilung (graues Dichte-Histogramm) die Dichte der theoretischen Verteilung (gestrichelte Linie) gut approximiert. ### Simulation mit nicht-normalen Daten Simulationen sind insbesondere dann hilfreich, wenn die Stichproben-Eigenschaften eines Schätzers *unbekannt* sind. Im nachfolgenden Beispiel wiederholen wir die Simulation für eine Log-normal-verteilte Zufallsvariable[^simulation-4] $X$ mit Erwartungswert $\mu_X = 50$ und plotten die geschätzten Dichtefunktionen der Stichprobenverteilungen von $\overline{X}$ für einen Vergleich. [^simulation-4]: Ist $X$ log-normal-verteilt mit Erwartungswert $\mu_X$ und Varianz $\sigma_X^2$, dann ist $\log X$ normal-verteilt mit $\mu_{\log X} = \log \mu_X - \sigma_X^2/2$ und $\sigma^2_{\log X} = \mu_X^2 [\exp(\sigma_X^2)-1]$. ```{webr-r} #| autorun: true library(tidyr) set.seed(1234) # Parameter n <- 50 N <- 1000 # Simulation durchführe, Ergebnisse sammeln df <- tibble( `Normal` = sim_est, `log-Normal` = map_dbl( .x = 1:N, .f = ~ mean( x = rlnorm( n = n, meanlog = log(50) - 1/2, sdlog = 1 ) ) ) ) # Gesch. Dichtefunktionen plotten df %>% pivot_longer( cols = everything(), names_to = "Population" ) %>% ggplot( mapping = aes( x = value, col = Population ) ) + geom_density() + theme_cowplot() + labs(title = "Stichprobenverteilung bei log-normaler Population") + theme(legend.position = "top") ``` Die funktionale Form der Stichprobenverteilung von $\overline{X}$ für eine log-normalverteilte Population ist unbekannt. Unsere Simulation erlaubt jedoch Rückschlüsse auf Eigenschaften dieser Verteilung: Die Verteilung ist rechtsschief und weist eine höhere Varianz auf als die Stichprobenverteilung von $\overline{X}$ für die normal-verteilte Population. ### Regression: Verzerrung durch Confounder {#sec-simbiasconf} Der DGP in @fig-simconf zeigt eine typische Situation, in der wir den kausalen Effekt von $X$ auf $Y$ nicht ohne Weiteres mit Regression identifizieren können: Die Variablen $U$ und $e$ sind unbeobachtbar. $e$ ist ein unsystematischer Fehler, der lediglich $Y$ beeinflusst. Die Variable $U$ hingegen beeinflusst sowohl $X$ als auch $Y$ und ist damit eine Backdoor-Variable: $U$ ist ein *Confounder*, sodass der kausale Effekt von $X$ auf $Y$ in einer Regression, die nicht für den Effekt von $U$ kontrolliert, verzerrt geschätzt wird. Wie stark die Verzerrung ist, hängt von den Eigenschaften des DGP ab. Simulation kann hilfreich sein, um die Verzerrung eines Schätzers in Abhängigkeit der Parameter des DGP zu untersuchen. ```{dot} //| fig-width: 6 //| fig-height: 4 //| fig-cap: "DGP mit unbeobachtbarem Confounder U" //| label: "fig-simconf" digraph CONFSIM { layout=neato fontname="Helvetica,Arial,sans-serif" node [fontname="Helvetica,Arial,sans-serif", shape="circle"] edge [fontname="Helvetica,Arial,sans-serif"] X [pos="0,0!"]; Y [pos="4,0!"]; U [pos="2,2!", color=gray, fontcolor=gray]; e [pos="6,0!", color=gray, fontcolor=gray]; U -> Y [color=gray, style=dashed, headlabel = "beta = -1", labeldistance=10, labelangle=-17]; U -> X [color=gray, style=dashed, headlabel = "beta = 1", labeldistance=10, labelangle=17]; X -> Y [headlabel = "beta = 2", labeldistance=12.5, labelangle=10]; e -> Y [color=gray, style=dashed] } ``` Für die Simulation von Stichproben gemäß @fig-simconf nutzen wir den linearen Prozess ```{=tex} \begin{align*} X =& \, a + \alpha_1 U,\\ Y =& \, \beta_1 X + \beta_2 U + e, \end{align*} ``` mit $a \sim \mathcal{U}[0,1]$, $U\sim N(0, 2^2)$ und $e\sim N(0, 10)$ und wählen die Parameter $\alpha_1 = 1$, $\beta_1 = 2$ und $\beta_2 = -1$.[^simulation-5] Der nächste Code-Block implementiert den DGP und zeigt die KQ-Schätzung des fehlspezifizierten Modells [^simulation-5]: $\mathcal{U}[a,b]$ meint die stetige [Gleichverteilung](https://de.wikipedia.org/wiki/Gleichverteilung) auf dem Intervall $[a,b]$ ```{=tex} \begin{align*} Y = \beta_1 X + \epsilon \end{align*} ``` anhand einer simulierten Stichprobe mit einem Umfang von 100 Beobachtungen. ```{webr-r} #| autorun: true # Parameter n <- 100 alpha_1 <- 1 beta_1 <- 2 beta_2 <- -.25 set.seed(1234) # Stichprobe ziehen dat <- tibble( U = rnorm(n, 0, 2), X = alpha_1 * U + runif(n, 0, 1), Y = beta_1 * X + beta_2 * U + rnorm(n, 10) ) # Fehlspezifiziertes Modell schätzen lm( formula = Y ~ X - 1, data = dat ) %>% summary() ``` Analog zu @sec-simspv ist diese Simulation mit nur einer Iteration wenig aussagekräftig für die Verzerrung der KQ-Schätzers. Für eine ausführlichere Analyse implementieren wir zunächst beide Schritte der Simulation (Stichprobe generieren, Modell schätzen) in einer Funktion `sim_function_conf()`, der wir sämtliche Parameter sowie das zu schätzende Regressionsmodell als Fomel übergeben können. Diese Funktion gibt die aus dem `lm`-Objekt ausgelesene KQ-Schätzung von $\beta_1$ zurück. Wir testen die Funktion `sim_function_conf()` durch Verwendung des selben Seeds. ```{webr-r} #| autorun: true # Funktion für Simulation definieren sim_function_conf <- function(n, alpha_1, beta_1, beta_2, formula) { # Simulieren tibble( U = rnorm(n, 0, 2), X = alpha_1 * U + runif(n, 0, 1), Y = beta_1 * X + beta_2 * U + rnorm(n, 10) ) %>% # KQ-Schätzung lm( formula = as.formula(formula), data = . ) %>% # Schätzwert auslesen coefficients() %>% .["X"] } # Seed wie zuvor set.seed(1234) # Test sim_function_conf( n = 100, alpha = 1, beta_1 = 2, beta_2 = -.25, formula = "Y ~ X - 1" ) ``` Der mit `sim_function_conf()` simulierte Zug aus der Stichprobenverteilung stimmt mit dem Ergebnis des vorherigen Zufallsexperiments überein. Wir können nun die Simulation durch iterative Berechnung von `sim_function_conf()` über die Indexmenge `1:N` durchführen. Wie zuvor visualisieren wir die Simulationsergebnisse mit einem Plot der geschätzten Dichtefunktion von $\widehat{\beta}_1$. ```{webr-r} #| autorun: true # Parameter n <- 100 N <- 1000 alpha_1 <- 1 beta_1 <- 2 beta_2 <- -.25 # Simulation durchführen sim_est <- tibble( hatbeta_1 = map_dbl( .x = 1:N, .f = ~ sim_function_conf( n = n, alpha_1 = alpha_1 , beta_1 = beta_1, beta_2 = beta_2 , formula = "Y ~ X - 1" ) ) ) # Gesch. Stichprobenverteilung plotten ggplot( data = sim_est, mapping = aes(x = hatbeta_1) ) + # gesch. Dichtefunktion geom_density() + # Effekt einzeichnen geom_vline( xintercept = beta_1, linetype = "dashed", col = "red" ) + labs(title = "Verzerrte Schätzung durch Confounder") + theme_cowplot() ``` Der grafische Vergleich der simulierten Stichprobenverteilung von $\widehat{\beta}_1$ mit dem wahren Koeffizienten $\beta_1 = 2$ (gestrichelte rote Linie) zeigt, dass $\widehat{\beta}_1$ aufgrund des ausgelassenen Confounders $U$ für die hier gewählte Spezifikation des DGP tatsächlich ein *verzerrter* Schätzer für den kausalen Effekt $\beta_1$ ist: Die geschätzte Dichtefunktion ist ein Indikator dafür, dass die Verteilungsfunktion von $\widehat{\beta}_1$ einen Großteil der Wahrscheinlichkeitsmasse im Wertebereich oberhalb von $\beta_1 = 2$ hat, sodass wir den kausalen Effekt im Mittel *überschätzen*. Wir können die Simulationsergebnisse in `sim_est` nutzen, um die Verzerrung[^simulation-6] von $\widehat\beta_1$ zu schätzen. [^simulation-6]: Die Verzerrung (engl. *bias*) eines Schätzers $\widehat{\beta}_1$ für $\beta_1$ ist $\textup{bias}(\widehat{\beta}_1) = \textup{E}(\widehat{\beta}_1 - \beta_1)$. ```{webr-r} #| autorun: true # Verzerrung des KQ-Schätzers schätzen mean(sim_est$hatbeta_1) - beta_1 ``` Diese Berechnung anhand von Zügen aus der Stichprobenverteilung von $\widehat{\beta}_1$ zeigt eine deutliche positive Verzerrung des KQ-Schätzers von $\beta_1$ aufgrund des Confounders $U$ von etwa 0.94. Anhand von `sim_function_conf()` können wir die Identifizierbarkeit des kausalen Effekts $\beta_1$ bei Kontrolle für $U$ im Modell ```{=tex} \begin{align*} Y = \beta_0 + \beta_1 X + \beta_2 U + \varepsilon \end{align*} ``` überprüfen, indem wir die Simulation unter Angabe dieses Modells über das Argument `formula` durchführen. ```{webr-r} #| autorun: true # Parameter n <- 100 N <- 1000 alpha_1 <- 1 beta_1 <- 2 beta_2 <- -.25 # Simulation durchführen sim_est <- tibble( hatbeta_1 = map_dbl( .x = 1:N, .f = ~ sim_function_conf( n = n, alpha_1 = alpha_1, beta_1 = beta_1, beta_2 = beta_2 , # Kontrolle für Confounder U formula = "Y ~ X + U" ) ) ) # Gesch. Stichprobenverteilung plotten ggplot( data = sim_est, mapping = aes(x = hatbeta_1) ) + geom_density() + geom_vline( xintercept = beta_1, linetype = "dashed", col = "red" ) + labs(title = "Erwartungstreue bei multipler Regression") + theme_cowplot() ``` Die geschätzte Stichprobenverteilung von $\widehat{\beta}_1$ bei Kontrolle für den Confounder $U$ zeigt keine Anzeichen für eine Verzerrung des Schätzers für den kausalen Effekt von $X$ auf $Y$. Die Monte-Carlo-Schätzung der Verzerrung von $\widehat{\beta}_1$ unterstützt das Ergebnis der grafischen Analyse. ```{webr-r} #| autorun: true # Verzerrung des KQ-Schätzers schätzen mean(sim_est$hatbeta_1) - beta_1 ``` ### Regression: Verzerrung durch Collider Der Graph in @fig-simcol zeigt einen DGP, bei dem die Variable $C$ ein *Collider* auf dem Pfad $X\rightarrow C \leftarrow Y$ ist. Für die Identifikation des Effekts von $X$ auf $Y$ ist $C$ damit eine irrelavante Variable, für die *nicht* kontrolliert werden darf: Statistische Verfahren, die den Effekt von $X$ auf $Y$ unter Kontrolle des Effekts von $C$ auf $Y$ schätzen sind verzerrt, weil ein Backdoor-Pfad durch $C$ besteht. ```{dot} //| fig-width: 6 //| fig-height: 4 //| fig-cap: "DGP mit Collider-Variable C" //| label: "fig-simcol" digraph COLSIM { layout=neato fontname="Helvetica,Arial,sans-serif" node [fontname="Helvetica,Arial,sans-serif", shape="circle"] edge [fontname="Helvetica,Arial,sans-serif"] X [pos="0,0!"]; Y [pos="4,0!"]; C [pos="2,2!"]; e [pos="6,0!", color=gray, fontcolor=gray]; Y -> C [headlabel = "beta = .25", labeldistance=10, labelangle=-17]; X -> C [headlabel = "beta = .25", labeldistance=10, labelangle=17]; X -> Y [headlabel = "beta_1 = 2", labeldistance=12.5, labelangle=10]; e -> Y [style=dashed, color=gray]; } ``` Ähnlich wie in @sec-simbiasconf können wir die Konsequenzen der Kontrolle für eine Collider-Variable in einem linearen Regressionsmodell anhand einer Simulationsstudie untersuchen. Wir implementieren hierzu den DGP \begin{align*} X \sim& \, N(0,1),\\ Y =& \, \beta_1 X + e,\\ C =& \, N(0,1) + .25X + .25Y, \end{align*} wobei $e \sim N(0,3^2)$ erneut ein unsystematischer Fehlerterm für $Y$ ist. Wie im DGP aus @sec-simbiasconf wählen wir $\beta_1 = 2$ für den wahren kausalen Effekt von $X$ auf $Y$. ```{webr-r} #| autorun: true # Simulationsfunktion: DGP mit Collider sim_function_col <- function(n, beta_1, formula) { # Daten erzeugen tibble( X = rnorm(n), Y = beta_1 * X + rnorm(n, sd = 3), C = .25 * X + .25 * Y + rnorm(n) ) %>% # Modell schätzen lm( formula = as.formula(formula), data = . ) %>% # Koeffizientenschätzung auslesen coefficients() %>% .["X"] } ``` In der nachfolgenden Simulation iterieren wir `sim_function_col()` für zwei verschiedene Modelle und speichern Ergebnisse in einer `tibble`: Die erste Spalte von `sim_est_col` enthält simulierte KQ-Schätzungen von $\widehat{\beta}_1$ in einem einfachen Regressionsmodell, dass *nicht* für den Collider $C$ kontrolliert (`formula = "Y ~ X"`). Die zweite Spalte hingegen sammelt Schätzungen für ein multiples Modell mit Kontrolle für $C$ (`formula = "Y ~ X + C"`). ```{webr-r} #| autorun: true # Parameter n <- 150 N <- 1000 beta_1 = 2 sim_est_col <- tibble( # Ohne Kontrolle für C hatbeta_1 = map_dbl( .x = 1:N, .f = ~ sim_function_col( n = n, beta_1 = beta_1, formula = "Y ~ X" ) ), # Kontrolle für C hatbeta_1_C = map_dbl( .x = 1:N, .f = ~ sim_function_col( n = n, beta_1 = beta_1, formula = "Y ~ X + C" ) ) ) # Überblick head(sim_est_col) ``` Für die Auswertung der Simulationsergebnisse hinsichtlich der Auswirkungen einer Kontrolle für den Collider $C$ plotten wir die geschätzen Dichtefunktionen für die Stichprobenverteilung von $\widehat{\beta}_1$ aus beiden Monte-Carlo-Experimenten gemeinsam und vergleich mit dem wahren Wert $\beta_1 = 2$. ```{webr-r} #| autorun: true # Vergleich von simulieren Stichprobenverteilungen # in Regressionen mit und ohne Collider ggplot(data = sim_est_col) + geom_density( mapping = aes(x = hatbeta_1, linetype = "Y ~ X"), ) + geom_density( mapping = aes(x = hatbeta_1_C, linetype = "Y ~ X + C"), ) + geom_vline( xintercept = beta_1, linetype = "dashed", col = "red" ) + scale_linetype_manual( "Modell", values = c( "Y ~ X" = "solid", "Y ~ X + C" = "dotdash" ) ) + labs(title = "Simulation von Collider-Bias in linearer Regression") + theme_cowplot() + theme(legend.position = "top") ``` Die Abbildung zeigt deutlich die (negative) Verzerrung des KQ-Schätzers für den interessierenden kausalen Effekt im Modell mit Kontrolle für den Collider $C$. Die Simulationsergebnisse für eine einfache Regression von $Y$ auf $X$ hingegen stimmen mit dem theretischen Resultat überein, dass der KQ-Schätzer für $\beta_1$ in diesem Modell erwartungstreu ist. Durch Zusammenfassung der Simulationsergebnisse mit `dplyr::summarize()` können wir den Collider-Bias im Modell mit Kontrolle für $C$ schätzen und die Unverzerrtheit von $\widehat{\beta}_1$ in der einfachen Regression von $Y$ auf $X$ prüfen. ```{webr-r} #| autorun: true # Bias aus Simulationsergebnissen schätzen sim_est_col %>% summarise( bias_simple = mean(hatbeta_1) - beta_1, bias_collider = mean(hatbeta_1_C) - beta_1 ) ``` ## Regression: Durchschnittlicher Behandlungseffekt Wir können Simulation einsetzen, um die Güte der Schätzung von Behandlungseffekten in linearen Regressionsmodellen zu untersuchen. Hierzu implementieren wir einen DGP gemäß @fig-simheteff. ```{dot} //| fig-width: 5 //| fig-height: 3 //| fig-cap: "DGP mit heterogenen Behandlungseffekten" //| label: "fig-simheteff" digraph HETEFFSIM { layout=neato fontname="Helvetica,Arial,sans-serif" node [fontname="Helvetica,Arial,sans-serif", shape="circle"] edge [fontname="Helvetica,Arial,sans-serif"] pad=.5 B [pos="0,0!"]; Y [pos="2,0!"]; X [pos="2,2!"]; e [pos="4,0!", color=gray, fontcolor=gray]; B -> Y [label = "β_1 heterogen"]; X -> Y [headlabel = "β_2 = 2", labeldistance=7, labelangle=-25]; e -> Y [style=dashed, color=gray]; } ``` Für die Komponenten wählen wir: \begin{align} \begin{split} \text{B} &\sim \text{Bernoulli}(0.5) \\ X &\sim N(0, 2) \\ \beta_1 &\sim \mathcal{U}(0, 3) \\ e &\sim N(0, 2) \\ \\ Y &= 1 + \beta_1 \cdot \text{B} + 2 \cdot X + e \end{split}\label{eq:dgpatesim} \end{align} In diesem linearen DGP bestimmt $B$ über die Zuteilung der Behandlung. $X$ ist eine Determinante der Outcome-Variable $Y$ ohne Backdoor-Pfad. Der Behandlungseffekt $\beta_1$ ist *heterogen*, da $\beta_1$ aus einer kontinuierlichen Gleichverteilung auf dem Intervall $[0, 3]$ gezogen wird. Der durchschnittliche Behandlungseffekt (ATE) entspricht hier dem *Erwartungswert* von $\beta_1$, was zu einem ATE von 1.5 führt, da ```{=tex} \begin{align*} \textup{ATE} =&\, \textup{E}(\beta_1)\\ =&\, \frac{1}{2}(a + b)\\ =&\, \frac{1}{2}(0 + 3) = 1.5. \end{align*} ``` Wir verifizieren nachfolgend mit Monte-Carlo-Simulation für Stichproben von DGP \eqref{eq:dgpatesim} mit $n=150$, dass - der KQ-Schätzer für $\beta_1$ im Modell $Y=\beta_0 + \beta_1 B + \epsilon$ ein erwartungstreuer Schätzer für den ATE ist und - kontrollieren für $X$ hier die Präzision der Schätzung des ATE verbessert. Im nachfolgenden Code-Block erzeugen wir die Simulationsergebnisse in jeder Iteration durch Berechnung der interessierenen Schätzer auf derselben simulierten Stichprobe und geben die Ergebnisse als `tibble()` mit einer Reihe und zwei benannten Einträgen zurück. Die Funktion `map_dfr()` fasst die Ergebnisse nach der Iteration über `1:N` als `tibble()` zusammen. Wir nutzen hier einen Ansatz, bei dem die Schritte für die Simulation dem Argument `.f` über den Syntax `~ { ... }` innerhalb einer *anonymen Funktion* ohne Argumente.[^simulation-7] [^simulation-7]: Anonyme Funktionen in R (auch *Lambda-Funktionen* genannt) sind Funktionen, die ohne Namen definiert werden und in der Regel sofort angwendet werden. `~` ist eine Kurznotation für anonyme Funktionen innerhalb von `purrr`-Funktionen. ```{webr-r} #| autorun: true # Parameter n <- 150 N <- 1000 set.seed(1234) # Iteration mit map_dfr sim_res_ATE <- map_dfr( .x = 1:N, .f = ~ { # Datensatz simulieren df <- tibble( # Behandlungsindikator treat = sample(0:1, size = n, replace = TRUE), # Determinante von Y X = rnorm(n = n, sd = 2), # Heterogener Behandlungseffekt TE = runif(n = n, 0, 3), # Outcome-Variable Y = 1 + TE * treat + 2 * X + rnorm(n = n, sd = 2) ) # ATE mit Regression schätzen `Modell ohne Kovariable X` = lm( formula = Y ~ treat, data = df ) %>% coefficients() %>% .["treat"] # Regression mit X `Modell mit Kovariable X` = lm( formula = Y ~ treat + X, data = df ) %>% coefficients() %>% .["treat"] return( # tibble zurückgeben tibble( `Modell ohne Kovariable X`, `Modell mit Kovariable X` ) ) } ) ``` Für die grafische Auswertung transformieren wir die Simulationsergebnisse in `sim_res_ATE` mit `pivot_longer()` in ein langes Format. ```{webr-r} #| autorun: true # Vergleich simulierter Dichtefunktionen sim_res_ATE %>% pivot_longer( cols = everything(), names_to = "variable", values_to = "hatbeta_1" ) %>% ggplot( mapping = aes( x = hatbeta_1, col = variable) ) + geom_density() + # ATE einzeichnen geom_vline( xintercept = 1.5, linetype = "dashed", col = "red" ) + scale_color_discrete("KQ-Schätzer im") + theme_cowplot() + labs(title = "Präzisere Schätzung mit mutlipler Regression") + theme(legend.position = "top") ``` Die grafische Analyse der Stichprobenverteilungen bestätigt, dass der KQ-Schätzer für $\beta_1$ sowohl im Modell mit Kontrolle für $X$ als auch im Modell ohne $X$ erwartungstreu für den ATE von 1.5 ist. Da $X$ gemäß DGP \eqref{eq:dgpatesim} einen größeren Teil der Variation in $Y$ verursacht, verbessert die Kontrolle für $X$ hier die Präzision der Schätzung des ATE mit KQ. Die Variabilität der Ansätze können wir durch Zusammenfassung der Spalten in `sim_res_ATE` mit `summarize()` und `across()` quantifizieren, wobei wir `sd` als zusammenfassende Statistik nutzen und `.names = "sd_{.col}"` ein entsprechendes Präfix für die Variablennamen einführt. ```{webr-r} #| autorun: true # Standardabweichung der Schätzungen sim_res_ATE %>% summarise( across( .cols = everything(), .fns = sd, .names = "sd_{.col}" ) ) ``` ## Verletzung von Modellannahmen: Heteroskedastizität Unter der Annahme von *Homoskedastizität* (konstante Varianz) der Fehlerterme des Modells ist die KQ-Schätzung von Behandlungseffekten anhand linearer Regression effizient, d.h. der KQ-Schätzer hat die geringste Varianz unter allen unverzerrten Schätzern des interessierenden Koeffizienten.[^simulation-8] Hängt die Varianz der Fehlerterme von den Regressoren ab, liegt *Heteroskedastizität* vor: Die Fehlerterm-Varianz unterscheided sich zwischen den Beobachtungen, in Abhängigkeit der Ausprägungen der Regressoren. Unter Heteroskedastie ist der KQ-Schätzer nicht länger der effizienteste Schätzer: In Verfahren wie der gewichteten KQ-Methode (WKQ)[^simulation-9] haben Beobachtungen mit geringer Fehlerterm-Varianz einen größeren Einfluss als Beobachtungen mit großer Varianz, sodass der WKQ-Schätzer präziser als der KQ-Schätzer sein kann. [^simulation-8]: Die ist die Kern-Aussage des [Gauss-Markov-Theorems](https://de.wikipedia.org/wiki/Satz_von_Gau%C3%9F-Markow). [^simulation-9]: Engl. *weighted least squares*. Weiterhin sind bei heteroskedastischen Fehlertermen die (aus historischen Gründen) standardmäßig von `summary()` berechneten Standardfehler *ungültig*, und können die Unsicherheit von KQ unterschätzen. Dies ist problematisch, da unter Verwendung von ungültigen Standardfehlern berechnete Inferenzsstatistiken ebenfalls ungültig sind. Als Konsequenz erhöht sich die Gefahr falscher Schlussfolgerungen hinsichtlich des zu schätzenden Behandlungseffekts anhand von Teststatistiken und Konfidenzintervallen. Wir zeigen nachfolgend, wie die Konsequenzen von Heteroskedastizität für den KQ-Schätzer mit Simulation untersucht werden können. Hierfür betrachten wir den DGP ```{=tex} \begin{align} \begin{split} X \sim&\, \chi^2_{50}\\ e \sim&\, N(0, \sigma_e^2(X))\\ Y =&\, \beta_1 X + e, \end{split}\label{eq:simhetdgp} \end{align} ``` wobei die Fehlerterm-Varianz $\sigma_e^2(X)$ eine Funktion von $X$ ist, \begin{align*} \sigma_e^2(X) = \textup{Var}(e\vert X) = \lvert X-50 \rvert^p. \end{align*} Der Parameter $p$ bestimmt die funktionale Form dieser bedingten Varianz. Der nächste Code-Chunk plottet die Funktion $f(x) = \lvert X-50 \rvert^p$ für verschiedene $p$ über den Wertebereich von $X$. Hierfür berechnen wir $f(x)$ für ein mit `expand_grid()` erstelltes Raster von Werten für $X$ und $p$. ```{webr-r} #| autorun: true # Raster erstellen, f(x) berechnen df <- expand_grid( x = seq(30, 70, by = 0.1), p = c(0, 1, 1.5, 1.75), ) %>% mutate(y = abs(x - 50)^p) # Funktionen plotten ggplot( data = df, mapping = aes(x = x, y = y, color = factor(p)) ) + geom_line() + labs( title = "Plot von |x-50|^p für verschiedene p", y = "|x-50|^p " ) + scale_color_discrete(name = "p") + theme_cowplot() + theme(legend.position = "top") ``` Für die Simulation heteroskedastischer Daten mit \eqref{eq:simhetdgp} verwenden wir $p = 1.5$ und $\beta_1 = 0$ (kein Effekt). Die nachfolgende Grafik zeigt eine simulierte Stichprobe mit 150 Beobachtungen. ```{webr-r} #| autorun: true n <- 150 # Daten erzeugen df_het <- tibble( # Regressor X = rchisq(n = n, df = 50), # Heteroskedastische Fehler e = abs(X - 50)^(1.5) * rnorm(n), # Outcome-Variable Y = e ) # Punkteplot ggplot( data = df_het, mapping = aes(x = X, y = Y) ) + geom_point(col = "steelblue", alpha = .75) + labs(title = "Heteroskedastische Daten") + theme_cowplot() ``` In den Simulationsdurchläufen schätzen wir jeweils $\beta_1$ im Modell $Y=\beta_1X+\epsilon$ mit KQ, berechnen den nur bei Homoskedastie gültigen Standardfehler und eine bei Heteroskedastie robuste Schätzung. Für lezteres verwenden wir `sandwich::vcovHC()` und mit `type = "HC1"`.[^simulation-10] Wir fassen diese Schritte in der Funktion `sim_function_het()` zusammen. [^simulation-10]: Der Typ `HC1`, auch bekannt als Huber-White-Schätzer, ist ein häufig genutzter Standardfehler und wird in späteren Kapiteln dieses Companions häufig angewendet. ```{webr-r} #| autorun: true sim_function_het <- function(n, p) { # Regressor X <- rchisq(n = n, df = 50) # Heteroskedastische Fehler e <- abs(X - 50)^(1.5) * rnorm(n, mean = 0, sd = 1) # Outcome-Variable Y <- e # Modell schätzen model <- lm(Y ~ X - 1) # HC0-Standardfehler mit summary() s <- summary(model) # Wurzel(rob Varianz-Schätzer) = robuster SE sw <- sqrt( sandwich::vcovHC(model, type = "HC1") ) # WKQ # Gewichte residuals <- resid(model) weights <- 1 / (residuals^2) # WKQ-Schätzung wkq_model <- lm( formula = Y ~ X - 1, weights = weights ) tibble( # Koeffizientenschätzung beta = s$coefficients["X", 1], # WKQ-Schätzer beta_wkq = wkq_model$coefficients["X"], # HC0-Standardfehler se = s$coefficients["X", 2], # HC1-Standardfehler se_rob = sw ) } ``` ```{webr-r} #| autorun: true # Parameter N <- 1000 n <- 150 p <- 1.5 set.seed(1234) # Simulation durchführen sim_res_het = map_dfr( .x = 1:N, .f = ~ sim_function_het(n = n, p) ) # Überblick head(sim_res_het) ``` ```{webr-r} #| autorun: true # Standardabweichung der Stichprobenverteilung sd(sim_res_het$beta) # Von summary() berechneter Standardfehler mean(sim_res_het$se) # Robuster Standardfehler mean(sim_res_het$se_rob) ``` Die quantitative Auswertung der Ergebnisse deuten darauf hin, dass der nur bei Homoskedastie gültige Standardfehler die Standardabweichung von $\widehat{\beta}_1$ deutlich unterschätzt. Der anhand der Simulationsergebnisse geschätzte Erwartungswert des robusten Standardfehlers hingegen liegt nahe bei der Standardabweichung der simulierten Stichprobenverteilung von $\widehat{\beta}_1$. Die Tendenz des nur bei Homoskedastie gültigen Standardfehlers, kleinere Schätzungen der Unsicherheit von $\widehat{\beta}_1$ zu liefern, können wir anhand der geschätzten Dichtefunktionen für `se` und `se_rob` grafisch zeigen. ```{webr-r} #| autorun: true ggplot(data = sim_res_het) + geom_density(mapping = aes(x = se, linetype = "Nicht robust")) + geom_density(mapping = aes(x = se_rob, linetype = "Robust")) + scale_linetype_manual( "Standardfehler", values = c("Nicht robust" = "solid", "Robust" = "dashed") ) + theme_cowplot() + labs(title = "Standardfehler bei Heteroskedastizität") + theme(legend.position = "top") ``` Als Konsequenz falscher Standardfehler kann die nominale $\alpha$-Fehlerrate von Signifikanztests verletzt sein: Der Hypothesentest lehnt dann eine korrekte Nullhypothese mit einer geringeren oder höheren Fehlerwahrscheinlichkeit (das nominale Signifikanzniveau $\alpha$) als gewünscht ab. Inwiefern die in der Simulation gewählte Form von Heteroskedastizität die $\alpha$-Fehlerrate eines t-Tests beeinflusst, können wir anhand der Simulationsergebnisse schätzen: Wir berechnen die Ablehnrate für t-Tests der (wahren) Null-Hypothese $H_0:\beta_1=0$ zum 5%-Niveau jeweils für robuste und nicht-robuste Standardfehler und vergleichen. Wir vergleichen die so berechneten Test-Statistiken mit dem entsprechenden kritischen Wert der t-Verteilung mit $n-1$ Freiheitsgraden. ```{webr-r} #| autorun: true # kritischer Wert der t_n-1-Verteilung tcrit <- qt(p = .975, df = n - 1) # Ablehnrate des t-Tests schätzen sim_res_het %>% summarise( # alpha-Fehler der Tests schätzen alpha_err_h = mean(abs(beta/se) > tcrit), alpha_err_r = mean(abs(beta/se_rob) > tcrit) ) ``` Die Simulationsergebnisse deuten darauf hin, dass ein Fehler erster Art für den Signifikanztest ohne korrigierte Standardfehler oberhalb des 5%-Niveaus liegt. Für die robuste Teststatistik hingegen liegt die Schätzung des $\alpha$-Fehlers nahe 5%. ## Der Bootstrap In den zuvor diskutierten Simulationsstudien haben wir iterativ Stichproben basierend auf einem bekannten DGP generiert, um die Eigenschaften der Stichprobenverteilung eines Schätzers (für bestimmte Parameterwerte) zu untersuchen. Der Bootstrap [@Efron1979] ist eine Methode zur Approximation der unbekannten Stichprobenverteilung eines Schätzers, bei der zufällige Stichproben mit Zurücklegen aus einer beobachteten Stichprobe gezogen werden.[^simulation-11] Bootstrapping approximiert die Verteilung eines Schätzers aus beobachteten Datenpunkten und *ohne Kenntnis des wahren DGP*. Der Bootstrap ist daher eine Simulationsmethode, die in empirischen Anwendungen für die Bestimmung der Unsicherheit (Standardfehler) eines Schätzers genutzt werden kann. Bootstrapping ist hilfreich für Hypothesentests, wenn die beobachtete Stichprobe klein oder untypisch ist, was dazu führen kann, dass Standard-Inferenz für Behandlungseffekte nicht zuverlässig anwendbar oder ungültig ist. Insbesondere wenn Modellannahmen verletzt sind wird Bootstrapping häufig in empirischen Studien eingesetzt, um Konfidenzintervalle, Standardfehler oder p-Werte zu berechnen. [^simulation-11]: Diese Vorgehensweise wird auch als *Resampling* bezeichnet. Die nachfolgenden Key Facts fassen die wesentlichen Eigenschaften des Bootstraps und die grundsätzliche Vorgehensweise zusammen. ::: callout-note #### Key Facts zum Bootstrap Der Bootstrap ist eine statistische Methode, die durch wiederholtes Ziehen von Stichproben mit Zurücklegen aus einer beobachteten Stichprobe die Stichprobenverteilung eines Schätzers approximieren kann. Bootstrapping ermöglicht die Berechnung gültiger Inferenzstatistiken, ohne starke Annahmen über den zugrundeliegenden unbekannten DGP zu treffen. **Vorgehensweise** 1. **Stichproben mit Zurücklegen ziehen**: Ziehe wiederholt ($B$ mal) Stichproben der Größe $N$ (Größe der ursprünglichen Stichprobe) mit Zurücklegen aus der originalen Stichprobe. 2. **Parameterschätzung für jede Bootstrap-Stichprobe**: Wende die Schätzfunktion $\widehat{\beta}$ auf jede der $B$ Bootstrap-Stichproben an und erhalte $B$ Bootstrap-Schätzwerte $\widehat{\beta}^1,\dots,\widehat{\beta}^B$ des interessierenden Koeffizienten $\beta$.[^simulation-12] 3. **Verteilungsparameter Schätzen**: Anhand der $B$ Bootstrap-Schätzungen aus Schritt 2 können die Eigenschaften der Verteilung des Schätzers $\widehat{\beta}$ (oder für Funktionen von $\widehat{\beta}$) geschätzt werden: - Das arithmetische Mittel der $B$ Bootstrap-Schätzungen ist ein konsistenter Schätzer des Erwartungswerts von $\widehat{\beta}$ - Der Standardfehler des Schätzers $\widehat{\beta}$ kann durch die Standardabweichung der $B$ Bootstrap-Schätzungen approximiert werden - Ein $1-\alpha\%$-Bootstrap-Konfidenzintervall für $\beta_1$ kann anhand des $\alpha/2\%$- und des $1-\alpha/2\%$-Perzentils der Bootstrap-Verteilung approximiert werden - Bootstrapping wird häufig genutzt, um die Unsicherheit eines Schätzers zu ermitteln. Dies ist insbesondere dann hilfreich, wenn keine analytische Darstellung für die Verteilung eines Schätzers existiert oder "starke" Verteilungsannahmen (z.B. Normalverteilung) unplausibel sind. - Ein Bootstrap-Schätzer für einen unbekannten Parameter[^simulation-13] ist im Allgemeinen *nicht erwartungstreu* für den interessierenden Effekt, selbst wenn ein erwartungstreuer Schätzer verwendet wird. Bootstrap-Schätzer sind jedoch in vielen Fällen unter schwachen Annahmen *konsistent* für den wahren Parameter. ::: [^simulation-12]: Häufig wird für $B$ eine ungerade Anzahl (z.B. 199) gewählt. Dies ermöglicht ein genaue Bestimmung des Medians der Bootstrap-Verteilung, was die Berechnung von Konfidenzintervallen erleichtert. [^simulation-13]: Das arithmetische Mittel von $B$ Bootstrap-Schätzwerten ist ein einfacher Bootstrap-Schätzer. Der nächste Abschnitt diskutiert Bootstrapping der Stichprobenverteilung von $\overline{X}$ anhand des Beispiels aus @sec-simspv. ### Interaktive Visualisierung des Bootstraps Die nachfolgende interaktive Visualisierung illustriert die Vorgehensweise des oben beschriebenen Bootstrap-Algorithmus für die Schätzung eines $95\%$-Konfidenzintervalls für den Erwartungswert einer normalen [Verteilung]{style="color: purple;"} basierend auf einer [Stichprobe]{style="color: red;"}, wenn das arithmetische Mittel als Schätzer verwendet wird. Diese [Pupulationsverteilung]{style="color: purple;"} ist $N(\mu = 50, \sigma^2 = 12^2)$. In jeder Iteration des Bootstraps wird eine [Bootstrap-Stichprobe]{style="color: orange;"} mit Zurücklegen aus der originalen [Stichprobe]{style="color: red;"} gezogen. Für jede [Bootstrap-Stichprobe]{style="color: orange;"} berechnen wir das [Arithmetische Mittel (grüne Punkte)]{style="color: darkgreen;"}. Für hinreichend viele [Bootstrap-Schätzungen]{style="color: darkgreen;"} erhalten wir eine repräsentative [Bootstrap-Verteilung]{style="color: darkgreen;"} unseres Schätzers.[^simulation-14] Basierend auf dieser Verteilung berechnen wir eine [Bootstrap-Schätzung des Erwartungswerts]{style="color: black;"} $\mu$ (schwarzer Punkt) und ein [$95\%$-Konfidenzintervall]{style="color: black;"} (schwarze Linie). [^simulation-14]: Die Grüne Linie zeigt eine Kerndichte-Schätzung der Bootstrap-Verteilung. Die Bandweite dieser Schätzung kann in der Applikation variiert werden. In der interaktiven Grafik können wir insbesondere folgende Eigenschaften des Bootstrap-Verfahrens überprüfen: - Die [Populationsverteilung]{style="color: purple;"} ist normal. Daher ist auch die Stichproben-Verteilung des arithmetischen Mittels $\overline{X}$ für *jede* Stichprobengröße normal. Die Illustration zeigt, dass die [Bootstrap-Verteilung]{style="color: darkgreen;"} diese Normalverteilung des Schätzers (die typische Glockenform) tatsächlich gut approximieren kann, wenn $B$ hinreichend groß gewählt wird. - Der Bootstrap-Schätzer von $\mu$ ist nicht erwartungstreu für $\mu$, da wir mit Zurücklegen aus einer [Stichprobe]{style="color: red;"} ziehen und so die zentrale Tendenz in den beobachteten Daten reproduzieren. Diese Verzerrung verringert sich mit der Größe der [originalen Stichprobe]{style="color: red;"}. - Größere [Stichproben]{style="color: red;"} verringern die Unsicherheit des interessierenden Schätzers $\overline{X}$: Die Varianz der Stichprobenverteilung von $\overline{X}$ nimmt mit der Stichprobengröße ab. Da der Bootstrap die Verteilung von $\overline{X}$ approximiert, spiegelt sich diese Eigenschaft auch in der [Bootstrap-Verteilung]{style="color: darkgreen;"}. <iframe class="obs-soft-box-shadow" width="100%" height="668" frameborder="0" src="https://observablehq.com/embed/62d24833f3dbd27a?cells=svg%2Cviewof+B%2Cviewof+BW%2Cviewof+speed%2Cviewof+n"> </iframe> ### Bootstrap mit R durchführen Mit dem R-Paket `boot` kann Bootstrapping komfortabel für eine selbst definierte Schätzfunktion durchgeführt werden. Diese Funktion (`boot_fun` im nachfolgenden Code-Beispiel) muss den Schätzer unter Angabe einer Indexmenge für die beobachteten Daten berechnen und zurückgeben.[^simulation-15] [^simulation-15]: `boot::boot()` zieht in jeder Bootstrap-Iteration mit Zurücklegen aus der Indexmenge der originalen Stichprobe, um die Bootstrap-Stichprobe festzulegen. Der nächste Code-Chunk zeigt, wie `boot::boot()` angewendet werden kann, um den Bootstrap von @Efron1979 für $\overline{X}$ anhand einer Stichprobe mit $n=250$ aus der in @sec-simspv betrachteten $N(\mu = 50, \sigma^2 = 12^2)$-Populationsverteilung für $B=999$ Iterationen zu berechnen. ```{webr-r} #| autorun: true library(boot) set.seed(1234) # Parameter festlegen n <- 250 B <- 999 # Stichprobe erzeugen stichprobe <- rnorm( n = n, mean = 50, sd = 12 ) # Schätzfunktion definieren boot_fun <- \(dat, indizes) mean(dat[indizes]) # Bootstrap durchführen ( boot_est <- boot( data = stichprobe, statistic = boot_fun, R = B ) ) ``` Der Output in `boot_est` beinhaltet die Schätzung von $\mu$ mit $\overline{X}$ anhand der originalen Stichprobe (`original`), eine Schätzung der Verzerrung von $\overline{X}$ (`bias`) sowie den Bootstrap-Standardfehler (`std. error`). Wir können die Berechnung dieser Maße mit der originalen Stichprobe sowie den Bootstrap-Berechnungen von $\overline{X}$ (`boot_est$t`) nachvollziehen. ```{webr-r} #| autorun: true # Stichprobenmittel (originale Stichprobe) mean(stichprobe) # äquiv. zu boot_est$t0 # Bootstrap-Schätzer für Erwartungswert mean(boot_est$t) # 'bias' mean(boot_est$t) - mean(stichprobe) # Bootstrap-Standardfehler sd(boot_est$t) ``` Die Bootstrap-Verteilung von $\overline{X}$ ist recht gut mit der tatsächlichen $N(\mu=50,\sigma^2 = 12^2/250)$-Stichprobenverteilung vergleichbar. ```{webr-r} #| autorun: true # Verteilung der Bootstrap-Mittelwerte plotten tibble(mu_boot = boot_est$t) %>% ggplot() + # Kerndichteschätzung der Bootstrap-Verteilung geom_density( mapping = aes(x = mu_boot) ) + # Wahre Stichprobenverteilung stat_function( fun = dnorm, args = list( mean = 50, sd = 12/sqrt(n) ), color = "red", lty = 2 ) + labs(title = "Vergleich von Bootstrap- und wahrer Stichprobenverteilung") + theme_cowplot() ``` Mit `boot::boot.ci()` erhalten wir ein symmetrisches 95%-Bootstrap-Konfidenzintervall für $\mu$, \begin{align} 95\%\textup{-KI}_\textup{boot} = \big[2\cdot \overline{X}_o - q_{.975,\,\textup{boot}}, \ 2\cdot \overline{X}_o - q_{.025,\,\textup{boot}}\big].\label{eq:95bootci} \end{align} ```{webr-r} #| autorun: true # 95%-Bootstrap-KI boot.ci( boot.out = boot_est, type = "basic" ) ``` Für die Berechnung per Hand mit der Formel \eqref{eq:95bootci} bestimmen wir die Quantile der Bootstrap-Verteilung mit der Funktion `quantile()`. ```{webr-r} #| autorun: true # 95%-Bootstrap-KI 'per Hand' berechnen ( KI_boot <- 2 * boot_est$t0 - quantile(boot_est$t, probs = c(.975, .025)) ) ``` Anhand dieses Konfidenzintervalls können wir Hypothesen zum 5%-Niveau testen. ```{webr-r} #| autorun: true # H_0: mu = 50 testen between(50, KI_boot[1], KI_boot[2]) ``` Wir können die (wahre) Nullhypothese $H_0: \mu = 50$ also zum 5%-Niveau nicht ablehnen. ## Zusammenfassung In diesem Kapitel haben wir Monte-Carlo-Simulationen als Hilfsmittel zur Analyse der Eigenschaften von Koeffizientenschätzern in ökonometrischen Modellen betrachtet. In Simulationen können die Parameter eines DGP vollständig kontrolliert und insbesondere variiert werden, um deren Auswirkungen auf das Verhalten eines Schätzers systematisch zu untersuchen: Simulation erlaubt es die Verteilung eines Schätzers durch iterative Berechnung auf simulierten Datensätzen *näherungsweise* zu bestimmen. Solche Approximationen sind hilfreich, da analytische Lösungen für die Verteilungseigenschaften von Schätzfunktionen oft nicht verfügbar sind. In methodischer Forschung erlauben Simulations-Studien es insbesondere, verschiedene anwendbare Schätzer gemäß ihrer Güte bei Verletzungen von Modellannahmen zu bewerten und so Empfehlungen hinsichtlich der Verlässlichkeit in empirischen Anwendungen zu treffen. Der *Bootstrap* ist eine Simulationsmethode zur Abschätzung der Unsicherheit eines Schätzers im empirischen Anwenungen ohne starke Verteilungsannahmen zu treffen. *Bootstrapping* ist besonders nützlich, wenn es plausible ist, dass herkömmliche Inferenz unzuverlässig ist oder versagt, weil die Approximationen mit einer asymptotischen Verteilung aufgrund zu kleiner Stichprobengröße schlecht ist, oder gar die asymptotische Verteilung unbekannt ist. Der *Bootstrap* wird in empirischen Studien häufig zur Berechnung verlässlicher Konfidenzintervalle für kausale Effekte angewendet. Die gezeigten Code-Beispiele illustrieren die Anwendung dieser Techniken in R mit Paketen für Iteration (`purrr`) und Resampling (`boot`).