6 Matching
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionWarning: package 'gt' was built under R version 4.3.3In randomisierten kontrollierten Studien stellt eine zufällige Zuweisung der Behandlung sicher, dass die Individuen beider Gruppen im Mittel vergleichbar sind, das heißt es gibt keine systematischen Unterschiede der Studiensubjekte hinsichtlich der Verteilung von Charakteristika zwischen den Gruppen. Dann ist es plausibel, eine beobachtete mittlere Differenz in der Outcome-Variable zwischen den Gruppen alleine auf die Behandlung zurückzuführen.
In der Praxis, insbesondere in ökonomischen Studien, sind randomisierte kontrollierte Studien aus ethischen und/oder finanziellen Gründen nicht durchführbar. Stattdessen werden nicht-experimentelle Daten1 genutzt, die jedoch nur selten eine unmittelbare Vergleichbarkeit von Behandlungs- und Kontrollgruppe gewährleisten.
1 Nicht-experimentelle Daten sind Beobachtungsdaten, bei denen die Zuordnung zu Behandlungs- und Kontrollgruppen nicht zufällig erfolgt, wodurch Confounding entstehen kann.
In diesem Kapitel betrachten wir Methoden, die in solchen Forschungsdesigns – bei hinreichenden Informationen über die Studiensubjekte – eine konsistente Schätzung kausaler Effekte ermöglichen können, indem eine statistische Vergleichbarkeit von Behandlungs- und Kontrollgruppe hergestellt wird. Dies kann durch eine gezielte Gewichtung von Beobachtungen anhand individueller Merkmale erfolgen. Andere etablierte Methoden schätzen den kausalen Effekt nach Selektion von vergleichbaren Teilmengen von Subjekten beider Gruppen aus der ursprünglichen Stichprobe, sogenanntes selektierendes Matching.
Da Matching Beobachtungen hinsichtlich beobachtbarer Merkmale vergleicht, kann die Wahrscheinlichkeit einer verzerrten Schätzung des kausalen Effekt durch falsche Modellspezifikationen geringer sein als für eine Schätzung des Effekts anhand multipler Regression. Weiterhin basieren Matching-Methoden nicht auf der Annahme eines linearen Zusammenhangs zwischen Kovariablen und der erklärenden Variable und können für die Schätzung unterschiedlicher Spezifikationen von Behandlungseffekten herangezogen werden.
Für die Gültigkeit eines Schätzers basierend auf Matching sind zwei Annahmen erforderlich.
Bedingte Unabhängigkeit. Seien \(Y^{(0)}_i\) und \(Y^{(1)}_i\) potentielle Ergebnisse der Outcome-Variable \(Y\) für ein Subjekt \(i\) mit \(B_i=0\) (keine Zuweisung zur Behandlung) beziehungsweise \(B_i=1\) (Behandlung) und \(X_i\) die beobachteten Kovariablen. Wir nehmen an, dass \[\begin{align} \left\{Y^{(0)}_i, Y^{(1)}_i\right\} \perp B_i\vert X_i, \label{eq:cia} \end{align}\] d.h. die Behandlungszuweisung/-selektion ist unabhängig von den potentielle Ergebnissen \(Y^{(0)}_i\) und \(Y^{(1)}_i\), wenn wir für die Kovariablen \(X\) kontrollieren.
Überlappung. Für jede mögliche Kombination von Kovariablen \(X_i\) gibt es eine positive Wahrscheinlichkeit \(<1\), sowohl zur Behandlungsgruppe (\(B_i = 1\)) als auch zur Kontrollgruppe (\(B_i = 0\)) zugewiesen zu werden, \[\begin{align} 0 < \text{P}(B_i=1\lvert X_i) < 1, \label{eq:overlap} \end{align}\] d.h. keine Beobachtung hat eine Behandlungswahrscheinlichkeit von exakt \(0\) oder \(1\).
Annahme 1 stellt sicher, dass die Zuweisung zur Behandlungsgruppe nach Kontrolle für die Kovariablen \(X\) als zufällig betrachtet werden kann. Somit ist es möglich den kausalen Effekt der Behandlung zu identifizieren, indem wir hinsichtlich der Kovariablen \(X\) ähnliche Subjekte (vgl. Annahme 2) aus Kontroll- und Behandlungsgruppe vergleichen.
Annahme 2 setzt voraus, dass es eine ausreichende Überlappung in den Verteilungen der Kovariablen zwischen Behandlungs- und Kontrollgruppe gibt. Dann ist sichergestellt, dass für jedes Subjekt in einer Gruppe ein hinsichtlich seiner Charakteristika vergleichbares Subjekt in der anderen Gruppe geben kann, sodass ein Vergleich möglich ist.
6.1 Balance: Vergleichbarkeit von Behandlungs- und Kontrollgruppe
Der Lehrstuhl für Ökonometrie an der Universität Duisburg-Essen betreibt einen Ökonometrie-Blog und interessiert sich für den kausalen Effekt der Einführung eines darkmode auf die Verweildauer der Leser auf der Webseite. Die Webseite ist zwar nicht-kommerziell, hat sich allerdings insbesondere für die Akquise internationaler Studierender für den Studiengang MSc. Econometrics als wichtiges Marketing-Instrument erwiesen. Ein ansprechendes Design wird daher als hoch-relevant erachtet.
Idealerweise sollte der Effekt des Design-Relaunches auf die Nutzungsintensität in einem kontrollierten randomisierten Experiment untersucht werden. Hierbei würden wir Nutzern zufällig das neue oder das alte Design zuweisen und den Effekt als Differenz des durchschnittlichen Verweildauer der Gruppen bestimmen. Eine solche Studie ist jedoch aus technischen und finanziellen Gründen nicht realisierbar, sodass die Auswirkungen des darkmode mit vorliegenden nicht-experimentellen Nutzungsstatistiken für die Webseite geschätzt werden sollen.
Die Nutzungsstatistiken sind im Datensatz darkmode.csv enthalten und sollen der Analyse des Effekts des darkmode (dark_mode) auf die Verweildauer der Leser auf der Webseite (read_time) dienen.
Tabelle 6.1 zeigt die Definitionen der Variablen in darkmode.csv.
| Variable | Beschreibung |
|---|---|
| read_time | Lesezeit (Minuten/Woche) |
| dark_mode | Indikator: Beobachtung nach Einführung darkmode |
| male | Indikator: Individuum männlich |
| age | Alter (in Jahren) |
| hours | Bisherige Verweildauer auf der Seite |
Für die Analyse lesen wir zunächst den Datensatz darkmode.csv mit readr::read_csv() ein und verschaffen uns einen Überblick über die verfügbaren Variablen.
dark_mode hat den Typ logical. Mit dplyr::mutate_all() können wir komfortabel alle Spalten in den Typ numeric transformieren.
Eine naive Schätzung des durchschnittlichen Behandlungseffekts (ATE) \(\widehat{\tau}^{\text{naiv}}\) erhalten wir als Mittelwertdifferenz von read_time für die Behandlungsgruppe (dark_mode == 1) und die Kontrollgruppe (dark_mode == 0) \[\begin{align}
\widehat{\tau}^{\text{naiv}} = \overline{\text{read\_time}}_{\text{Behandlung}} - \overline{\text{read\_time}}_{\text{Kontrolle}}.\label{eq:naivATEdarkmode}
\end{align}\]
Diese Berechnung ist schnell mit R durchgeführt.
Die Schätzung ergibt einen negativen Behandlungseffekt, mit der Interpretation, dass das neue Design zu einer Reduktion der Lesezeit um etwa 0.44 Minuten pro Woche führt. Dieses Ergebnis ist allerdings zweifelhaft, weil eine Isolierung des Behandlungseffekts aufgrund von Backdoor-Pfaden im DGP vermutlich nicht gewährleistet ist. Indikatoren hierfür sind systematische Unterschiede hinsichtlich von (möglicherweise unbeobachtbaren) Charakteristika von Kontrollgruppe und Behandlungsgruppe.
Da die User sich beim Aufrufen der Seite aktiv für oder gegen den das neue Design entscheiden müssen (und somit selektieren, ob Sie in der Behandlungs- oder Kontrollgruppe landen), liegt wahrscheinlich Confounding vor: Unsere Hypothese ist zunächst, dass männliche User eine durchschnittlich längere Lesezeit aufweisen und mit größerer Wahrscheinlichkeit auf das neue Design wechseln als nicht-männliche Leser. Dann ist male eine Backdoor-Variable. Diese Situation ist unter der Annahme, dass nur diese Faktoren den DGP bestimmen, in Abbildung 6.1 dargestellt.
Der DGP in Abbildung 6.1 führt zu einer verzerrten Schätzung des kausalen Effekts von dark_mode auf read_time mit \(\eqref{eq:naivATEdarkmode}\), wenn das Verhältnis von männlichen und nicht-männlichen Usern in Bahandlungs- und Kontrollgruppe nicht ausgeglichen ist. Wir überprüfen dies mit R.
Die Zusammenfassung anteile_m zeigt, dass der Anteil männlicher User in der Behandlungsgruppe deutlich höher ist als in der Kontrollgruppe.
6.1.1 Matching durch Gewichtung
Matching eliminiert die Variation von male zwischen den Gruppen. Eine Möglichkeit hierfür ist die Gewichtung der Beobachtungen in der Kontrollgruppe entsprechend der Anteile von Männern und Nicht-Männern in der Behandlungsgruppe, sodass die Vergleichbarkeit mit der Behandlungsgruppe hinsichtlich des Geschlechts gewährleistet ist. Dies wird in der Literatur als Balance bezeichnet. Der Behandlungseffekt wird dann analog zu \(\eqref{eq:naivATEdarkmode}\) geschätzt.
Die Gewichte für Beobachtungen in der Kontrollgruppe \(w_i\) werden berechnet als \[\begin{align}
w_i =
\begin{cases}
\text{ant\_m}_B/\text{anz\_m}_{K}, & \text{falls } \text{male}_i = 1\\
\text{ant\_nm}_B/\text{anz\_nm}_{K}, & \text{sonst.}\\
\end{cases}\label{eq:darkmodeweights}
\end{align}\] Anhand der Formel für einen gewichteten Durchschnitt, \[\begin{align}
\overline{X}_w = \frac{\sum_i w_i \cdot X_i}{\sum_i w_i},
\end{align}\] berechnen wir die gewichteten Mittelwerte für male und read_time in der Kontrollgruppe.
Ein Vergleich des gewichteten Mittelwertes von male in der Kontrollgruppe mit dem Mittelwert in der Behandlungsgruppe (male_k) zeigt, dass die Gewichte die Variation in male zwischen beiden Gruppen eliminieren, sodass die Backdoor durch male geschlossen ist. Mit wmean_read_time_K haben wir einen entsprechend gewichteten Mittelwert der Verweildauer für die Kontrollgruppe berechnet. Wir schätzen den Behandlungseffekt nun als \[\begin{align}
\widehat{\tau}^{\text{w}} = \overline{\text{read\_time}}_{B} - \overline{\text{read\_time}}_{w,K}.\label{eq:weightedATEdarkmode}
\end{align}\]
Entgegen der naiven Schätzung andhand von \(\eqref{eq:naivATEdarkmode}\) erhalten wir nach Matching für male eine positive Schätzung des Behandlungseffekts.
Die Schätzung des Behandlungseffekts anhand von \(\eqref{eq:weightedATEdarkmode}\) entspricht dem geschätzten Koeffizienten \(\widehat{\beta}_1\) aus einer gewichteten KQ-Regression im Modell \[\begin{align*} \text{read\_time} = \beta_0 + \beta_1 \text{dark\_mode} + u, \end{align*}\] wobei die Beobachtungen der Kontrollgruppe wie in \(\eqref{eq:darkmodeweights}\) gewichtet werden und \(w_i=1\) für Beobachtungen der Behandlungsgruppe ist. Wir überprüfen dies mit R.
Der geschätzte Koeffizient von dark_mode entspricht \(\widehat{\tau}^w\).
Da male eine binäre Variable ist, reduziert sich eine Beurteilung der Vergleichbarkeit der Verteilungen von male in Behandlungs- und Kontrollgruppe auf einen simplen Vergleich des Männeranteils beider Gruppen. In der Praxis gibt es oft eine Vielzahl potentieller Backdoor-Variablen, die zudem kontinuierlich verteilt sind. Es scheint plausibel, dass das Alter der Nutzer sowohl die Akzeptanz des Design-Updates als auch die Lesezeit beeinflusst. Die bisherige Verweildauer ist mindestens eine plausible Determinante der Lesezeit.
Der erweiterte DGP ist in Abbildung Abbildung 6.2 dargestellt, wobei ein zusätzlicher Backdoor-Pfad durch age verläuft.
Die Beurteilung der Balance von Kontrollgruppe und Behandlungsgruppe kann durch eine grafische Gegenüberstellung der empirischen Verteilungen der Kovariablen beider Gruppen erfolgen. Wir visualisieren die empirischen Verteilungen mit ggplot2. Hierzu standardisieren wir age und hours zunächst mit base::scale().
Für die Kovariablen age und hours eignen sich die geschätzten Dichtefunktionen für einen Vergleich der Verteilungen in Behandlungs- und Kontrollgruppe.
Die graphische Analyse zeigt deutliche Unterschiede in den Verteilungen von age zwischen Kontroll- und Behandlungsgruppe. Für eine Beurteilung mit deskriptiven Statistiken wird häufig eine sogenannte Balance Table herangezogen. Wir berechnen diese für age, hours und male mit cobalt::bal.tab()
Die Einträge M.0.Un und M.1.Un zeigen die jeweiligen Stichprobenmittelwerte der Variablen für Kontroll- und Behandlungsgruppe. Diff.Un gibt eine standardisierte Mittelwertdifferenz \(\textup{SMD}\) an, wobei \[\begin{align*}
\textup{SMD}_j := \left(\overline{X}_{j,B} - \overline{X}_{j,K}\right) \bigg/ \sqrt{\frac{1}{2}\left(\widehat{\text{Var}}(X_{j,B}) + \widehat{\text{Var}}(X_{j,K})\right)},
\end{align*}\] mit Stichprobenmitteln \(\overline{X}_{j,B}\) und \(\overline{X}_{j,K}\) und Stichprobenvarianzen \(\widehat{\text{Var}}(X_{j,B})\) und \(\widehat{\text{Var}}(X_{j,K})\) für eine kontinuierliche Kovariable \(j\).2 Obwohl es keinen einheitlichen Schwellenwert für die standardisierte Differenz gibt, der ein erhebliches Ungleichgewicht anzeigt, gilt für kontinuierliche Variablen eine standardisierte (absolute) Differenz von weniger als \(0.1\) als Hinweis auf einen vernachlässigbaren Unterschied zwischen den Gruppen.
2 Siehe P. Austin (2011) für einen Überblick zu Balance-Statistiken.
Die Balance Table weist also auf einen vernachlässigbaren Unterschied für hours hin und bestätigt den aus den Grafiken abgeleiteten Eindruck einer relevanten Differenzen für age.
6.1.2 Entropy Balancing
Entropy Balancing (Hainmueller 2012) ist eine weitere Methode zur Gewährleistung der Vergleichbarkeit von Behandlungs- und Kontrollgruppe anhand von Gewichten. Das Verfahren nutzt Konzepte aus der Informationstheorie um die Gewichte für Subjekte in der Kontrollgruppe so anzupassen, dass die Verteilung der Matchingvariablen in der Kontrollgruppe die Verteilung in der Behandlungsgruppe möglichst gut approximiert. Dies geschieht unter der Restriktion, dass bestimmte empirische Momente (meist Mittelwerte und Varianzen) der Matchingvariablen exakt übereinstimmen. Mathematisch werden die Gewichte für Kontrollgruppenbeobachtungen durch Minimierung der Kullback-Leibler-Divergenz zwischen den Verteilungen gefunden, wobei die Divergenz ein Maß für den Unterschied von Wahrscheinlichkeitsverteilungen ist.
Entropy Balancing ist im R-Paket WeightIt implementiert. Wir zeigen, wie die benötigten Gewichte für eine Schätzungen des ATT im Website-Beispiel mit WeightIt::wightit() bestimmt werden können. Über das Argument moments legen wir fest, dass die Gewichte unter der Restriktion übereinstimmender Mittelwerte aller Matching-Variablen zwischen Behandlungs- und Kontrollgruppe erfolgen soll.
Wir schätzen den Behandlungseffekt nach Entropy Balancing mit gewichteter Regression.
Beachte, dass die von summary() berechneten Standardfehler bei Entropy Balancing ungültig sind. In Abschnitt Kapitel 6.4 erläutern wir die Berechnung von Standardfehlern für Matching-Schätzer mit dem Bootstrap.
6.1.3 Mehrere Matching-Variablen und der Propensity Score
Bei mehreren Backdoor-Variablen kann eine Gewichtung anhand der Behandlungswahrscheinlichkeit (Treatment Propensity) erfolgen. Die Idee hierbei ist, dass der DGP wie in Abbildung 6.3 dargestellt werden kann.
Hierbei beeinflussen die Backdoor-Variablen age und male die Behandlungsvariable dark_mode lediglich durch die Behandlungswahrscheinlichkeit Treatment Propensity. Diese Darstellung zeigt, das die mehrdimensionale Information bzgl. der Ähnlichkeit von Subjekten hinsichtlich der beobachteten Kovariablen in einer einzigen Variable zusammengefasst werden kann. Die Backdoor-Pfade können daher geschlossen werden, indem wir Subjekte anhand von Treatment Propensity derart gewichten, dass beide Gruppen hinsichtlich der Verteilung der Behandlungswahrscheinlichkeit vergleichbar sind. Betrachte erneut \(\eqref{eq:cia}\) und beachte, dass \[\begin{align} Y_i = Y_i^{(1)} D_i + Y_i^{(0)} (1-D_i). \end{align}\] Rosenbaum und Rubin (1983) zeigen, dass es hinsichtlich \(\eqref{eq:cia}\) äquivalent ist für die Treatment Propensity \(P_i(X_i):=P(B_i=1\vert X_i = x)\) zu kontrollieren, d.h. \[\begin{align} \left\{Y_i^{(1)},Y_i^{(0)}\right\} \perp B_i\vert X_i \quad\Leftrightarrow\quad \left\{Y_i^{(1)},Y_i^{(0)}\right\} \perp B_i\vert P_i(X_i). \end{align}\]
Der Behandlungseffekt kann so als Differenz von gewichteten Gruppenmittelwerten berechnet werden, mit inversem Wahrscheinlichkeitsgewicht (IPW) \(w_{i,B} = 1/P_i(X_i)\) für Beobachtungen in der Behandlungsgruppe und \(w_{i,K} = 1/(1-P_i(X_i))\) für Beobachtungen in der Kontrollgruppe, \[\begin{align} \tau^{\text{IPW}} = \frac{1}{n}\sum_{i=1}^n \left[\frac{B_i Y_i}{P_i(X_i)} - \frac{(1-B_i)Y_i}{1-P_i(X_i)} \right].\label{eq:tauipw} \end{align}\]
Grundsätzlich ist TreatmentPropensity eine nicht beobachtbare Variable und muss daher aus den Daten geschätzt werden. Eine geschätzte Behandlungswahrscheinlichkeiten \(\widehat{P}_i(X_i)\) wird als Propensity Score bezeichnet. In der Praxis erfolgt die Schätzung von Propensity Scores meist mit logistischer Regression. Ein erwartungstreuer Schätzer des ATE ist \[\begin{align} \widehat{\tau}^{\text{IPW}} = \frac{1}{n}\sum_{i=1}^n \left[\frac{B_i Y_i}{\widehat{P}_i(X_i)} - \frac{(1-B_i)Y_i}{1-\widehat{P}_i(X_i)} \right].\label{eq:hattauipw} \end{align}\] Hirano, Imbens, und Ridder (2003) diskutieren Alternativen zu \(\eqref{eq:hattauipw}\) für die Schätzung anderer Typen von Behandlungseffekten.
Wir schätzen nachfolgend die Propensity Scores für unser Anwendungsbeispiel, erläutern die Berechnung der Gewichte sowie die Schätzung von Behandlungseffekten mit gewichteter Regression. Hierbei betrachten wir eine Variante von \(\eqref{eq:hattauipw}\) mit normalisierten Gewichten \(\tilde{w}_{i,B} = w_{i,B}/\sum_i w_{i,B}\) und \(\tilde{w}_{i,K} = w_{i,K}/\sum_i w_{i,K}\) die sich jeweils zu 1 summieren.3 Dies ergibt den Hájek-Schätzer4 \[\begin{align} \widehat{\tau}_N^{\text{IPW}} = \frac{\sum_i\tilde{w}_{i,B}Y_i}{\sum_i\tilde{w}_{i,B}} - \frac{\sum_i\tilde{w}_{i,K}Y_i}{\sum_i\tilde{w}_{i,K}}.\label{eq:hattauhajek} \end{align}\]
3 Eine Normalisierung der Gewichte reduziert die Varianz des Schätzers, vgl. Hirano, Imbens, und Ridder (2003)
4 Siehe Hájek (1971).
Zunächst Schätzen wir ein logistisches Regressionsmodell mit age, male und hours als erklärende Variablen für dark_mode.
Die Propensity Scores erhalten wir als angepasste Werte aus der Regression darkmode_ps_logit mit fitted(). Wir erweitern den Datensatz mit den Ergebnissen.
Zur Beurteilung der Überlappung (vgl. Annahme \(\eqref{eq:overlap}\) können wir die Verteilung der Propensity Scores nach Behandlungs-Indikator mit Histogrammen visualisieren.
Ein Vergleich der Histogramme zeigt, dass die Überlappung der Propensity Scores in der linken Flanken der Verteilungen der Kontrollgruppe und in der rechten Flanke der Behandlungsgruppe schlechter wird. Wir entfernen zunächst Beobachtungen aus der Stichprobe deren Propensity Scores wenig bzw. keine Überlappung aufweisen.
IPWs anhand der Propensity Scores können leicht mit der Vorschrift \[\begin{align} \text{IPW} = \frac{\text{dark\_mode}}{\text{PS}} + \frac{1 - \text{dark\_mode}}{1 - \text{PS}}, \end{align}\] berechnet werden.
Eine Schätzung des durchschnittlichen Behandlungseffekts gemäß \(\eqref{eq:hattauhajek}\) implementieren wir mit dplyr.
Diese Schätzung des Behandlungseffekts ist äquivalent zur gewichteten KQ-Schätzung anhand eines einfachen linearen Regressionsmodells.
Unsere Schätzung des ATE ist der geschätzte Koeffizient von dark_mode. Die ausgegebenen Standardfehler und Inferenzstatistiken sind jedoch ungültig aufgrund der Gewichtung mit IPWs, den inversen geschätzten Wahrscheinlichkeiten für eine Behandlung. Der Grund hierfür ist, dass die Berechnung der Standardfehler in summary() die zusätzliche Unsicherheit durch die geschätzen Propensity Scores nicht berücksichtigt! Später im Kapitel erläutern wir die Berechnung gültiger Standardfehler für IPW-Schätzer basierend auf Propensity Scores mit dem Bootstrap.
6.2 Selektierende Matching-Verfahren
Das grundsätzliche Konzept von selektierendem Matching wird in der nachstehenden interaktiven Grafik veranschaulicht. Hier betrachten wir beobachtete Ausprägungen von zwei (unabhängig und identisch verteilten) Matching-Variablen für Subjekte in der Behandlungsgruppe (blau) sowie Kontrollgruppe (rot). Als Matches qualifizieren sich sämtliche Beobachtungen der anderen Gruppe, deren Euklidische Distanz zu dem ausgewählten Punkt das über den Slider eingestellte maximale Caliper nicht überschreitet.5 Diese Region wird durch den gestrichelten Kreis gekennzeichnet und kann über den zugehörigen Slider angepasst werden. Per Klick auf eine Beobachtung werden Matches aus der anderen Gruppe durch eine verbindende Linie und farbliches Hervorheben kenntlich gemacht. Mit dem Slider für \(k\) wird festgelegt, dass lediglich die nahesten \(k\) qualifizierten Beobachtungen als Matches behandelt werden. Die Grafik illustriert inbesondere, dass Beobachtungen in Abhängigkeit von \(k\) und dem Caliper (falls gewünscht) mehrfach (s.g. Matching mit Zurücklegen) oder auch gar nicht gematcht werden können.
5 Es handelt sich hierbei um einen Spezialfall von Matching anhand der Mahalanobis-Distanz, siehe Kapitel 6.2.3.
Linking to librsvg 2.56.3
Für die nachfolgenden Code-Beispiele verwenden wir das R-Paket MatchIt. MatchIt::matchit() nutzt standardmäßig Eins-zu-Eins-Matching (ohne Zurücklegen) von Beobachtungen der Treatment-Gruppe mit Beobachtungen der Kontrollgruppe.6 Die für das Matching zu verwendenden Variablen werden über das Argument formula als Funktion des Behandlungsindikators definiert. matchit() bereitet das Objekt für eine Schätzung des ATT mit einer geeigneten Funktionen, s. ?matchit und hier insb. die Erläuterungen der Argumente replace = F, ratio = 1 und estimand = "ATT" für Details. Mit cobalt::balt.tab() erhalten wir eine balance table für den gematchten Datensatz.
6 Dieses Schema zielt auf eine Schätzung des ATT ab.
Wir zeigen als nächstes, wie MatchIt::matchit() für Matching anhand der Regressoren age, hours, und male in unserem Website-Beispiel für unterschiedliche Varianten durchgeführt werden kann.
6.2.1 Exaktes Matching
Exaktes Matching ordnet einem Subjekt aus der Behandlungsgruppe ein oder mehrere Subjekte aus der Kontrollgruppe zu, wenn die boebachteten Ausprägung der Matching-Variablen exakt übereinstimmen. Hierbei muss die ‘Distanz’ zwischen den Ausprägung der Matching-Variablen folglich \(0\) sein. Dieses Verfahren findet meist bei ausschließlich diskret verteilten Merkmalen Anwendung. Bei kontinuierlich verteilten Merkmalen (vgl. die obige interaktive Grafik) sind exakte Matches zwar theoretisch unmöglich, ergeben sich jedoch in der Praxis aus der Datenerfassung, bspw. durch Rundungsfehler. In matchit() erhalten wir exaktes Ein-zu-eins-Matching mit method = "exact".
Aufgrund der kontinulierlich verteilten Variable hours gibt es in unserem Website-Beispiel keine exakten Matches. Dieses Verfahren ist hier folglich ungeeignet.
6.2.2 Coarsened Exact Matching
Bei dieser Methode werden kontinuierliche Matching-Variablen grob (Engl. coarse) klassiert, ähnlich wie bei einem Histogram. Diese Diskretisierung ermöglicht es exakte Übereinstimmungen zwischen Behandlungs- und Kontrollgruppenbeobachtungen hinsichtlich ihrer klassierten Ausprägungen zu finden. Sowohl Behandlungs- als auch Kontrollbeobachtungen die mindestents einen exakten Match haben, werden Teil des gematchten Datensatzes. In matchit() wird Coarsened Exact Matching mit method = "cem" durchgeführt. Über das Argument cutpoints geben wir an, dass hours in 6 Klassen und age in 4 Klassen eingeteilt werden soll.7 Mit k1k = TRUE erfolgt Eins-zu-eins-Matching: Bei mehreren exakten Matches wird die Beobachtung mit der geringsten Mahalanobis-Distanz (für die unklassierten Matching-Variablen) gewählt.
7 Diese Werte wurden ad-hoc gewählt da sie zu einem guten Ergebnis führen.
Mit Coarsened Exact Matching erhalten wir einen Datensatz mit 82 Beobachtungen und guter Balance.
6.2.3 Matching mit der Mahalanobis-Distanz
Die Euklidische Distanz misst den direkten Abstand zwischen zwei Punkten und ist nicht invariant gegenüber Transformationen, insbesondere bei unterschiedlichen Skalierungen und bei Korrelation der Matching-Variablen. Die Mahalanobis-Distanz hingegen ist ein standardisiertes Distanzmaß, das unter Berücksichtigung der Varianz-Kovarianz-Struktur der Daten angibt, wie viele Standardabweichungen zwei Datenpunkte voneinander entfernt sind. Die Mahalanobis-Distanz ist invariant gegenüber linearen Transformationen (Skalierung, Translation und Rotation) der Daten und bietet ein genaueres Maß für die Unähnlichkeit zweier Beobachtungen hinsichtlich ihrer Ausprägungen der Matching-Variablen.
Betrachte die Datenpunkte \(P_1=(X_1,Y_1)'\) und \(P_2=(X_2,Y_2)'\) für die Matching-Variablen \(X\) und \(Y\). Die Mahalanobis-Distanz zwischen \(P_1\) und \(P_2\) ist definiert als \[\begin{align*} d_M(P_1,\,P_2) = \sqrt{(P_1 - P_2)'\boldsymbol{S}^{-1} (P_1 - P_2)}, \end{align*}\] wobei \(\boldsymbol{S}\) die Varianz-Kovarianz-Matrix von \(X\) und \(Y\) ist. Die Mahalanobis-Distanz \(d_M(P_1,\,P_2)\) ist also die Euklidische Distanz zwischen den standardisierten Datenpunkten.
In empirischen Anwendungen ersetzen wir die (unbekannten) Komponenten der Varianz-Kovarianz-Matrix durch Stichprobenvarianten. Dies ergibt die Formel
\[\begin{align*} \widehat{d}_M(P_1,\,P_2) = \sqrt{ \begin{pmatrix} X_1 - X_2\\ Y_1 - Y_2 \end{pmatrix}' \begin{pmatrix} \widehat{\text{Var}}(X^2) & \widehat{\text{Cov}}(X, Y) \\ \widehat{\text{Cov}}(X, Y) & \widehat{\text{Var}}(Y^2) \end{pmatrix}^{-1} \begin{pmatrix} X_1 - X_2\\ Y_1 - Y_2 \end{pmatrix} }. \end{align*}\]
Die nachstehende interaktive Grafik zeigt Beobachtungen zweier Matching-Variablen, die aus einer bivariaten Normalverteilung mit positiver Korrelation generiert werden. Diese bivariate Verteilung ist identisch für Beobachtungen aus der Kontrollgruppe (rot) und Beobachtungen aus der Behandlungsgruppe (blau).8 Für eine ausgewählte Beobachtung aus der Behandlungsgruppe (schwarzer Rand) werden potentielle Matches in der Kontrollgruppe innerhalb der vorgegebenen Mahalanobis-Distanz in Cyan kenntlich gemacht (und umgekehrt). Beachte, dass die Mahalanobis-Distanz die Varianzen und Kovarianzen der Daten berücksichtigt, sodass die gematchten Beobachtungen in einem elliptischen Bereich um die betrachtete Beobachtung liegen. Die Form dieser Ellipse verändert sich mit der Korrelation zwischen den Matching-Variablen, die über den Slider Correlation festgelegt ist. Matching anhand der euklidische Distanz ignoriert die Skalierung und und die Abhängigkeit der Matching-Variablen: Die gematchten Beobachtungen liegen innerhalb eines entsprechenden Kreises (gestrichelte Linie) um die interessierende Beobachtung. Je nach Kovarianz-Struktur können die Mengen der gematchten Beobachtungen sich deutlich unterscheiden. Beachte, dass die euklidische Distanz hier der Spezialfall der Mahalanobis-Distanz bei unkorrelierten Daten ist.
8 Insbesondere ist die Varianz beider Matching-Variablen identisch.
Für Eins-zu-Eins-Matching im Website-Beispiel anhand der Mahalanobis-Distanz mit matchit() setzen wir distance = "mahalanobis" und wählen method = "nearest". Mit diesen Parametern wird jeder Beobachtung aus der Behandlungsgruppe die gemäß \(d_M\) am ehesten vergleichbarste Beobachtung aus der Kontrollgruppe zugewiesen, wobei keine mehrfachen Matches zulässig sind.
Die Ergebnisse zeigen, dass für sämtliche \(149\) Beobachtungen aus der Behandlungsgruppe ein individueller Match in der Kontrollgruppe gefunden werden konnte. Es werden lediglich \(2\) Beobachtungen der \(151\) Beobachtungen in der Kontrollgruppe nicht gematcht.
Entsprechend zeigt die Balance-Table eine ähnliche Diskrepanz beider Gruppen hinsichtlich der Matching-Variablen an.
Mahalanobis-Distanz mit Caliper .25 für Propensity Scores basierend auf logistischer Regression
Für eine strengeres Matching-Kriterium kann ein Caliper, d.h. eine maximal zulässige Distanz, herangezogen werden. Die Mahalanobis-Distanz hat jedoch keine einheitliche Skala: Ob eine Distanz als groß oder klein betrachten werden kann, hängt von der Anzahl der Matching-Variablen und dem Überlappungsgrad zwischen den Gruppen ab. Daher wird die Beschränkung durch einen Caliper nicht auf \(\widehat{d}_M\) sondern auf Propensity Scores angewendet.
Im nächsten Code-Beispiel spezifizieren wir mit distance = "glm", dass Propensity Scores gemäß der Vorschrift in formula geschätzt werden. Mit mahvars = ~ age + male + hours legen wir die Matching-Variablen für die Berechnung von \(\widehat{d}_M\) fest. caliper = .25 legt fest, dass lediglich Beobachtungen der Kontrollgruppe bei einer absoluten Differenz der Propensity Scores von höchstens \(0.25\) Standardabweichungen als Match für eine Beobachtung in der Behandlungsgruppe qualifiziert sind.
Die Balance-Table zeigt einen deutlichen Effekt der Beschränkung qualifizierter Beobachtungen durch caliper = .25: Aufgrund der oberen Grenze für die Propensity-Score-Differenz von XX wird für lediglich \(104\) Beobachtungen aus der Behandlungsgruppe ein individueller Match in der Kontrollgruppe gefunden.9 Weiterhin finden wir eine verbesserte Balance für den gematchten Datensatz.
9 Die durch caliper implizierte Obergrenze ergibt sich als .25 * sd(fitted(darkmode_ps_logit))).
6.2.4 Propensity Score Matching
Eine gängige Variante ist Matching ausschließlich anhand von Propensity Scores innerhalb eines Calipers.
Laut Balance-Table führt Eins-zu-Eins-Matching basierend auf Propensity Scores zu einem Datensatz mit \(104\) gematchten Beobachtungen in der Behandlungsgruppe. Hinsichtlich der standardisierten Mittelwertdifferenz (Diff.Adj) erzielt diese Methode die beste Balance unter den betrachteten Ansätzen.
Vergleich der Balance verschiedener Verfahren mit Love-Plot
Standardisierte Mittelwertdifferenzen für verschiedene Matching-Verfahren können grafisch mit einem Love-Plot (Love 2004) veranschaulicht werden. Hierzu nutzen wir cobalt::love.plot() und übergeben die mit matchit() generierten Objekte im Argument weights.
Die Grafik zeigt, dass Coarsened Exact Matching (CEM) unter allen betrachteten Verfahren die Stichprobe mit der besten Balance ergibt. Diesen gematchten Datensatz erhalten wir mit MatchIt::match.data().
darkmode_matched enthält Gewichte (weights) für die jeweilige Gruppe zu denen gemachte Beobachtungen gehören (subclass). Dies ist relevant, falls Beobachtungen mehrfach gematcht werden. Wegen Eins-zu-eins-Matching ohne Zurücklegen gibt es in unserem Beispiel XX Beobachtungspaare und sämtliche Gewichte sind 1. Die Berücksichtigung der Gewicht in den nachfolgenden Aufrufen von Schätzfunktionen (bspw.lm()) ist daher nicht nötig und erfolgt lediglich zur Illustration der grundsätzlichen Vorgehensweise.
Eine Wiederholung der grafischen Analyse in Kapitel 6.1 zeigt eine deutlich verbesserte Vergleichbarkeit hinsichtlich der Verteilung der Matching-Variablen in darkmode_matched.
Wir beobachten eine bessere Balance bei age und hours. Inbesondere ist male für Kontroll- und Behandlungsgruppe ausgeglichen.
6.3 Schätzung und Inferenz für den Behandlungseffekt nach Matching
Wir schätzen nun den Behandlungseffekt von dark_mode auf read_time für die mit CEM und Propensity Score Matching ermittelten Datensätzen und vergleichen die Ergebniss anschließend mit einer Regressionsschätzung ohne Matching.
Wir kombinieren die Matching-Verfahren mit linearer Regression, d.h. wir Schätzen den Behandlungseffekt anhand es gematchten Datensatzes als Mittelwertdifferenz nach zusätzlicher Kontrolle für die Matching-Variablen. Diese Kombination von Matching mit Regression wird in der Literatur als Regression Adjustment bezeichnet und ist insbesondere hilfreich, wenn Backdoors mit Matching geschlossen werden sollen, der kausale Effekt jedoch nur unter Verwendung einer nicht-trivialen Regressionsfunktion ermittelt werden kann. Zum Beispiel kann bei einer kontinuierlichen Behandlungsvariable und einem nicht-linearen Zusammenhang mit \(Y\) der kausale Effekt nicht durch einen bloßen Mittelwertvergleich erfasst werden, sondern erfordert eine adäquate Modellierung dieses Zusammenhangs in der Regressionsfunktion. Die zusätzliche Kontrolle für Matching-Variablen kann die Varianz der Schätzung verringern und das Risiko einer verzerrten Schätzung abmildern, falls nach Matching noch Unterschiede in der Balance von Behandlungs- und Kontrollgruppe vorliegen.
Beachte, dass für die gematchten Datensätze jeweils ein durchschnittlicher Behandlungseffekt für die Beobachtungen mit erfolgter Behandlung ermittelt wird: In sämtlichen oben gezeigten Matchibg-Verfahren werden mit estimand = "ATT" vergleichbarere Kontrollbeobachtungen für die behandelten Beobachtungen ermittelt. Wir schätzen den Effekt der Behandlung, indem wir die Ergebnisse von behandelten Personen mit denen von gematchten Personen vergleichen, die keine Behandlung erhalten haben. Diese Vergleichsgruppe dient als Ersatz für den hypothetischen Zustand der Behandlungsgruppe, wenn keine Behandlung erfolgt wäre. Dies entspricht der Definition eines ATT — ein average treatment effect on the treated.
6.3.1 Cluster-robuste Standardfehler
Für Matching-Verfahren sind die mit summary() berechneten Standardfehler (und damit auch Konfidenzintervalle, t-Statistiken und p-Werte) für den Behandlungseffekt grundsätzlich ungültig. Je nach Matching-Verfahren liegen unterschiedliche Quellen von Schätzunsicherheit vor, die bei der Berechnung von Standardfehlern zusätzlich zu der “üblichen” Stichproben-Variabilität berücksichtig werden müssen. Gründe hierfür sind der Matching-Prozess ansich oder weitere Unsicherheit durch die Schätzung zusätzlicher Parameter, etwa bei der Berechnung von Propensity Scores mit logistischer Regression. Eine weiterere Ursache zusätzlicher Variation durch den Matching-Prozess, die wir bisher nicht näher betrachtet haben ensteht durch Zurücklegen, d.h. wenn Beobachtungen mehrfach gematcht werden können. Auch dieser Faktor wird in der von summary() verwendeten Formel für den Standardfehler des Effekt-Schätzers nicht berücksichtigt.
Die Standardfehlerberechnung für Matching-Schätzer von Behandlungseffekten ist ein Gegenstand aktueller methodischer Forschung. P. C. Austin und Small (2014) und Abadie und Spiess (2022) belegen die Gültigkeit von cluster-robusten Standardfehlern mit Clustering auf Ebene der Beobachtungsgruppen (subclass im output von match.data()) bei Matching ohne Zurücklegen. Für Matching anhand von Propensity Scores (auch mit Zurücklegen) zeigt Imbens (2016), dass ignorieren der zusätzlichen Unsicherheit durch die Schätzung der Propensity Scores zu konservativer Inferenz für den ATE anhand eines cluster-robusten Standardfehlerschätzers führt, jedoch ungültige Inferenz für die Schätzung des ATT bedeuten kann. Ähnlich zu P. C. Austin und Small (2014) deuten die Ergebnisse der Simulationsstudie von Bodory u. a. (2020) jedoch auf grundsätzlich bessere Eigenschaften der Schätzung hin, wenn die Standardfehler nicht für die Schätzung der Propensity Scores adjustiert werden.
Weiterhin ist die Kontrolle für Kovariablen mit Erklärungskraft für die Outcome-Variable und für die Matching-Variablen mit Regression Adjustment (vgl. Kapitel 6.3) für die Schätzung des Behandlungseffekts nach Matching eine etablierte Strategie, vgl. Hill und Reiter (2006) und Abadie und Spiess (2022). So können die Varianz der Schätzung und das Risiko einer Verzerrung der Standardfehler aufgrund verbleibender Imbalance von Behandlungs- und Kontrollgruppe nach Matching verringert werden.
Zur Demonstration von (cluster)-robuster Inferenz und für eine tabellarische Zusammenfassung der Ergebnisse nutzen wir die Pakete marginaleffects und modelsummary. Mit marginaleffects::avg_comparisons() können p-Werte und Kofindenzintervalle unter Berücksichtigung von robuster Standardfehlern und der Gewichte aus dem Matching-Verfahren berechnet werden.
Wir fassen die Ergebnisse mit modelsummary::modelsummary() tabellarisch zusammen.
6.4 Bootstrap-Schätzung kausaler Effekte bei Matching
Ein Bootstrap-Verfahren generiert mit Resampling (wiederholtes Ziehen mit Zurücklegen) aus dem Original-Datensatz (viele) künstliche Datensätze, für die der Schätzer (d.h. das gesamte Verfahren inkl. Matching!) jeweils berechnet wird. Die Verteilung der so gewonnenen Bootstrap-Schätzwerte approximiert die wahre, unbekannte Stichprobenverteilung des Schätzers des Behandlungseffekts. Mit dieser simulierten Verteilung können wir Inferenz betreiben: Wir können einen Bootstrap-Punktschätzer des Behandlungseffekts (Stichprobenmittel der Bootstrap-Schätzungen) sowie Standardfehler (Standardabweichung der der Bootstrap-Schätzungen) und p-Werte berechnen.
Der Bootstrap kann hilfreich sein, wenn unklar ist, wie Standardfehler für die Unsicherheit des Matching-Prozesses zu adjustieren sind, um gültige Inferenzsstatistiken zu erhalten. Abadie und Imbens (2008) zeigen analytisch, dass der Standard-Bootstrap die Stichprobenverteilung für Schätzer kausaler Effekte anhand von gematchten Datensätzen (d.h. bei Zuordnung/Selektion von Beobachtungen mit Matching) nicht korrekt abbilden kann. Grundsätzlich problematisch hierbei ist, wenn der Bootstrap eine verzerrte Schätzung produziert und/oder zu kleine Standardfehler liefert. Abadie und Imbens (2008) belegen die Tendenz des Bootstraps zu konservative (d.h. zu große) Standardfehler zu produzieren. Simulationsnachweise (Bodory u. a. 2020; Hill und Reiter 2006; P. C. Austin und Small 2014; P. C. Austin und Stuart 2017) finden, dass Bootstrap-Standardfehler u.a. bei Propensity Score Matching mit Zurücklegen leicht konservativ sind somit das gewünschte nominale Signifikanzniveau eines Bootstrap-Hypothesentests nicht überschritten wird, weshalb der Standard-Bootstrap trotz seiner Schwächen in der empirischen Forschung oft angewendet wird.
Wir betrachen als nächstes einen Bootstrap-Algorithmus für Inferenz bezüglich eines kausalen Effekts nach Matching und demonstrieren die Schätzung anhand unseres Website-Beispiels für den ATT nach Propensity-Score-Matching.
Wir Implementieren nun einen Bootstrap-Schätzer des ATT im Website-Beispiel für Propensity-Score-Matching. Hierzu definieren wir eine R-Funktion boot_fun() für die Schritte 1 und 2 im obigen Algorithmus.
Wir berechnen nun eine Bootstrap-Schätzung des ATT von dark_mode auf readingtime nach Propensity Score Matching mit einem caliper von 0.25 sowie den zugehörigen Standardfehler und ein 95%-KI mit der zuvor definierten Funktion boot_fun.
Den Bootstrap-Schätzer des ATT sowie den Bootstrap-Standardfehler berechnen wir mit mean() und sd() anhand der \(B=199\) Bootstrap-Replikationen in boot_out$t.
Beachte, dass der Bootstrap-Schätzer des Behandlungseffekts nicht unmittelbar von boot() ausgegeben wird: Der Eintrag original ist die Schätzung anhand der ursprünglichen gesamten Stichprobe und bias ist die Differenz zwischen dieser Schätzung und dem Mittelwert der Bootstrap-Schätzungen.
Wir können prüfen, dass die Berechnung des Standardfehlers dem Stichprobenstandardabweichung der Bootstrap-Schätzungen entspricht.
Ein 95%-Konfidenzintervall für den kausalen Effekt erhalten wir mit boot::boot.ci().10
10 type = "bca" (bias-corrected accelerated) ist eine gängige Implementierung für die Berechnung des Bootstrap-Konfidenz-Intervalls.
Beachte, dass der Bootstrap-Standardfehler sowie das Bootstrap-Konfidenzintervall nahe der mit avg_comarisons berechneten Werte für sum_PSC sind.
Bootstrap-Standardfehler für IPW-Schätzer des ATE berechnen
Die Bootstrap-Funktion boot_fun kann leicht für eine Schätzung des Standardfehlers für den IPW-Schätzer des ATE aus Kapitel 6.1.3 angepasst werden. Statt einer Matching-Prozedur berechnen wir hierzu für \(B\) Bootstrap-Stichproben den Schätzer \(\widehat{\tau}^\text{IPW}\) mit Trimming der Propensity Scores.
In diesem Fall ist der Bootstrap-Standardfehler gut mit dem anhand von summary(model_ipw) berechneten Standardfehler vergleichbar. Ein 95%-Konfidenzintervall für den ATE erhalten wir wie zuvor mit boot.ci().
6.5 Regression, Matching und Doubly-Robust-Schätzung
Als Doubly-Robust-Schätzer bezeichnet man Methoden, die trotz Fehlspezifikationen im Matching-Verfahren oder bei falscher Form der Regressionsfunktion für die Outcome-Variable eine robuste Schätzung des kausalen Effekts ermöglichen.11 Unter der Voraussetzung, dass die verwendeten Matching-Variablen sämtliche Backdoors schließen, ist mit Doubly-Robust-Schätzung eine konsistente Schätzung des Behandlungseffekts unter abgeschwächten Annahmen gewährleistet. Dies macht Doubly-Robust-Schätzer besonders nützlich in empirischen Anwendungen, wenn Modellspezifikationen herausfordernd sind.
11 Bspw. kann eine falsche funktionale Form bei logistischer Regression eine verzerrte Schätzung von Propensity Scores und damit eine unzureichende Balance bedeuten.
Der von Wooldridge (2010) vorgeschlagene Doubly-Robust-Schätzer (Inverse Probability Weighted Regression Adjustment, IPWRA) für den ATE erreicht seine vorteilhaften Eigenschaften durch eine geschickte Kombination von IPW und Regression Adjustment.
Wir implementieren den Doubly-Robust-Schätzer des ATE von Wooldridge (2010) für das Website-Beispiel in der Funktion IPWRA() under Adaption des Schemas von IPW_boot().
Schätzung des ATE mit IPWRA():
Wie bei \(\widehat{\tau}^\text{IPW}\) gibt es keine formale Darstellung des Standardfehlers für den IPWRA-Schätzer, sodass auch hier der Bootstrap für Inferenz genutzt werden muss.
Wir berechnen die interessierenden Statistiken analog zu der Vorgehensweise in Kapitel 6.4.
Die Bootstrap-Schätzung des ATE mit IPWRA ist mit den Ergebnissen für die Bootstrap-Variante von \(\widehat{\tau}^\text{IPW}\) vergleichbar, hat allerdings einen etwas kleineren Standardfehler.
6.6 Zusammenfassung
In diesem Kapitel haben wir gezeigt, wie Matching-Methoden genutzt werden, um kausale Effekte aus nicht-experimentellen Daten zu schätzen. Selektierende Ansätze wie exaktes Matching, Coarsened Exact Matching und Propensity Score Matching stellen die Vergleichbarkeit von Behandlungs- und Kontrollgruppen her, indem gezielt Beobachtungen mit ähnlich verteilten Kovariablen ausgewählt werden. Gewichtendes Matching passt Beobachtungen durch Gewichte an, um die Verteilung von Counfoundern für beide Gruppen statistisch zu anzugleichen, ohne eine explizite Selektion von Beobachtungen vorzunehmen. Fortgeschrittene Verfahren kombinieren Matching mit Regression Adjustment, um verbleibende Unterschiede mit parametrischer Modellierung auszugleichen, während Doubly Robust Estimation selbst bei Fehlspezifikationen im Outcome-Modell oder des Matchings konsistente Schätzungen liefert. Bootstrap-Verfahren ermöglichen die Berechnung verlässlicher Standardfehler und verbessern die Robustheit der Inferenz, gerade für komplexe Matching-Ansätze.
Wir haben gezeigt, wie diese Methoden effizient in R implementiert werden können, insbesondere mit Paketen wie MatchIt und weightit für das Matching sowie boot für Bootstrap-basierte Ansätze.