PS-Gewichtung
1 Der Datensatz
Wir möchten evaluieren, ob ein entwickeltes Englischtraining die
Englischkenntnisse von Studierenden erhöht. Dazu wurde eine
quasi-experimentelle Studie durchgeführt, bei der die Studierenden
selbst entscheiden konnten, ob sie am Training teilnehmen
(Interventionsgruppe; IG) oder nicht (Kontrollgruppe; KG). Nach dem
Training wurden in beiden Gruppen die Englischkenntnisse gemessen. Als
Kovariaten wurden die Englischvorkenntnisse und die Vorliebe für
Englisch erhoben. Für die Analysen nutzen wir den simulierten Datensatz
(englisch.rda) mit dem Objektnamen englisch
mit den folgenden Variablen:
- treat: Teilnahme am Englischtraining (0 = kein Englischtraining, 1 = Englischtraining) (Interventionsvariable)
- eng_post: Englischkenntnisse nach dem Training (Skala von 0 bis 5) (Outcome-Variable)
- eng_pre: Englischvorkenntnisse (Skala von 0 bis 5) (Kontinuierliche Kovariate)
- eng_vorliebe: Vorliebe für Englisch (0 = keine Vorliebe für Englisch, 1 = Vorliebe für Englisch) (Dichotome Kovariate)
2 Erklärung der R-Funktionen
Für die PS-Gewichtung verwenden wir das Paket WeightIt.
Mit der im Paket enthaltenen Funktion weightit() werden die
Gewichte zur Balancierung der Kovariaten bestimmt. Die folgende Tabelle
bietet eine Übersicht über die wichtigsten Argumente der Funktion.
| Argument | Bedeutung |
|---|---|
weightit( |
|
treat ~ x1 + x2, |
Angabe einer Interventionsvariable (treat)
und einer oder mehrerer Kovariaten (x1 + x2 + ...) |
data, |
Datensatz |
method, |
Wahl der Methode, die dazu genutzt wird, die Gewichte
zu schätzen; z.B. method = "ps" für PS-Gewichtung
(default) |
estimand) |
Angabe, welcher kausale Effekt geschätzt werden soll
(estimand = "ATE" für den durchschnittlichen kausalen
Effekt (default), estimand = "ATT" für den
durchschnittlichen kausalen Effekt der Behandelten,
estimand = "ATC" den durchschnittlichen kausalen Effekt der
Unbehandelten) |
Mit der Funktion bal.tab() aus dem Paket
cobalt können Kennwerte zur Prüfung der Balance gewichteter
Daten ausgegeben werden. Die relevanten Argumente der Funktion werden in
der folgenden Tabelle erklärt.
| Argument | Bedeutung |
|---|---|
bal.tab( |
|
x, |
weightit-Objekt |
stats, |
Kennwerte, die dargestellt werden sollen; z.B.
stats = "mean.diffs" oder stats = "m" für
Mittelwertdifferenzen, stats = "variance.ratios" oder
stats = "v" für Varianzverhältnisse |
continuous |
Angabe, ob bei kontinuierlichen Variablen
standardisierte Mittelwertsdifferenzen (continuous = "std";
default) oder unstandardisierte Mittelwertsdifferenzen
(continuous = "raw") ausgegeben werden sollen |
binary |
Angabe, ob bei dichotomen Variablen standardisierte
Mittelwertsdifferenzen (continuous = "std") oder
unstandardisierte Mittelwertsdifferenzen
(continuous = "raw"; default) ausgegeben werden sollen |
Mit der Funktion love.plot() aus dem Paket
cobalt können Unterschiede in Verteilungskennwerten auf den
Kovariaten vor und nach der Gewichtung visualisiert werden. Die
relevanten Argumente der Funktion sind in der folgenden Tabelle
erklärt.
| Argument | Bedeutung |
|---|---|
love.plot( |
|
x, |
weightit-Objekt |
stats, |
Kennwerte, die dargestellt werden sollen; z.B.
stats = "mean.diffs" für Mittelwertdifferenzen,
stats = "variance.ratios" für Varianzverhältnisse |
thresholds, |
Markierung von Grenzwerten für eine akzeptable Balance
(optional); z.B. nach Steiner et al. (2010) für Mittelwertsdifferenzen
thresholds = c(-0.1, 0.1) |
sample.names, |
frei wählbare Namen für die Beschriftungen der
ungewichteten und gewichteten Stichprobe in der Grafik; z.B.
sample.names = c("ungewichtet", "gewichtet") |
stars) |
Angabe von Variablenarten, die in der Grafik mit einem
Stern markiert werden sollen; z.B. stars = "raw" für
Variablen mit nicht-standardisierten Mittelwertdifferenzen |
Mit der Funktion bal.plot() aus dem Paket
cobalt können die Verteilungen der Kovariaten zwischen den
Interventionsgruppen vor und nach der PS-Gewichtung visualisiert werden.
Die relevanten Argumente der Funktion sind in der folgenden Tabelle
erklärt.
| Argument | Bedeutung |
|---|---|
bal.plot( |
|
x, |
matchit-Objekt |
var.name, |
Variable, für die die Kovariatenbalance dargestellt werden soll |
which, |
Angabe, welche Verteilungen dargestellt werden sollen;
z.B. which = "both" für beide Stichproben;
which = "unadjusted" nur für Gesamtstichprobe;
which = "adjusted" nur für gewichtete Stichprobe |
type, |
Art der Darstellungsform; z.B.
type = "histogram" für ein Histogramm;
type = "density" für Dichte |
mirror, |
Spiegelung gruppenspezifischer Verteilung an der
x-Achse bei mirror = TRUE |
sample.names) |
frei wählbare Namen für die ungewichtete und gewichtete
Stichprobe in der Grafik; z.B.
sample.names = c("Gesamtstichprobe", "Gewichtete Stichprobe") |
3 Vorbereitung
Zunächst installieren (falls nötig) und laden wir die nötigen Pakete:
# install.packages("WeightIt") # Installation vor der ersten Nutzung des Pakets WeightIt
library(WeightIt) # Laden des Pakets WeightIt
# install.packages("cobalt") # Installation vor der ersten Nutzung des Pakets cobalt
library(cobalt) # Laden des Pakets cobaltund lesen den Datensatz ein:
load("englisch.rda") # Einlesen des Datensatzes (wenn der Datensatz nicht im
# Arbeitsverzeichnis liegt, muss der Pfad entsprechend dem
# Speicherort des Datensatzes angepasst werden)4 Overlap herstellen
Um Overlap herzustellen, schätzen wir zunächst die Propensity-Scores für jede Person mit Hilfe einer logistischen Regression der Interventionsvariable auf die beiden Kovariaten und speichern die geschätzten Propensity-Scores im Datensatz.
# logistische Regression zum Schätzen der PS:
ps <- glm(treat ~ eng_pre + eng_vorliebe, data = englisch, family = "binomial")
# Speichern der PS als Variable im bestehenden Datensatz:
ps1 <- ps$fitted
englisch$ps <- ps1Um Personen auszuschließen, die außerhalb der Region of Common Support liegen und somit Overlap herzustellen, erstellen wir einen Datensatz, in dem nur noch Personen enthalten sind, deren Propensity-Scores innerhalb des Wertebereichs liegen, für den es sowohl in der Interventions- als auch in der Kontrollgruppe Werte gibt. Dazu identifizieren wir den kleinsten sowie den größten PS, bei dem es in beiden Interventionsgruppen noch Personen gibt und schließen alle Personen außerhalb dieser Region aus.
# Extraktion des größeren der beiden kleinsten Propensity-Scores in der Interventions- und der Kontrollgruppe sowie
lower_border <- max(tapply(ps1, englisch$treat, min))
# Extraktion des kleineren der beiden größten Propensity-Scores in der Interventions- und der Kontrollgruppe:
upper_border <- min(tapply(ps1, englisch$treat, max))
# Ausschluss der Personen außerhalb der Region of Common Support bzw.
# Einbezug nur von Personen innerhalb der Region of Common Support:
englisch.overlap <- englisch[ps1 >= lower_border & ps1 <= upper_border, ]5 Gewichte bestimmen
Nun können Gewichte bestimmt werden. Wir beziehen uns dafür
ausschließlich auf den Datensatz in dem Overlap hergestellt wurde und
wählen hier beispielhaft den durchschnittlichen kausalen Effekt (ATE)
als Zieleffekt. Mit der weightit werden sowohl
Propensity-Scores geschätzt (im Vergleich zum Abschnitt zuvor jedoch nur
für die Daten innerhalb der Region of Common Support) als auch Gewichte
erstellt.
ps_weighted <- weightit(
treat ~ eng_pre + eng_vorliebe, # Propensity-Score Modell
data = englisch.overlap, # verwendeter Datensatz
method = "ps", # Gewichtung auf Basis der Propensity-Scores
estimand = "ATE") # ATE als ZieleffektÜber den summary()-Befehl können wir uns die Ergebnisse
der Gewichtung ausgeben lassen. Die einzelnen Abschnitte des Outputs
werden nun nacheinander gezeigt und erläutert.
summary(ps_weighted)...
## - Weight ranges:
##
## Min Max
## treated 1.10 |---------------------------| 9.76
## control 1.12 |-------------------------| 9.15
...Unter Weight ranges werden
Minima und Maxima der Gewichte pro Gruppe angegeben. In der
Interventionsgruppe reichen die Gewichte von einem Minimum von
1,10 bis zu einem Maximum von 9,76. In der
Kontrollgruppe ist das Minimum 1,12 und das Maximum
9,15.
...
## - Units with the 5 most extreme weights by group:
##
## 243 224 282 269 255
## treated 4.879 5.0023 5.3069 6.5514 9.7567
## 110 208 212 180 206
## control 5.1715 6.5884 7.3825 8.076 9.1454
##
...Unter Units with 5 greatest
weights by group werden die fünf Personen pro Gruppe mit den
höchsten Gewichten angegeben. In der Interventionsgruppe ist
beispielsweise die Person mit dem höchsten Gewicht (von
9,76) die Person in Zeile 255 des Datensatzes.
Die Größe der Gewichte sollte angeschaut werden, um einen Eindruck davon
zu bekommen, in welchem Ausmaß extreme Gewichte vorliegen. Sind Gewichte
sehr extrem, können diese das Ergebnis stark beeinflussen. Um dies zu
umgehen, besteht die Möglichkeit, Gewichte zu trimmen, d.h. sie auf
einen weniger extremen Wert zu setzen (siehe auch weiterführende
Literatur).
...
## - Effective Sample Sizes:
##
## Control Treated
## Unweighted 221 210
## Weighted 172 168
...Effective Sample Sizes
beschreibt die effektive Stichprobengröße der ungewichteten und der
gewichteten Stichprobe. Es gibt 221 Personen in der
Kontrollgruppe und 210 Personen in der Interventionsgruppe.
Nach der Gewichtung ist die effektive Stichprobengröße und damit die
Präzision der Schätzung geringer.
6 Prüfung der Balance
Wir können die Balance sowohl anhand der Verteilungskennwerte als auch anhand der Verteilungen prüfen.
6.1 Kennwerte
Die Kennwerte zur Prüfung der Balance können wir uns mit der Funktion
bal.tab() aus dem Paket cobalt ausgeben
lassen.
bal.tab(ps_weighted, # weightit-Objekt
stats = c("m", "v"), # Kennwerte, die ausgegeben werden sollen
# hier Mittelwertsdifferenzen ("m") und Varianzverhältnisse ("v")
continuous = "std", # standardisierte Mittelwertsdifferenzen für kontinuierliche
# Variablen
binary = "raw") # unstandardisierte Mittelwertsdifferenzen für binäre Variablen## Call
## weightit(formula = treat ~ eng_pre + eng_vorliebe, data = englisch.overlap,
## method = "ps", estimand = "ATE")
##
## Balance Measures
## Type Diff.Adj V.Ratio.Adj
## prop.score Distance -0.025 0.934
## eng_pre Contin. -0.036 0.723
## eng_vorliebe Binary -0.007 .
##
## Effective sample sizes
## Control Treated
## Unadjusted 221 210
## Adjusted 172 168Unter Balance Measures finden wir die relevanten Kennwerte. Diff.Adj sind die standardisierten (für metrische Variablen) und unstandardisierten (für binäre Variablen) Mittelwertsdifferenzen sowie V.Ratio.Adj die Varianzverhältnisse (für metrische Variablen).
Wir sehen, dass sich die Mittelwertsdifferenzen der
Englischvorkenntnisse (Diff = -0,036) und
Propensity-Scores (Diff = -0,025) innerhalb
der von Steiner et al. (2010) vorgeschlagenen Grenzen für die Balance
metrischer Variablen (Betrag der standardisierten Mittelwertsdifferenz
kleiner als 0,1 Standardabweichungen) befinden. Der Unterschied im
Anteil von Personen mit Vorliebe für Englisch zwischen beiden
Interventionsgruppen ist verschwindend gering (Diff =
-0.007, d.h. 0,7%).
Das Varianzverhätnis der Propensity-Scores beträgt
0,934, das der Englischvorkenntnisse 0,723.
Das Varianzverhältnis der Englischvorkenntnisse liegt somit außerhalb
der von Steiner et al. (2010) vorgeschlagenen Grenzen für akzeptable
Balance (Varianzverhältnisse zwischen 0,8 und 1,25). Sie liegt jedoch
noch innerhalb der weniger strengen Grenzen (zwischen 0,5 und 2) von
Stuart and Rubin (2008). Für die Variable Vorliebe für Englisch wird
kein Varianzverhältnis ausgegeben, da es sich hier um eine binäre
Variable handelt.
Die Mittelwertsdifferenzen und Varianzverhältnisse können wir grafisch darstellen:
love.plot(ps_weighted, # weightit-Objekt
stats = "mean.diffs", # Mittelwertsdifferenzen
thresholds = c(-0.1, 0.1), # Einzeichnen der Grenzwerte nach
# Steiner et al. (2010)
sample.names = c("ungewichtet", # Namen für Stichproben
"gewichtet"),
stars = "raw" # Markierung nicht-standardisierter
# Mittelwertsdifferenzen
)
Abbildung 1: (standardisierte) Mittelwertsdifferenzen der Kovariaten und Propensity-Scores vor (rote Punkte) und nach (türkisfarbene Punkte) der PS-Gewichtung. Die durchgezogene senkrechte Linie kennzeichnet eine Mittelwertsdifferenz von 0 und somit gleiche Mittelwerte in der Interventionsgruppe und Kontrollgruppe. Die gestrichelten Linien kennzeichnen die Grenzen für akzeptable Mittelwertsdifferenzen nach Steiner et al. (2010). (* kennzeichnet Variablen für die die unstandardisierten Mittelwertsdifferenzen angegeben wurden.)
Die Grafik gibt für die Propensity-Scores und die Vorkenntnisse die standardisierten Mittelwertsdifferenzen und für die Vorliebe für Englisch (mit Stern gekennzeichnet) die unstandardisierten Mittelwertsdifferenzen zwischen Interventionsgruppe und Kontrollgruppe an. Die roten Punkte zeigen die (un-)standardisierten Mittelwertsdifferenzen vor der PS-Gewichtung in der Gesamtstichprobe. Wir sehen vor der Gewichtung starke Imbalance, alle Mittelwertsdifferenzen liegen deutlich außerhalb der Grenzen nach Steiner et al. (2010), also außerhalb der gestrichelten Linien. Die türkisfarbenen Punkte stellen die Mittelwertsunterschiede nach der Gewichtung dar. Diese deuten auf eine gute Balance hin, sie liegen innerhalb der Grenzen nach Steiner et al. (2010).
love.plot(ps_weighted, # weightit-Objekt
stats = "variance.ratios", # Varianzverhältnisse
thresholds = c(0.8, 1.25), # Einzeichnen der Grenzwerte nach
# Steiner et al. (2010)
sample.names = c("ungewichtet", # Namen für Stichproben
"gewichtet")
)
Abbildung 2: Varianzverhältnisse der (kontinuierlichen) Kovariate und Propensity-Scores vor (rote Punkte) und nach (türkisfarbene Punkte) der PS-Gewichtung. Die durchgezogene senkrechte Linie kennzeichnet ein Varianzverhältnis von 1 und somit gleiche Varianzen in der Interventionsgruppe und Kontrollgruppe. Die gestrichelten Linien kennzeichnen die Grenzen für akzeptable Varianzverhältnisse nach Steiner et al. (2010).
Diese Grafik gibt für die Propensity-Scores und die Vorkenntnisse die Varianzverhältnisse zwischen Interventions- und Kontrollgruppe an. Für die Englischvorliebe (eng_vorliebe) werden keine Varianzverhältnisse ausgegeben, da es sich um eine dichotome Variable handelt. Die roten Punkte stellen die Varianzverhältnisse vor der Gewichtung dar. Die Vorkenntnisse (eng_pre) sind vor der Gewichtung nicht hinreichend balanciert, wohingegen die Propensity-Scores (prop.score) vor der Gewichtung innerhalb der Grenzen nach Steiner et al. (2010) liegen und somit als hinreichend balanciert betrachtet werden können. Die türkisfarbenen Punkte stellen die Varianzverhältnisse nach der Gewichtung dar. Zwar liegen die Propensity-Scores nach der Gewichtung nach wie vor innerhalb der Grenzen nach Steiner et al. (2010), die Vorkenntnisse sind jedoch nach wie vor nicht ausreichend balanciert (sie sind sogar schlechter balanciert). Verwendet man die weniger strengen Grenzen nach Stuart und Rubin (2008) liegen die Varianzverhältnisse sowohl vor als auch nach der Gewichtung im akzeptablen Bereich (zwischen 0,5 und 2).
6.2 Verteilung der Kovariaten
Mit dem ´bal.plot´ Befehl können Unterschiede in der gesamten Verteilung der Kovariaten nach der Gewichtung untersucht werden.
bal.plot(ps_weighted, # matchit-Objekt
var.name = "eng_pre", # Variable, für die Kovariatenverteilungen
# dargestellt werden sollen
which = "both", # Kovariatenverteilung vor und nach Matching
type = "histogram", # Darstellung als Histogramm
mirror = TRUE, # Spiegelung der gruppenspezifischen Verteilungen an x-Achse
sample.names = c("ungewichtete Stichprobe", # Namen für Stichproben
"gewichtete Stichprobe"))
Abbildung 3: Verteilung der kontinuierlichen Kovariate eng_pre in der Kontrollgruppe (rote Balken) und in der Interventionsgruppe (türkisfarbene Balken) in der Gesamtstichprobe vor (links) und nach der PS-Gewichtung (rechts).
In Abbildung 3 ist zu sehen, dass die Verteilung der Englischvorkenntnisse sich vor der PS-Gewichtung (links) deutlich zwischen Interventions- und Kontrollgruppe unterscheidet. Personen in der Interventionsgruppe haben beispielsweise höhere Vorkenntnisse als Personen in der Kontrollgruppe. Nach der Gewichtung (rechts) ist die Verteilung der Vorkenntnisse in beiden Interventionsgruppen sehr ähnlich, was für eine gute Balance spricht.
Insgesamt sprechen die Ergebnisse nicht für eine sehr gute, aber doch akzeptable Balance.
7 Schätzung kausaler Effekte
Bei der PS-Gewichtung werden Selektionseffekte kontrolliert, indem einzelne Beobachtungen nicht gleichermaßen in die Analyse einbezogen, sondern entsprechend ihrer Über- oder Unterrepräsentation in der jeweiligen Interventionsgruppe gewichtet werden. Dabei werden Beobachtungen unterrepräsentierter Personen in einer Gruppe stärker und Beobachtungen überrepräsentierter Personen in einer Gruppe weniger stark gewichtet.
Um den Effekt unseres Englisch-Trainings auf die Englischkentnisse zu schätzen, speichern wir zunächst die oben bestimmten Gewichte in unserem Datensatz.
englisch.overlap$ATE_weights <- ps_weighted$weightsDa der durchschnittliche kausale Effekt als Zieleffekt (s.o.
estimand = "ATE") festgelegt wurde, können wir diesen nun
mit dem lm()-Befehl schätzen:
m1 <- lm(eng_post ~ treat, # Regression vom Englisch-Nachtest
# auf die Interventionsvariable
data = englisch.overlap, # Datensatz
weights = ATE_weights) # die geschätzten Gewichte zur Berechnung des ATE
summary(m1) # Ergebnisse##
## Call:
## lm(formula = eng_post ~ treat, data = englisch.overlap, weights = ATE_weights)
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -2.231 -0.547 -0.066 0.513 4.008
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.0091 0.0378 79.67 < 0.0000000000000002 ***
## treat 0.2345 0.0536 4.38 0.000015 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.788 on 429 degrees of freedom
## Multiple R-squared: 0.0427, Adjusted R-squared: 0.0405
## F-statistic: 19.2 on 1 and 429 DF, p-value: 0.0000152Der geschätzte ATE des Englischtrainings kann unter
Coefficients abgelesen werden
und liegt bei 0,235 (SE = 0,054). Personen in
der Region of Common Support haben im Mittel ein um 0,235 Punkte
besseres Ergebnis im Englisch-Test, wenn sie am Training teilnehmen, als
wenn sie nicht am Training teilnehmen. Dieser Effekt ist statistisch
signifikant von null verschieden (\(t(1)\) = 4,380; \(p\) < ,001).
Die hier angegebene Inferenzstatistik berücksichtigt nicht, dass auch die PS-Schätzung aus dem ersten Schritt mit Unsicherheit verbunden ist. Es existieren Verfahren, die diese Unsicherheit in der Inferenz für die kausalen Effekte berücksichtigen.
8 Weiterführende Literatur
Weiterführende Literatur zur Umsetzung von PS-Methoden in R finden Sie hier:
- Leite, W. (2016). Practical propensity score methods using R. Sage Publications.
Weitere Informationen zu den Paketen WeightIt und cobalt sind hier zu finden:
Greifer, N. (2022, 28 Juni). WeightIt: Weighting for Covariate Balance in Observational Studies. NGreifer. https://ngreifer.github.io/WeightIt/.
Greifer, N. (2022, 03. November). Covariate balance tables and plots: A guide to the cobalt package. Cran R- Project. https://cran.r-project.org/web/packages/cobalt/vignettes/cobalt_A0_basic_use.html
Benchmarks für Balance finden Sie beispielsweise hier:
Steiner, P. M., Cook, T. D., Shadish, W. R., & Clark, M. H. (2010). The importance of covariate selection in controlling for selection bias in observational studies. Psychological methods, 15(3), 250. https://doi.org/10.1037/a001871
Stuart, E. A., & Rubin, D. B. (2008). Best practices in quasi-experimental designs: matching methods for causal inference. In J. Osborne (Hrsg.), Best Practices in Quantitative Methods (155-176). Sage Publications. https://doi.org/10.4135/9781412995627