7 Inferenzstatistik & Gewichtung

library(tidyverse)

pend_kap7 <- haven::read_dta("./orig/PENDDAT_cf_W13.dta",
                         col_select = c("pnr","welle","palter","zpsex","azges1","statakt")) %>% 
  mutate(across(everything(), ~ifelse(.x < 0,NA,.x)) )

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.

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

t.test() zur Berechnung von t-Tests
cor.test() zur Berechnung von Korrelationen und weiterer Zusammenhangsmaße
{survey} zur Gewichtung

Neue Pakete:

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(pend_kap7$azges1, mu = 35)


    One Sample t-test

data:  pend_kap7$azges1
t = -10.612, df = 8368, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 35
95 percent confidence interval:
 33.57846 34.02173
sample estimates:
mean of x 
  33.8001

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 (zpsex=1, daher group1) und Frauen (zpsex=2, daher group2).

t.test(pend_kap7$azges1~pend_kap7$zpsex)


    Welch Two Sample t-test

data:  pend_kap7$azges1 by pend_kap7$zpsex
t = 36.395, df = 8318.8, p-value < 2.2e-16
alternative hypothesis: true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
 7.206945 8.027473
sample estimates:
mean in group 1 mean in group 2 
       37.74659        30.12938

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

Assuming that the null hypothesis is true and the study is repeated an infinite number times by drawing random samples from the same populations(s), less than 5% of these results will be more extreme than the current result.²

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(pend_kap7$palter~pend_kap7$zpsex,alternative = "two.sided")
t.test(pend_kap7$palter~pend_kap7$zpsex,alternative = "less")
t.test(pend_kap7$palter~pend_kap7$zpsex,alternative = "greater")

7.2 Korrelation

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

cor.test(pend_kap7$palter,pend_kap7$azges1,method = "pearson")


    Pearson's product-moment correlation

data:  pend_kap7$palter and pend_kap7$azges1
t = -3.8563, df = 8367, p-value = 0.000116
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.06348945 -0.02071526
sample estimates:
        cor 
-0.04212165

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

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

cor.test(pend_kap7$palter,pend_kap7$azges1,method = "spearman")


    Spearman's rank correlation rho

data:  pend_kap7$palter and pend_kap7$azges1
S = 1.0003e+11, p-value = 0.02845
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
-0.02394976

7.2.1 Übung

7.3 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 Mittel, ungleich verteilten Auswahlwahrscheinlichkeiten zu begegnen. Die Gewichte für den Personendatz des PASS CampusFiles sind im Datensatz pweights_cf_W13.dta zu finden. Für eine Gewichtung müssen wir diesen erst an die Personendatei ranspielen - mit einem left_join()

wgt_df <- haven::read_dta("./orig/pweights_cf_W13.dta")
head(wgt_df)

# A tibble: 6 × 7
         pnr welle                   sample             wqp ppbleib   psu strpsu
       <dbl> <dbl+lbl>               <dbl+lbl>        <dbl>   <dbl> <dbl>  <dbl>
1 1000001901 1 [Welle 1 (2006/2007)] 2 [Microm-Stic… 20075.   NA       80     64
2 1000001901 3 [Welle 3 (2008/2009)] 2 [Microm-Stic… 32959.   NA       80     64
3 1000001902 1 [Welle 1 (2006/2007)] 2 [Microm-Stic… 23077.   NA       80     64
4 1000002001 1 [Welle 1 (2006/2007)] 1 [BA-Stichpro…  1639.    1.87    80     64
5 1000002001 2 [Welle 2 (2007/2008)] 1 [BA-Stichpro…  1668.    1.12    80     64
6 1000002001 3 [Welle 3 (2008/2009)] 1 [BA-Stichpro…  1748.    1.36    80     64

pend_kap7w <- pend_kap7 %>% left_join(wgt_df, by = join_by(pnr,welle))

Die einfachste Variante für eine Gewichtung ist die Option wt= in count():

pend_kap7w %>% 
  count(zpsex,statakt)

# A tibble: 8 × 3
  zpsex statakt     n
  <dbl>   <dbl> <int>
1     1       1  4685
2     1       2  3240
3     1       3  2047
4     1      NA  3555
5     2       1  4785
6     2       2  2899
7     2       3  3434
8     2      NA  3779

pend_kap7w %>% 
  count(zpsex,statakt,wt = wqp)

# A tibble: 8 × 3
  zpsex statakt          n
  <dbl>   <dbl>      <dbl>
1     1       1 206326062.
2     1       2  26678834.
3     1       3  64867832.
4     1      NA 143378863.
5     2       1 177223245.
6     2       2  21205451.
7     2       3 109600208.
8     2      NA 163820124.

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)
pend_kap7_wgt <- svydesign(id      = ~pnr,
                           weights = ~wqp,
                           data    = pend_kap7w)

svymean(~palter, pend_kap7_wgt, na.rm = TRUE)

         mean     SE
palter 51.171 0.5755

mean(pend_kap7$palter, na.rm = TRUE)

[1] 46.49612

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

svytable(~zpsex+statakt,pend_kap7_wgt)

     statakt
zpsex         1         2         3
    1 206326062  26678834  64867832
    2 177223245  21205451 109600208

table(pend_kap7$zpsex,pend_kap7$statakt)

   
       1    2    3
  1 4685 3240 2047
  2 4785 2899 3434

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

7.3.1 Übung

7.4 Übungen

7.4.1 Übung 1

pend_ue7 <- haven::read_dta("./orig/PENDDAT_cf_W13.dta",
                            col_select = c("pnr","welle","zpsex","palter","netges")) %>% 
  mutate(across(everything(), ~ifelse(.x < 0,NA,.x))) # missings  mit NA überschreiben

Testen Sie die Hypothese, dass ein signifikanter Unterschied im Alter (palter) zwischen Männern und Frauen besteht (zpsex). Die Missings beider Variablen wurden im mutate() der Befehle oben schon mit NA überschrieben - 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 die Korrelation zwischen Alter (palter) und Einkommen (netges).

7.4.2 Übung 2

wgt_df <- haven::read_dta("./orig/pweights_cf_W13.dta")
head(wgt_df)

# A tibble: 6 × 7
         pnr welle                   sample             wqp ppbleib   psu strpsu
       <dbl> <dbl+lbl>               <dbl+lbl>        <dbl>   <dbl> <dbl>  <dbl>
1 1000001901 1 [Welle 1 (2006/2007)] 2 [Microm-Stic… 20075.   NA       80     64
2 1000001901 3 [Welle 3 (2008/2009)] 2 [Microm-Stic… 32959.   NA       80     64
3 1000001902 1 [Welle 1 (2006/2007)] 2 [Microm-Stic… 23077.   NA       80     64
4 1000002001 1 [Welle 1 (2006/2007)] 1 [BA-Stichpro…  1639.    1.87    80     64
5 1000002001 2 [Welle 2 (2007/2008)] 1 [BA-Stichpro…  1668.    1.12    80     64
6 1000002001 3 [Welle 3 (2008/2009)] 1 [BA-Stichpro…  1748.    1.36    80     64

pend_ue7 <- pend_ue7 %>% left_join(wgt_df, by = join_by(pnr,welle))

Legen Sie die Gewichtung auf Basis von wqp an.
Berechnen Sie den Mittelwert für den Nettoverdienst netges mit und ohne Gewichtung.

Tastaturbefehle: Alt Gr + * auf Windows. Auf macOS Alt + N und anschließend ein Leerzeichen, damit die Tilde erscheint.↩︎
Failing Grade: 89% of Introduction-to-Psychology Textbooks That Define or Explain Statistical Significance Do So Incorrectly. Advances in Methods and Practices in Psychological Science, 2515245919858072.↩︎
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.↩︎