7 Inferenzstatistik & Zusammenhangsmaße

library(tidyverse)

etb18_kap7 <- haven::read_dta("./data/BIBBBAuA_2018_suf1.0.dta",
                         n_max = 1000,
                         col_select = c("intnr","zpalter","S1","F231","az","F600_12",
                                        "m1202","F204","gew2018_hr17","gew2018","F518_SUF",
                                        "F411_01","S3","Mig")) %>% 
  mutate(F600_12 = ifelse(F600_12 > 4,NA,F600_12),
         zpalter = ifelse(zpalter>100,NA,zpalter),
         m1202 = ifelse(m1202 > 4,NA,m1202))

Bisher haben wir die Angaben aus unserem Datensatz immer als fix angesehen. Ziel einer statistischen Auswertung ist aber meistens, Aussagen über die Grundgesamtheit oder Population zu treffen. Im Fall der ETB 2018 wären das also alle Erwerbstätigen in Deutschland.

Für die grundlegenden inferenzstatistischen Tests helfen uns im wesentlichen 3 Funktionen:

t.test() zur Berechnung von t-Tests
cor.test() zur Berechnung von Korrelationen und weiterer Zusammenhangsmaße
{effectsize} zur Berechnung von bspw. Cohen’s D

Neue Pakete:

install.packages("effectsize") # u.a. Cohen's D berechnen 
install.packages("correlation") # Korrelationmatrix erstellen
install.packages("survey") # Gewichtung (weiter unten)

7.1 t-Tests

Eines der zentralen Werkzeuge der grundlegenden Inferenzstatistik ist der t-Test. In R steht uns dafür t.test() zur Verfügung. Mit der Option mu = können wir einen Wert für die $H_{A}$ angeben:

t.test(etb18_kap7$zpalter, mu = 45)


    One Sample t-test

data:  etb18_kap7$zpalter
t = 9.1026, df = 995, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 45
95 percent confidence interval:
 47.50841 48.88717
sample estimates:
mean of x 
 48.19779

Ein weiterer typischer Anwendungsfall für t-Tests ist der Gruppenvergleich, dazu geben wir in t.test() die zu testende Variable und und nach einer ~¹ die Gruppenvariable an. Wir testen hier also auf Altersunterschiede zwischen Männern (S1=1, daher group1) und Frauen (S1=2, daher group2).

t.test(etb18_kap7$zpalter~etb18_kap7$S1)


    Welch Two Sample t-test

data:  etb18_kap7$zpalter by etb18_kap7$S1
t = -2.549, df = 974.9, p-value = 0.01095
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 -3.1679965 -0.4119377
sample estimates:
mean in group 1 mean in group 2 
       47.28484        49.07480

Es hat sich als Konvention etabliert, von einem signifikanten Unterschied zu sprechen wenn die Irrtumswahrscheinlichkeit unter 5% liegt. Das bedeutet:

Da hier der p-Wert sehr viel kleiner ist als 0.05 ist (p-value < 2.2e-16)², können wir von einen statistisch signifikanten Unterschied sprechen.

Standardmäßig bekommen wir einen beidseitigen Test ("two.sided"), wir können aber auch einen links- ("less") oder rechtsseitigen ("greater") Test anfordern mehr dazu:

t.test(etb18_kap7$zpalter~etb18_kap7$S1,alternative = "two.sided")
t.test(etb18_kap7$zpalter~etb18_kap7$S1,alternative = "less")
t.test(etb18_kap7$zpalter~etb18_kap7$S1,alternative = "greater")

Für Effektgrößen wie Cohen’s D empfiehlt sich das Paket {effectsize}:

library(effectsize)

Warning: Paket 'effectsize' wurde unter R Version 4.2.2 erstellt

cohens_d(etb18_kap7$zpalter~etb18_kap7$S1)

Warning: Missing values detected. NAs dropped.

Cohen's d |         95% CI
--------------------------
-0.16     | [-0.29, -0.04]

- Estimated using pooled SD.

7.1.1 Übung

7.2 Korrelation

Den Korrelationskoeffizienten können wir in R mit cor.test() berechnen:

F231: Wie viele Stunden arbeiten Sie i.d.R. im Durchschnitt pro Woche von zu Hause aus?

cor.test(etb18_kap7$zpalter,etb18_kap7$F231,method = "pearson")


    Pearson's product-moment correlation

data:  etb18_kap7$zpalter and etb18_kap7$F231
t = 0.48542, df = 286, p-value = 0.6277
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.08717658  0.14379451
sample estimates:
       cor 
0.02869194

Es handelt sich mit 0.0287 also um einen geringen Zusammenhang. Der p-Wert gibt uns auch hier wieder Auskunft über die stat. Signifikanz: mit 0.62775 liegt der p-Wert deutlich unter 0,05.

Um eine Korrelationsmatrix zu erhalten hilft das Paket {correlation}:

install.packages("correlation")

library(correlation)

Warning: Paket 'correlation' wurde unter R Version 4.2.2 erstellt

etb18_kap7 %>% select(zpalter,F231,az) %>% 
  correlation() %>% 
  summary(.)

# Correlation Matrix (pearson-method)

Parameter |    az | F231
------------------------
zpalter   | -0.05 | 0.03
F231      | -0.10 |     

p-value adjustment method: Holm (1979)

corr_df <- 
  etb18_kap7 %>%
    group_by(S1) %>%
    select(zpalter,F231,az) %>% 
    correlation()

Adding missing grouping variables: `S1`

data.frame(corr_df)

  Group Parameter1 Parameter2            r   CI      CI_low    CI_high
1     1    zpalter       F231 -0.008016615 0.95 -0.16106862 0.14541192
2     1    zpalter         az -0.006654458 0.95 -0.09536128 0.08215722
3     1       F231         az -0.107424635 0.95 -0.25645259 0.04659229
4     2    zpalter       F231  0.029864101 0.95 -0.14722774 0.20510074
5     2    zpalter         az -0.038373166 0.95 -0.12495280 0.04878647
6     2       F231         az  0.042311811 0.95 -0.13429355 0.21631207
           t df_error         p              Method n_Obs
1 -0.1020381      162 1.0000000 Pearson correlation   164
2 -0.1467035      486 1.0000000 Pearson correlation   488
3 -1.3752506      162 0.5128577 Pearson correlation   164
4  0.3300070      122 1.0000000 Pearson correlation   124
5 -0.8638192      506 1.0000000 Pearson correlation   508
6  0.4696813      123 1.0000000 Pearson correlation   125

7.3 Weitere Zusammenhangsmaße:

7.3.1 Rangkorrelation & Kendall’s $\tau$

Ein klassisches ordinales Merkmal ist der höchste Ausbildungsabschluss in m1202. Wir sehen uns den (möglichen) Zusammenhang zwischen dem höchsten Ausbildungsabschluss m1202 und F600_12 an:

v	l
F600_12	Häufigkeit: unter Lärm arbeiten
1	häufig
2	manchmal
3	selten
4	nie

Für den Spearman-Rangkorrelationskoeffizienten können wir method = "spearman" nutzen:

cor.test(etb18_kap7$m1202,etb18_kap7$F600_12,method = "spearman")
cor.test(etb18_kap7$m1202,etb18_kap7$F600_12,method = "spearman")

Ein weiteres Zusammenhangsmaß für ordinale Variablen sind Konkordanzmaße wie Kendall’s $\tau$, mit der Option method = "kendall":

cor.test(etb18_kap7$m1202,etb18_kap7$F600_12, method = "kendall")

Hier wird Kendall’s $\tau_b$ ausgegeben - siehe Anhang für Kendall’s $\tau_a$

7.3.2 $\chi^2$ & Cramér’s $v$

Mit chisq.test() bekommen wir $\chi^2$ ausgegeben, diese müssen wir jedoch auf ein xtabs()-Objekt anwenden (F204 Mehrarbeitsvergütung und S1 Geschlecht)

tab1 <- xtabs(~ F204 + S1, data = etb18_kap7)
chisq.test(tab1)


    Pearson's Chi-squared test

data:  tab1
X-squared = 12.541, df = 4, p-value = 0.01375

Für Cramér’s $v$ können wir wieder auf {effectsize} zurückgreifen:

library(effectsize)
cramers_v(etb18_kap7$F204,etb18_kap7$S1)

Cramer's V (adj.) |       95% CI
--------------------------------
0.14              | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

Weitere Maße

7.3.3 Übung

7.4 Gewichtung

Bei der Datenanalyse ist man oft mit einer Stichprobe aus einer größeren Population konfrontiert und man möchte aus dieser Stichprobe Rückschlüsse auf die größere Population ziehen. Die meisten statistischen Verfahren für diese sog. „Inferenzstatistik“ beruhen dabei auf der Annahme, dass die Stichprobe eine einfache Zufallsstichprobe ist. Bei einer solchen Stichprobe gelangen alle Elemente der Grundgesamtheit mit der gleichen Wahrscheinlichkeit in die Stichprobe. In der Praxis sind solche Stichproben aber die große Ausnahme. Häufig haben bestimmte Gruppen von Personen höhere Auswahlwahrscheinlichkeiten als andere. Kohler/Kreuter, S.81

Gewichte sind ein häufig verwendetes Gegenmittel. Die einfachste Variante für eine Gewichtung ist die Option wt= in count():

etb18_kap7 %>% 
  count(S1,m1202,wt = gew2018)

# A tibble: 10 × 3
   S1           m1202      n
   <dbl+lbl>    <dbl>  <dbl>
 1 1 [männlich]    -1   1.06
 2 1 [männlich]     1  61.6 
 3 1 [männlich]     2 266.  
 4 1 [männlich]     3  68.0 
 5 1 [männlich]     4 144.  
 6 2 [weiblich]    -1   1.01
 7 2 [weiblich]     1  37.4 
 8 2 [weiblich]     2 246.  
 9 2 [weiblich]     3  29.8 
10 2 [weiblich]     4 127.

Für umfangreichere Anwendungen stehen in R stehen die Pakete {survey} und das darauf aufbauende {srvyr} zur Verfügung.

install.packages("survey")

Zunächst verwenden wir svydesign(), um die Gewichtung festzulegen. Im weiteren stellt {survey} dann zahlreiche Funktionen zur Verfügung, die eine gewichtete Variante der basis-Befehle sind - bspw. svymean() und svytable():

library(survey)
etb18_kap7_weighted <- svydesign(id      = ~intnr,
                            weights = ~gew2018,
                            data    = etb18_kap7)

svymean(~zpalter, etb18_kap7_weighted, na.rm = TRUE)

          mean     SE
zpalter 44.234 0.6175

mean(etb18_kap7$zpalter, na.rm = TRUE)

[1] 48.19779

Für Tabellen gibt es in {survey} auch eine Funktion:

svytable(~S1+m1202,etb18_kap7_weighted)

   m1202
S1       -1       1       2       3       4
  1   1.063  61.553 265.619  68.011 143.747
  2   1.012  37.438 245.726  29.845 126.756

xtabs(~S1+m1202,etb18_kap7)

   m1202
S1   -1   1   2   3   4
  1   1  32 221  54 181
  2   2  26 253  41 189

Für Regressionsmodelle gibt es bspw. survey::svyglm()

7.5 Übungen

7.5.1 Übung 1

etb18_ue7 <- haven::read_dta("./data/BIBBBAuA_2018_suf1.0.dta",
                         n_max = 1000,
                         col_select = c("intnr","S1","az","F518_SUF","m1202","F411_01","gew2018")) %>% 
  mutate(across(matches("m1202|F411_01"), ~ifelse(.x > 4|.x<0,NA,.x)),
         F518_SUF = ifelse(F518_SUF>99990,NA,F518_SUF)) # missings in m1202 mit NA überschreiben

Testen Sie die Hypothese, dass ein signifikanter Unterschied in der Arbeitszeit (az) zwischen Männern und Frauen besteht (S1). Beide Variablen haben keine Missings - Sie können direkt los legen.

Sehen Sie sich die Informationen an, welche t.test() erstellt: legen Sie das Ergebnis des Tests als Objekt ab (test1 <- ...) und sehen Sie sich die Informationen unter test1$ an (drücken Sie die ↹ Taste).

Berechnen Sie das Cohen’s d für diesen Zusammenhang.

7.5.2 Übung 2

Untersuchen Sie die Korrelation zwischen der Wochenarbeitszeit (az) und dem Einkommen (F518_SUF) der Befragten.
Die Missings in F518_SUF sind mit dem Befehl oben schon ausgeschlossen, az hat keine Missings.
Berechnen Sie die Rangkorrelation für den Zusammenhang zwischen der Häufigkeit von starkem Termin- oder Leistungsdruck F411_01 und der Ausbildungsvariable m1202.

var	.	value	name
F411_01	Wie häufig kommt es vor, dass Sie unter starkem Termin- oder Leistungsdruck arbe	1	häufig
		2	manchmal
		3	selten
		4	nie
		9	keine Angabe
m1202	Höchster Ausbildungsabschluss	-1	keine Angabe
		1	Ohne Berufsabschluss
		2	duale o. schulische Berufsausbildung/einf.,mittl. Beamte
		3	Aufstiegsfortbildung (Meister, Techniker, kfm. AFB u.ä.)
		4	Fachhochschule, Universität/ geh., höhere Beamte

7.5.3 Übung 3

Legen Sie die Gewichtung auf Basis von gew2018 an.
Berechnen Sie den Mittelwert für den Bruttoverdienst F518_SUF mit und ohne Gewichtung.

7.6 Anhang

7.6.1 Weitere Korrelationsmaße

Kendall’s $\tau_a$, welches im Nenner alle Paarvergleiche berücksichtigt, können wir mit KendallTauA() aus dem Paket DescTools berechnen:

install.packages("DescTools")

library(DescTools)
KendallTauA(etb18_kap7$m1202,etb18_kap7$F600_12)

[1] 0.08970494

KendallTauB(etb18_kap7$m1202,etb18_kap7$F600_12) # entspricht in method = "kendall" cor.test()

[1] 0.1367877

cor.test(etb18_kap7$m1202,etb18_kap7$F600_12,method = "kendall")


    Kendall's rank correlation tau

data:  etb18_kap7$m1202 and etb18_kap7$F600_12
z = 4.9325, p-value = 8.118e-07
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.1367877

Der Wert von Kendall’s $\tau_a$ ist niedriger als von Kendall’s $\tau_b$, da hier der Nenner durch die Berücksichtigung aller möglichen Paarvergleiche größer wird, der Zähler aber für beide Varianten von Kendall’s $\tau$ gleich definiert ist.

Eine andere Alternative, welche auch beim Vorliegen von Bindungen den vollen Wertebereich [-1,1] erreichen kann ist Goodman & Kruskal’s $\gamma$. Dieses können wir mit dem Befehl GoodmanKruskalGamma aus dem Paket {DescTools} berechnen:

library(DescTools)
GoodmanKruskalGamma(etb18_kap7$m1202,etb18_kap7$F600_12)

[1] 0.2064219

Auch Goodman & Kruskal’s $\gamma$ deutet auf einen negativen Zusammenhang hin, hier ist die Stärke aber deutlich höher. Dies ist auf die Berücksichtigung der Bindungen zurückzuführen: hier werden alle Bindungen ausgeschlossen, also auch Paarvgleiche mit Bindungen nur auf einer Variable. Es reduziert sich also der Nenner, somit ergibt sich im Ergebnis ein höherer Koeffizient für Goodman & Kruskal’s $\gamma$ als für Kendall’s $\tau_b$.

Hier findet sich eine Liste weiterer Kennzahlen, die sich mit {correlate} berechnen lassen.

Tastaturbefehle: Alt Gr + * auf Windows. Auf macOS Alt + N und anschließend ein Leerzeichen, damit die Tilde erscheint.↩︎
2.2e-16 steht für 2.2 aber mit 16 Nullen vorweg. Das ist Rs Art zu sagen, dass der Wert sehr sehr klein ist.↩︎