library(foreign)
allbus <- read.spss("./Daten/Allbus2021/ZA5280_v2-0-0.sav",
to.data.frame = TRUE,
use.value.labels = FALSE,
reencode = TRUE)5. Wie nutzen wir einfache Gewichte in R? Das survey Paket
Einstieg
Hallo und herzlich willkommen zum fünften Teil der Videoserie zum Thema “Gewichtung mit R”. Im vorigen Teil vier haben wir schon ein paar Funktion aus dem sjmisc und dem weights Paket kennengelernt, mit denen wir Daten einfach durch die Angabe der Gewichtungsvariable gewichten konnten.
In diesem Teil wollen wir uns nun das survey Paket genauer anschauen. Das survey Paket ist der Standard für die Gewichtung von Daten in R. Es bietet nicht nur die meisten Funktionen, sondern ist am meisten ausgereift, was die Entwicklung betrifft. Der Aufbau der Funktionen ist sehr kohärent und die Integration in einen typischen Analyseablauf gelingt ausgezeichnet.
Der stärkste Vorteil des survey Paketes liegt aber darin, das man das komplexe Sampling-Design angeben kann, aus dem die Daten stammen. Erst damit wird es möglich auch die statistischen Folgen des Stichprobendesigns korrekt in den Analysen zu berücksichtigen - also auch die Varianzen und damit Standardfehler richtig zu schätzen.
Dazu schauen wir uns zuerst die grundlegende Funktionsweise des survey Paketes und seiner Funktionen an. Danach gibt es zwei Beispielen, an denen deutlich wird, wie wichtig die Berücksichtigung des Stichprobendesigns ist, um belastbare Ergebnisse zu erhalten. Ein Beispiel wird die Daten des Allbus nutzen, das andere die des ESS.
Am Ende des Videos sollen Sie in der Lage sein, ein komplexes Stichprobendesign korrekt als survey-Designobjekt in R definieren zu können um damit verschiedene deskriptive, bivariate und multivariate Analysen durchführen zu können. Außerdem sollen Sie aus diesem Video mitnehmen, dass das Stichprobendesign nicht nur die Punktschätzer verzerren kann, sondern auch Auswirkungen auf die Signifikanztests hat, weil die Schätzung der Varianzen und damit der Standardfehler beeinflusst wird.
Wechseln wir also nun in R und schauen und das survey Paket an…
Das survey Paket in R
Bevor wir loslegen können, importieren wir wieder den Allbus und den ESS Datensatz, bzw. beide ESS Datensätze und spielen sie zusammen.
library(foreign)
ess.i <- read.spss("./Daten/ESS10/ESS10.sav",
to.data.frame = TRUE,
use.value.labels = FALSE,
reencode = TRUE)
ess.sc <- read.spss("./Daten/ESS10/ESS10SC.sav",
to.data.frame = TRUE,
use.value.labels = FALSE,
reencode = TRUE)
# Zusammenführen
vars <- names(ess.i)[names(ess.i) %in% names(ess.sc)]
ess <- rbind(ess.sc[,vars], ess.i[,vars])
ess <- ess[!(ess$cntry %in% c("IL", "RS")),]Dann laden wir installieren wir das survey Paket und laden es aus der Bibliothek:
#install.pakckages("survey")
library(survey)Grundlegende Funktionsweise
Man kann das survey Paket für die einfache Gewichtung von Daten mit Hilfe einer Gewichtungsvariable nutzen, oder man kann mit dem Paket zusätzlich zur Gewichtung das komplexe Stichprobendesign von Daten berücksichtigen.
In beiden Fällen ist der Ablauf gleich:
Man definiert zuerst ein sogenanntes survey-Design Objekt mit der Funktion svydesign(). Dieses Objekt enthält die Daten und die Informationen zur Gewichtung und zum Stichprobendesign.
Dann nutzt man - je nach dem, welche Analyse man durchführen möchte - die passende Funktion aus dem survey Paket, und wendet sie auf das survey-Design Objekt an.
Diese Grundlogi schauen wir uns als erstes an einem einfachen Beispiel an. Wie schon in Teil vier wollen wir aus dem Allbus die Sonntagsfrage nach der Wahlabsicht als Tabelle mit gewichteten Daten ausgeben lassen. Das war die Variable pv01. Und die Gewichtungsvariable im Allbus war die Variable wghtpew.
Als erstes definieren wir das survey-Design Objekt mit der Funktion svydesign().
In der einfachsten Form, wenn das survey-Design nur die Gewichtungsvariable enthalten soll, dann sieht der Befehl so aus (wichtig ist, das das survey Paket alle Variablenangaben mit einer Tilde ~ erwartet:1
allbus_svy <- svydesign(ids=~1, # Definiert die Cluster. Wenn ~1 keine Cluster
weights=~wghtpew, # Gewichtungsvariable
data=allbus) # DatensatzJetzt können wir die Funktion svytable() nutzen, um eine Häufigkeitstabelle auf Basis des survey-Design Objektes zu erstellen. Da im survey-Design Objekt die Gewichtung definiert wurde, erhalten wir nun eine gewichtete Häufigkeitstabelle:
svytable(~pv01, design=allbus_svy)pv01
1 2 3 4 6 42 90 91
1073.4533 605.0665 486.5088 1018.8546 252.2724 228.1871 153.6724 212.9900 Das wir hier Dezimalstellen haben, obwohl es absolute Häufigkeiten sind, liegt an der Gewichtung. Natürlich gibt es nicht 1073.4533011 Befragte. Wen das stört: man kann Angeben, auf welche Fallzahl standardisiert werden soll und das gerundet wird:
svytable(~pv01,
design=allbus_svy,
Ntotal=sum(!is.na(allbus$pv01)),
round=T)pv01
1 2 3 4 6 42 90 91
1072 604 486 1018 252 228 153 213 Wenn man auf 100 Befragte standardisiert, erhält man natürlich das selbe Ergebnis, wie wenn man eine relative Häufigkeitstabelle mit prop.table() anfordert:
svytable(~pv01,
design=allbus_svy,
Ntotal=100)pv01
1 2 3 4 6 42 90 91
26.629917 15.010315 12.069169 25.275449 6.258299 5.660799 3.812259 5.283793 prop.table(svytable(~pv01, design=allbus_svy))*100pv01
1 2 3 4 6 42 90 91
26.629917 15.010315 12.069169 25.275449 6.258299 5.660799 3.812259 5.283793 Wir können auch Kreuztabellen erstellen. Hier sehen wir auch, warum das survey Paket die Tilde ~ vor der Variablenangabe erwartet - die Angaben werden als formula Angabe interpretiert. So können wir Interaktionen definieren:
svytable(~pv01+sex, design=allbus_svy) sex
pv01 1 2 3
1 547.4252344 523.2610320 0.0000000
2 318.2593017 285.0534399 0.0000000
3 279.4806222 207.0282094 0.0000000
4 457.9816846 558.3785739 2.4943506
6 134.6130991 115.3988497 0.5066198
42 161.8945168 65.2793475 0.0000000
90 80.3243937 71.5941665 0.0000000
91 96.9266016 116.0633544 0.0000000Hätten wir oben bei der Definition des survey-Design Objekt keine Gewichtung angegeben, hätten wir als Resultat eine ungewichtete Tabelle erhalten (und erhalten deshalb eine Warnung - schließlich nutzen wir extra die svytable() Funktion um gewichten zu können):
allbus_svy_ohne_weight <- svydesign(ids=~1,
data=allbus)
prop.table(svytable(~pv01, design=allbus_svy_ohne_weight))pv01
1 2 3 4 6 42 90
0.25906607 0.14778937 0.11798311 0.23447591 0.07575758 0.06830601 0.04048684
91
0.05613512 Verschiedene Funktionen
Das survey Paket hält eine große Vielzahl an Funktionen für unterschiedliche statistische Verfahren bereit.
Das können wir uns mal in der Beschreibung des Paketes anschauen… Entweder in RStudio, indem wir im Reiter “Pakete” auf das survey Paket klicken, oder inde, wir die pdf Dokumentation aufrufen: https://cran.r-project.org/web/packages/survey/survey.pdf
Das Paket enthält viele Funktionen, die wir überhaupt nicht direkt nutzen wollen, oder die sehr spezielle Probleme lösen. Was für uns wichtig ist, sind die Funktionen, die mit svy...() beginnen. Also scrollen wir mal zum Buchstaben “S”…
Hier sehen Sie eine Reihe von Funktionen, die ganz unterschiedliche statistische Tests mit dem survey-Design Objekt durchführen können, die also die im survey-Design Objekt angegebene Gewichtung und das Stichprobendesign berücksichtigen.
Zum Beispiel hier, einen \(\chi^2\) Test mit der Funktion svychisq() oder hier einen gewichteten Mittelwert mir der Funktion svymean(), oder auch einen t-Test mit svyttest().
Alle funktionieren nach der gleichen Logik. Schauen wir uns nochmal ein einfaches Beispiel an, der Mittelwert für Religiösität im Allbus (Variable rb07, 1 ‘nicht religiös’ … 10 ‘religiös’):
svymean(~rb07, allbus_svy, na.rm=T) mean SE
rb07 4.2915 0.0539Und jetzt machen wir mal einen Mittelwertvergleich mittels t-Test, ob die Bayern tatsächlich religiöser sind als der Rest von Deutschland. Dazu bereite ich kurz eine dichotome Variable vor, Bayern ja/nein:
# Bayern ist Merklamsausprägung 90
attr(allbus$land,"value.labels") THUERINGEN SACHSEN-ANHALT SACHSEN
"160" "150" "140"
MECKLENB.-VORPOMMERN BRANDENBURG EHEM. BERLIN-OST
"130" "120" "112"
EHEM. BERLIN-WEST SAARLAND BAYERN
"111" "100" "90"
BADEN-WUERTTEMBERG RHEINLAND-PFALZ HESSEN
"80" "70" "60"
NORDRHEIN-WESTFALEN BREMEN NIEDERSACHSEN
"50" "40" "30"
HAMBURG SCHLESWIG-HOLSTEIN
"20" "10" # Gruppierungsvariable vorbereiten:
allbus$bayern <- car::recode(allbus$land,
"10:80 = FALSE;
90 = TRUE;
100:160 = FALSE;
else = NA")Und diese nutze ich nun als Gruppierungsvariable in meiner svyttest() Funktion:
svyttest(rb07~bayern, allbus_svy, na.rm=T)Tja, und wie geplant: Der Code läuft nicht, sondern produziert eine Fehlermeldung! Warum ist das so: Wir haben die Variable bayern im Objekt allbus erstellt. Wir arbeiten in der svyttest() Funktion mit dem zuvor erstellten survey-Design Objekt. In diesem ist die bayern Variable ja gar nicht enthalten.
Also wir merken uns: wenn Änderungen an den Variablen durchgeführt werden, muss das survey-Design Objekt aktualisiert werden!
allbus_svy <- svydesign(ids=~1, # Definiert die Cluster. Wenn ~1 keine Cluster
weights=~wghtpew, # Gewichtungsvariable
data=allbus) # Datensatz
svyttest(rb07~bayern, allbus_svy, na.rm=T)
Design-based t-test
data: rb07 ~ bayern
t = 5.2395, df = 3553, p-value = 1.703e-07
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
0.4890504 1.0739078
sample estimates:
difference in mean
0.7814791 Regression mit dem survey Paket
Kommen wir nun noch zu einer letzten svy...() Funktion, die ich Ihnen zeigen möchte.
Es handelt sich um den Regressionsbefehl.
Regressionen mit gewichteten Daten lassen sich auch mit der base-R Funktion lm() für linear model, also für die OLS Regression berechnen.
Hier sehen wir den OLS Regressionsbefehl, mit Gewichtung um die Variable wghtpew.
fit.lm <- lm(ps03 ~ pt03 + isced97 + incc + age + eastwest,
data=allbus,
weight=wghtpew)
summary(fit.lm)
Call:
lm(formula = ps03 ~ pt03 + isced97 + incc + age + eastwest, data = allbus,
weights = wghtpew)
Weighted Residuals:
Min 1Q Median 3Q Max
-2.8220 -0.6086 -0.1136 0.4918 3.7924
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.408099 0.105812 41.660 < 2e-16 ***
pt03 -0.393719 0.011370 -34.627 < 2e-16 ***
isced97 -0.025499 0.017216 -1.481 0.139
incc -0.014828 0.003561 -4.164 3.22e-05 ***
age -0.004984 0.000996 -5.004 5.95e-07 ***
eastwest 0.320186 0.045058 7.106 1.49e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.9069 on 2918 degrees of freedom
(2418 observations deleted due to missingness)
Multiple R-squared: 0.3406, Adjusted R-squared: 0.3395
F-statistic: 301.5 on 5 and 2918 DF, p-value: < 2.2e-16Regressionen mit gewichteten Daten lassen sich aber auch im survey Paket umsetzen. Dazu nutzen wir die Funktion svyglm():
fit <- svyglm(ps03 ~ pt03 + isced97 + incc + age + eastwest,
design=allbus_svy)
summary(fit)
Call:
svyglm(formula = ps03 ~ pt03 + isced97 + incc + age + eastwest,
design = allbus_svy)
Survey design:
svydesign(ids = ~1, weights = ~wghtpew, data = allbus)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.408099 0.121524 36.273 < 2e-16 ***
pt03 -0.393719 0.013667 -28.807 < 2e-16 ***
isced97 -0.025499 0.018322 -1.392 0.164
incc -0.014828 0.003734 -3.971 7.32e-05 ***
age -0.004984 0.001065 -4.678 3.03e-06 ***
eastwest 0.320186 0.039441 8.118 6.91e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.8210609)
Number of Fisher Scoring iterations: 2\(R^2\) berechnen mit Hilfe des performance Paketes:
performance::r2(fit)# R2 for Survey Regression
R2: 0.341
adj. R2: 0.339Wir wir sehen, sind die Ergebnisse identisch zu den Ergebnissen der Regression mit dem lm() Befehl.
Sie werden schon gemerkt haben, dass ist nicht das klassische OLS Modell, dass wir in R mit der Funktion lm() für linear model aufrufen. Stattdessen ist es ein sogenanntes Generalisiertes Lineares Modell GLM. Aber auch svyglm() schätzt, wenn man keine weiteren Argumente angibt, eine lineare Regression - nur eben mit dem Maximum Likelihood Schätzer.
Wir können die svyglm() aber auch nutzen um eine Vielzahl anderer Modelle zu schätzen. Zum Beispiel logistische Regression, Multinomiale Regression oder Poisson-Regression. Man muss nur die passende Link-Funktion angeben. Wer dazu mehr wissen will, schaut sich am besten die Hilfe zur svyglm() Funktion an:
?svyglmSIe werden sich jetzt sicher Fragen: Wenn der lm() Befehl die gleichen Ergebnisse liefert, wie die komplizierte svyglm() Funktion? Warum sollten wir sie dann nutzen?
Das ist eine sehr gute Frage. Aber die noch bessere Antwort gibt es im nächsten Teil 6.
Literatur
- Manderscheid, K., 2017. Sozialwissenschaftliche Datenanalyse mit R. Springer Fachmedien Wiesbaden, Wiesbaden. Kap. 11, S. 225-238. https://doi.org/10.1007/978-3-658-15902-3
Fußnoten
Als weiterführende Info: das ist deshalb der Fall, da die Angaben als Formel interpretiert werden und so auch Interaktionen definiert werden können↩︎