Musterlösung zur Hausaufgabe über die Winterferien

Aufgabenstellung hier

Die angegebene Bepunktung ist vorläufig!

Aufgabe 1

Die benötigte Dateien DPQ_J.xpt findet sich auf der NHANES-Webseite, hier.

Wir laden zunächst Tidyverse, sowie das haven-Paket, dass zum Laden der Dateien im XPT-Format benötigt wird.

library( tidyverse )
library( haven )

Nun laden wir die Tabelle:

read_xpt( "data/DPQ_J.xpt" ) -> dpq

head(dpq)

# A tibble: 6 × 11
   SEQN DPQ010 DPQ020 DPQ030 DPQ040 DPQ050 DPQ060 DPQ070 DPQ080 DPQ090 DPQ100
  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
1 93705      0      0      0      0      0      0      0      0      0     NA
2 93706      0      0      0      0      0      0      0      0      0     NA
3 93708      0      0      0      0      0      0      0      0      0     NA
4 93709     NA     NA     NA     NA     NA     NA     NA     NA     NA     NA
5 93711      1      0      1      0      0      0      0      0      0      0
6 93712      0      0      1      0      0      0      0      0      0      0

Die Tabelle hat folgende Spalten:

colnames(dpq)

 [1] "SEQN"   "DPQ010" "DPQ020" "DPQ030" "DPQ040" "DPQ050" "DPQ060" "DPQ070"
 [9] "DPQ080" "DPQ090" "DPQ100"

Zunächst wandeln wir die Tabelle in ein langes Format um, dann entfernen wir die nicht benötigte Frage 10, und anschließen ersetzen wir alle Antworten mit Codes > 3 durch NA.

dpq %>%
pivot_longer( -SEQN, names_to="question", values_to="answer" ) %>%
filter( question != "DPQ100" ) %>%
mutate( answer = ifelse( answer>3, NA, answer ) ) -> dpq_long

head( dpq_long, 15 )

# A tibble: 15 × 3
    SEQN question answer
   <dbl> <chr>     <dbl>
 1 93705 DPQ010        0
 2 93705 DPQ020        0
 3 93705 DPQ030        0
 4 93705 DPQ040        0
 5 93705 DPQ050        0
 6 93705 DPQ060        0
 7 93705 DPQ070        0
 8 93705 DPQ080        0
 9 93705 DPQ090        0
10 93706 DPQ010        0
11 93706 DPQ020        0
12 93706 DPQ030        0
13 93706 DPQ040        0
14 93706 DPQ050        0
15 93706 DPQ060        0

Nun können wir die die PUnktzahlen für jeden Probanden aufaddieren

dpq_long %>%
group_by( SEQN ) %>%
summarise( dpq_score = sum(answer) ) %>%
rename( subjectId = SEQN ) -> dpq_scores

head( dpq_scores )

# A tibble: 6 × 2
  subjectId dpq_score
      <dbl>     <dbl>
1     93705         0
2     93706         0
3     93708         0
4     93709        NA
5     93711         2
6     93712         1

Bei der Gelegenheit haben wir auch die Spalte SEQN in SubjectID umbenannt

Bepunktung

• Laden der Tabelle: 1 Punkt
• pivot_longer: 1 Punkte
• Frage 10 korrekt entfernt: 1 Punkt
• group_by und summarise korrekt: 1 Punkt
• Code erklärt (zumindest wie die NAs behandelt wurden): 1 Punkt

Häufige Fehler

• Die Antwort für DPQ100 durch NA zu ersetzen ist nicht korrekt, da die NAs dann das sum “infiziert”, es sei denn, man entfern alle Zeilen mit NA.
• Die Zeilen mit NA zu entfernen ist aber auch nicht korrekt, dannd ann werden nicht beantwortete Fragen als 0 gezählt, statt (wie es sein sollte) zu eine Punktsumme von NA zu führen

Aufgabe 2

Wir laden die Datei DEMO_J.XPT, die man hier auf den NHANES-Webseiten findet. Dabei wählen wir gleich die Spalten aus, die wir benötigen:

read_xpt( "data/DEMO_J.XPT" ) %>%
select( subjectId = SEQN, gender_code = RIAGENDR, age = RIDAGEYR, marital_status_code = DMDMARTL ) -> demo

head(demo)

# A tibble: 6 × 4
  subjectId gender_code   age marital_status_code
      <dbl>       <dbl> <dbl>               <dbl>
1     93703           2     2                  NA
2     93704           1     2                  NA
3     93705           2    66                   3
4     93706           1    18                  NA
5     93707           1    13                  NA
6     93708           2    66                   1

Nun erstellen wir eine Tabelle mit den Gender-Codes, wobei wir die Beduutung der Codes dem Codebuch hier entnehmen.

tibble( 
  gender_code = c( 1, 2 ),
  gender = c( "male", "female" )
) -> gender_codes

gender_codes

# A tibble: 2 × 2
  gender_code gender
        <dbl> <chr> 
1           1 male  
2           2 female

Die Codes für marital status finden wir hier. Wir übertragen die Codes:

tibble(
  marital_status_code = 1:6,
  marital_status = c( "married", "widowed", "divorced", "separated", 
      "never married", "living with partner" ) ) -> marital_status_codes
  
marital_status_codes  

# A tibble: 6 × 2
  marital_status_code marital_status     
                <int> <chr>              
1                   1 married            
2                   2 widowed            
3                   3 divorced           
4                   4 separated          
5                   5 never married      
6                   6 living with partner

Die Codes für “refused” (77) und “don’t know” (99) brauchen wir nicht, da sie ohnehin als NA gewertet werden sollen.

Wir fügen die beiden Code-Tabellen, sowie die DPQ-Punktsummen, der demo-Tabelle hinzu. Danach entfernen wir die nun nicht mehr benötigten Spalten mit den Codes:

demo %>%
left_join( gender_codes, by="gender_code", relationship="many-to-one" ) %>%
left_join( marital_status_codes, by="marital_status_code", relationship="many-to-one" ) %>%
left_join( dpq_scores, by="subjectId", relationship="one-to-one" ) %>%
select( -gender_code, -marital_status_code ) %>%
filter( age >= 18 ) -> full_tbl

head( full_tbl )

# A tibble: 6 × 5
  subjectId   age gender marital_status dpq_score
      <dbl> <dbl> <chr>  <chr>              <dbl>
1     93705    66 female divorced               0
2     93706    18 male   <NA>                   0
3     93708    66 female married                0
4     93709    75 female widowed               NA
5     93711    56 male   married                2
6     93712    18 male   <NA>                   1

Als letzten Schritt haben wir hier auch noch die nicht erwachsenen Probanden entfernt, wie in der Aufgabe gefordert. (Wenn dieser Schritt, erhält man dennoch dasselbe Ergebnis, da die DPQ-Scores für alle Minderjährigen NA sind, weil der DPQ-Frageboden bei diesen nicht erhoben wurde.)

Bepunktung

• Codes für Geschlecht und marital status korrekt aus Codebuch übertragen: 1 Punkt

• Codes für Geschlecht und marital status durch Begriffe ersetzt (oder ergänzt): 1 Punkt

• DPQ-Punktsummen angefügt: 1 Punkt

• sichergestellt, dass nur Probanden ab 18 in der Tabelle verbleiben: 1 Punkt

• Code zumindest kurz erklärt: ½ Punkt

• Kein Punktabzug, wenn relationship in left_join fehlt (nicht erforderlich)

• Kein Punktabzug, wenn unter-18-jährige nicht enfernt wurden, aber dafür bemerkt wurde, dass diese alle NA sind

Aufgabe 3

Wir berechnen die Mittelwerte der DPQ-Summen für Männer und Frauen:

full_tbl %>%
group_by( gender ) %>%
summarise( mean( dpq_score, na.rm=TRUE ) )

# A tibble: 2 × 2
  gender `mean(dpq_score, na.rm = TRUE)`
  <chr>                            <dbl>
1 female                            3.69
2 male                              2.77

Um Männer und Frauen zu vergleichen, führen wir einen t-Test durch:

t.test( dpq_score ~ gender, full_tbl )

    Welch Two Sample t-test

data:  dpq_score by gender
t = 7.8251, df = 5048.6, p-value = 6.133e-15
alternative hypothesis: true difference in means between group female and group male is not equal to 0
95 percent confidence interval:
 0.6939073 1.1578215
sample estimates:
mean in group female   mean in group male 
            3.693822             2.767958 

Der t-Test ist hoch signifikant (\(p<10^{-14}\)), der Unterschied ist also sicher echt.

Anmerkung: Da der t-Test ind er letzten Zeile auch die Mittelwerte ausgibt, ist die explizite Berechnung der Mittelwerte im Chunk davor eigentlich nicht erforderlich.

Um die Standardfehler der Mittelwerte pro Geschlecht zu erhalten, benutzen wir ebenfalls t.test, aber rufen die t-Test-Funktion un mit den Werten für jedes Geschlecht getrennt auf. Der Übersicht halber behalten wir nur die benötigten Spalten des Ergebnisses.

full_tbl %>%
group_by( gender ) %>%
summarise( broom::tidy( t.test( dpq_score ) ) ) %>%
select( gender, mean=estimate, conf.low, conf.high ) -> dpq_means

dpq_means

# A tibble: 2 × 4
  gender  mean conf.low conf.high
  <chr>  <dbl>    <dbl>     <dbl>
1 female  3.69     3.52      3.86
2 male    2.77     2.61      2.93

Alternativ kann man das 95%-KI auch so berechnen:

full_tbl %>%
filter( !is.na( dpq_score ) ) %>%
group_by( gender ) %>%
summarise( 
  n = n(),
  mean = mean( dpq_score ),
  sd = sd( dpq_score ) ) %>%
mutate(
  conf.low  = mean + qt( .025, n-1 ) * sd / sqrt(n),
  conf.high = mean + qt( .975, n-1 ) * sd / sqrt(n) )

# A tibble: 2 × 6
  gender     n  mean    sd conf.low conf.high
  <chr>  <int> <dbl> <dbl>    <dbl>     <dbl>
1 female  2590  3.69  4.43     3.52      3.86
2 male    2478  2.77  3.99     2.61      2.93

Da die Anzahl der Probanden sehr groß ist, ist es auch in Ordnung, qt( .025, n-1 ) einfach durch -2 und qt( .975, n-1 ) durch 2 zu ersetzen.

Nun die Säulendiagramme. Hier kann man geom_histogram verwenden.

ggplot( full_tbl ) +
geom_histogram( aes( x=dpq_score, y=after_stat(density) ), binwidth=1 ) +
facet_grid( rows=vars(gender) ) +
geom_vline( aes( xintercept=mean ), col="pink", data=dpq_means ) +
geom_vline( aes( xintercept=conf.low ), col="red", data=dpq_means ) +
geom_vline( aes( xintercept=conf.high ), col="red", data=dpq_means )  

Wir haben hier geom_histogram das Argument binwidth=1 gegeben, so dass jede Säule einer Einheit auf der Punkteskala entspricht. Das after_stat(density) ist nicht wirklich erforderlich, erleichtert aber den Vergleich zwischen den Geschlechtern.

Statt geom_histogram kann man auch geom_bar verwenden:

ggplot( full_tbl ) +
geom_bar( aes( x=dpq_score ) ) +
facet_grid( rows=vars(gender) )

Bepunktung

• Mittelwerte berechnet (ggf implizit in t-Test): 1 Punkt
• t-Test durchgeführt: 1 Punkt
• 95-KIs der Mittelwerte berechnet: 1 Punkt
• Histogramm geplottet, aufgeteilt nach Geschlecht: 1 Punkt
• Histogramm hat eine Säule pro Punktwert: 1 Punkt
• Mittelwert und KIs eingezeichnet: 1 Punkt
• Code kurz erklärt: ½ Punkt
• Ergebnis t-Test kommentiert: 1 Punkt
• Histogramm interpretiert: 1 Punkt

Aufgabe 4

Den Anteil von Probanden mit mindestens mittelschwerer Depression, aufgeschlüsselt nach Geschlecht und marital status, kann man wie folgt berechnen:

full_tbl %>%
filter( !is.na( dpq_score ) ) %>%
group_by( gender, marital_status ) %>%
summarise( fraction_depressed = mean( dpq_score >= 10 ) )

# A tibble: 14 × 3
# Groups:   gender [2]
   gender marital_status      fraction_depressed
   <chr>  <chr>                            <dbl>
 1 female divorced                        0.161 
 2 female living with partner             0.135 
 3 female married                         0.0717
 4 female never married                   0.108 
 5 female separated                       0.185 
 6 female widowed                         0.131 
 7 female <NA>                            0.118 
 8 male   divorced                        0.127 
 9 male   living with partner             0.0529
10 male   married                         0.0510
11 male   never married                   0.101 
12 male   separated                       0.0986
13 male   widowed                         0.151 
14 male   <NA>                            0.0661

Wenn man aber auch KIs braucht, geht man besser folgendermaßen vor:

full_tbl %>%
filter( !is.na( dpq_score ) ) %>%
group_by( gender, marital_status ) %>%
summarise( broom::tidy( binom.test( sum( dpq_score >= 10 ), n() ) ) ) %>%
select( gender, marital_status, fraction=estimate, conf.low, conf.high, n_depr=statistic, n=parameter ) -> depr_fracs

depr_fracs

# A tibble: 14 × 7
# Groups:   gender [2]
   gender marital_status      fraction conf.low conf.high n_depr     n
   <chr>  <chr>                  <dbl>    <dbl>     <dbl>  <dbl> <dbl>
 1 female divorced              0.161    0.122     0.206      50   311
 2 female living with partner   0.135    0.0926    0.186      30   223
 3 female married               0.0717   0.0573    0.0883     81  1130
 4 female never married         0.108    0.0805    0.142      46   424
 5 female separated             0.185    0.117     0.271      20   108
 6 female widowed               0.131    0.0930    0.178      35   267
 7 female <NA>                  0.118    0.0676    0.187      15   127
 8 male   divorced              0.127    0.0876    0.175      31   245
 9 male   living with partner   0.0529   0.0276    0.0905     12   227
10 male   married               0.0510   0.0396    0.0645     65  1275
11 male   never married         0.101    0.0746    0.133      45   446
12 male   separated             0.0986   0.0406    0.193       7    71
13 male   widowed               0.151    0.0848    0.240      14    93
14 male   <NA>                  0.0661   0.0290    0.126       8   121

HIer haben wir binom.test verwendet, aber nicht, um einen Binomial-Test durchzuführen, sondern nur, um ein binomiales KI zu erhalten. Dazu haben wir binom.test die Anzahl der Probanden mit DPQ-Score ab 10 (sum(depr_score>=10)) übergeben, sowie die Gesamtzahld er Probanden (n()) in der jeweiligen Gruppe.

Wir erzeugen damit Säulendiagramme:

depr_fracs %>%
ggplot() +
  geom_col( aes( x = marital_status, y = fraction, fill = gender ), position=position_dodge(.9), alpha=.6 ) +
  geom_errorbar( aes( x = marital_status, ymin = conf.low, ymax= conf.high, col = gender ), position=position_dodge(.9), width=.6 ) +
  ylab( "fraction with DPQ score ≥ 10")

Anmerkung: Oft wurde der marital status “refused” nicht (wie hier) durch NA ersetzt sondern beibehalten, was zu einer Säule mit extrem großem KI führt. Diese Säule sollte idealerweise entfernt werden, oder die y-Ache angepasst.

Bepunktung

• Anteil mit DPQ-Score >= 10 berechnet, aufgeschlüsselt nach Geschlecht und marital status: 1 Punkt
  • ½ Punkt Abzug, fals >10 statt >=10
• KIs berechnet: 2 Punkte
• Säulendiagramm erstellt: 1 Punkt
• Säulen für die beiden Geschlechter gut dargestellt: 1 Punkt
• Fehlerbalken eingetragen: 1 Punkt
• Code erklärt: ½ Punkt
• Plot interpretiert: 1 Punkt

Weitere Punkte

Punkte zur äußeren Form des gesamten Notebooks:

• Korrekt zu HTML formattiert mit allen Plots: 2 Punkte
• Für alle Ergebnis-Tabellen Tabellen-Anfang ausgegeben: 1 Punkt
• Notebook übersichtlich formattiert/gegliedert: 1 Punkt