Hausaufgabe vom 14.11.2024

Aufgabe 1

Berechnen Sie nochmals, wie schon in der vorigen Hausaufgabe, die durchschnittliche Körpergröße der erwachsenen Probanden aus der NHANES-Studie, aufgeschlüsselt nach Ethnie und Geschlecht.

suppressPackageStartupMessages(library(tidyverse))

read_csv("data_on_git/nhanes.csv") %>%  
  filter(!is.na(height), age >= 18) %>%
  group_by(gender, ethnicity) %>%
  summarise(Avg_Adult_Height = mean(height))
Rows: 8704 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): gender, ethnicity
dbl (4): subjectId, age, height, weight

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
`summarise()` has grouped output by 'gender'. You can override using the `.groups` argument.
# A tibble: 12 × 3
# Groups:   gender [2]
   gender ethnicity      Avg_Adult_Height
   <chr>  <chr>                     <dbl>
 1 female Mexican                    157.
 2 female NH Asian                   156.
 3 female NH Black                   162.
 4 female NH White                   161.
 5 female Other Hispanic             157.
 6 female Other/Mixed                162.
 7 male   Mexican                    170.
 8 male   NH Asian                   170.
 9 male   NH Black                   176.
10 male   NH White                   175.
11 male   Other Hispanic             169.
12 male   Other/Mixed                176.

Pivotieren Sie die Tabelle dann so, dass die Mittelwerte für die beiden Geschlechter nebeneinander statt untereinander stehen.

read_csv("data_on_git/nhanes.csv") %>%  
  filter(!is.na(height), age >= 18) %>%
  group_by(gender, ethnicity) %>%
  summarise(Avg_Adult_Height = mean(height)) %>%
  pivot_wider(id_cols="ethnicity", names_from="gender", values_from="Avg_Adult_Height")
Rows: 8704 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): gender, ethnicity
dbl (4): subjectId, age, height, weight

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
`summarise()` has grouped output by 'gender'. You can override using the `.groups` argument.
# A tibble: 6 × 3
  ethnicity      female  male
  <chr>           <dbl> <dbl>
1 Mexican          157.  170.
2 NH Asian         156.  170.
3 NH Black         162.  176.
4 NH White         161.  175.
5 Other Hispanic   157.  169.
6 Other/Mixed      162.  176.

Aufgabe 2

In einem Experiment soll der Einfluss verschiedener Wirkstoffe auf das Wachstum von Hefezellen gemessen werden. Dazu wird eine Suspension von Hefezellen in Wasser angesetzt. Die Suspension wird auf eine 48-Well-Platte verteilt; jedes Well erhält die gleiche Menge an Hefe. Anschließend werden die Wirkstoffe in die Wells gegeben. Jeweils 6 Wells, nämlich die, die in derselben “Spalte” der Platte liegen, erhalten denselben Wirkstoff. Die folgende Tabelle gibt an, welcher Wirkstoff in welche der 8 Spalten gegeben wurde:

wirkstoffe <- tibble(
  spalte = 1:8,
  wirkstoff = c( "Wasser", "DMSO", "Giftin", "HappyGrow", "KillEmAll", 
                 "BoringStuff", "Boostol", "Wasser." )
)

wirkstoffe
# A tibble: 8 × 2
  spalte wirkstoff  
   <int> <chr>      
1      1 Wasser     
2      2 DMSO       
3      3 Giftin     
4      4 HappyGrow  
5      5 KillEmAll  
6      6 BoringStuff
7      7 Boostol    
8      8 Wasser.

Anschließend wird in jedes Well eine für das Wachstum von Hefe geeignete Nährlösung hinzu gegeben, und die Platte dann für 3 Stunden bei 30 °C inkubiert.

Dann wird die Platte in einen. sog. Plattenscannner gegeben. Ein solche Scanner beleuchtet jedes Well von unten mit Licht einer vorgegebenen Wellenlänge und misst von oben, wie viel Licht durch das Well hindurch dringen kann. Je mehr Hefezellen im Well sind, desto trüber wird die Suspension, und desto größer ist die gemessene Licht-Absorption (optische Dichte, “OD”).

Der Plattenscanner erzeugt diese Excel-Datei als Ergebnis: example_plate_scan.xlsx. Laden Sie die Datei mit Excel (oder einem anderen Tabellenkalkulationsprogramm) und sehen Sie sie sich an.

Einlesen der Datei

Sie können diese Datei in eine CSV-Datei umwandeln und dann in R einlesen, oder Sie verwenden die Funktion read_excel aus dem Paket readxl, die Excel-Dateien direkt einlesen kann. In beiden Fällen müssen Sie die Zeilen über der eigentlichen Datentabelle entfernen. Bei read_xl können Sie hierzu das optionale Argument skip verwenden. So liest read_excel( "file.xls", skip=7) die Tabelle auf dem ersten Arbeitsblatt (Sheet) der Excel-Datei file.xls und überspringt dabei die ersten 7 Zeilen.

Die eingelesene Tabelle sollte 9 Spalten haben: 8 für die Spalten der Mikrotiterplatte und davor eine Spalte mit den Buchstaben, die die Zeilen nummerieren. Geben Sie dieser Spalte einen geeigneten Namen, z.B. zeile. Hierzu können Sie das Tidyverse-Verb rename verwenden.

library(readxl)

read_excel("data_on_git/example_plate_scan.xlsx", skip=5) %>%
  rename("row"="...1") -> plate_scan
New names:
• `` -> `...1`
plate_scan
# A tibble: 6 × 9
  row     `1`   `2`   `3`   `4`   `5`   `6`   `7`   `8`
  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 A     0.916 0.856 0.17   2.58 0.033 0.872  1.38 0.587
2 B     1.11  0.813 0.207  2.78 0.054 1.12   1.76 0.929
3 C     0.993 1.06  0.155  4.34 0.043 1.30   2.19 1.02 
4 D     0.926 1.11  0.221  3.19 0.054 1.36   1.70 1.16 
5 E     1.04  0.992 0.256  2.51 0.048 0.95   2.06 0.971
6 F     0.7   0.792 0.176  2.46 0.041 0.901  1.41 0.963

Pivotieren

Pivotieren Sie die Tabelle in ein “tidy-table”-Format, d.h., sorgen Sie dafür, das jedes Well eine eigene Tabellen-Zeile hat. Sie sollten nun 3 Spalten haben, für Platten-Zeile, Platten-Spalte, und Messwert.

plate_scan %>%
  pivot_longer(cols= -row, names_to="column", values_to="OD") ->plate_scan_tidy

plate_scan_tidy
# A tibble: 48 × 3
   row   column    OD
   <chr> <chr>  <dbl>
 1 A     1      0.916
 2 A     2      0.856
 3 A     3      0.17 
 4 A     4      2.58 
 5 A     5      0.033
 6 A     6      0.872
 7 A     7      1.38 
 8 A     8      0.587
 9 B     1      1.11 
10 B     2      0.813
# ℹ 38 more rows

Streudiagramm

Erzeugen Sie nun ein Streudiagramm, dass die gemessen optische Dichte für jeden Wirkstoff darstellt. Es könnte in etwa so aussehen:

plate_scan_tidy %>%
  transform(column=as.numeric(column)) %>%
  left_join(wirkstoffe, by=c("column"="spalte")) %>%
  mutate(Wirkstoff = fct_inorder(wirkstoff) ) %>%
  ggplot +
  geom_point(aes(x=Wirkstoff, y=OD)) +
  scale_y_log10()

library("ggbeeswarm")

plate_scan_tidy %>%
  transform(column=as.numeric(column)) %>%
  left_join(wirkstoffe, by=c("column"="spalte")) %>%
  mutate(Wirkstoff = fct_inorder(wirkstoff)) %>%
  ggplot +
  geom_beeswarm(aes(x=Wirkstoff, y=OD)) +
  scale_y_log10()

Mittelwerte

Ermitteln Sie nun noch die Mittelwerte der OD-Werte für jeden Wirkstoff.

plate_scan_tidy %>%
  transform(column = as.numeric(column)) %>%
  left_join(wirkstoffe, by=c("column"="spalte")) %>%
  mutate(Wirkstoff = fct_inorder(wirkstoff)) %>%
  select(-column, -wirkstoff) %>%
  group_by(Wirkstoff) %>%
  summarise(Avg_OD = mean(OD)) -> plate_scan_avg

plate_scan_avg
# A tibble: 8 × 2
  Wirkstoff   Avg_OD
  <fct>        <dbl>
1 Wasser      0.947 
2 DMSO        0.937 
3 Giftin      0.198 
4 HappyGrow   2.98  
5 KillEmAll   0.0455
6 BoringStuff 1.08  
7 Boostol     1.75  
8 Wasser.     0.938

Nützlich wäre auch, diese Mittelwerte durch den Mittelwert für Wasser zu teilen, um leichter zu sehen, welche Wirkstoffe das Wachstum fördern und welche es hemmen.

plate_scan_avg %>%
  mutate(WasserAvgOD = ((Avg_OD[1] + Avg_OD[8])/2)) %>%
  group_by(Wirkstoff) %>%
  summarise(Avg_OD_corr = Avg_OD/WasserAvgOD)
# A tibble: 8 × 2
  Wirkstoff   Avg_OD_corr
  <fct>             <dbl>
1 Wasser           1.00  
2 DMSO             0.993 
3 Giftin           0.209 
4 HappyGrow        3.16  
5 KillEmAll        0.0483
6 BoringStuff      1.15  
7 Boostol          1.86  
8 Wasser.          0.995

Platten-Plot

In diesem Plot habe ich einige Feinheiten geändert, so dass ihr Plot anders aussehen könnte. Ich habe die Punkte extra groß gemacht, damit sie wie Wells aussehen, und ich habe die y-Achse umgekehrt (weil Zeilen in Mikrotiterplatten von oben nach unten numeriert werden, ggplot die y-Achse aber von unten nach oben ansteigen lässt), die Farbskala logarithmisch gemacht, und die Palette geändert.

Einen solchen Plot benutzt man gerne, um sicher zu stellen, dass es keine Pippetierfehler oder Randeffekte gab.

plate_scan_tidy %>%
  transform(column= as.numeric(column)) %>%
  left_join(wirkstoffe, by=c("column"="spalte")) %>%
  mutate(Wirkstoff = fct_inorder(wirkstoff)) %>%
  ggplot + 
  geom_point(aes(x=Wirkstoff, y=row, col=OD), size=15) +
  scale_y_discrete(limits = rev(unique(sort(plate_scan_tidy$row))))+
  scale_colour_gradient(name ="OD", trans="log10",low="beige", high="brown")

Aufgabe 3: Stichproben-Mittelwert

Mit dem folgenden Code reduzieren wir die NHANES-Tabelle zunächste auf nur die erwachsenen Männer, für die die Körpergröße vorliegt.

suppressPackageStartupMessages(library(tidyverse))

read_csv("data_on_git/nhanes.csv") -> nhanes
Rows: 8704 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): gender, ethnicity
dbl (4): subjectId, age, height, weight

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
nhanes %>%
  filter(!is.na(gender), !is.na(age), !is.na(height), gender =="male", age >= 18) -> nhanes_men_only

nhanes_men_only
# A tibble: 2,630 × 6
   subjectId gender   age height weight ethnicity  
       <dbl> <chr>  <dbl>  <dbl>  <dbl> <chr>      
 1     93706 male      18   176.   66.3 NH Asian   
 2     93711 male      56   171.   62.1 NH Asian   
 3     93712 male      18   173.   58.9 Mexican    
 4     93713 male      67   179.   74.9 NH White   
 5     93715 male      71   171.   65.6 Other/Mixed
 6     93716 male      61   159.   77.7 NH Asian   
 7     93717 male      22   174.   74.4 NH White   
 8     93718 male      45   157.   54.4 NH Black   
 9     93723 male      64   170.   64.9 NH White   
10     93727 male      70   162.   62.7 NH Asian   
# ℹ 2,620 more rows

Wir bestimmen den Mittelwert der Körpergröße aller dieser Männer

nhanes_men_only %>%
  summarise(avg_man_height = mean(height)) -> population_mean

population_mean
# A tibble: 1 × 1
  avg_man_height
           <dbl>
1           173.

Dann wählen wir mit sample_n 10 Zeilen zufällig aus. Jedesmal, wenn Sie den Code ausführen, werden andere Zeilen ausgewählt.Bestimmen Sie mit summarise und mean den Mittelwert dieser 10 Männer.

nhanes_men_only %>%
  sample_n(10) %>%
  summarise(sample_mean = mean(height))
# A tibble: 1 × 1
  sample_mean
        <dbl>
1        172.
nhanes_men_only %>%
  sample_n(10) %>%
  summarise(sample_mean = mean(height))
# A tibble: 1 × 1
  sample_mean
        <dbl>
1        174.
nhanes_men_only %>%
  sample_n(10) %>%
  summarise(sample_mean = mean(height))
# A tibble: 1 × 1
  sample_mean
        <dbl>
1        170.

Mit map_dfr können wir diesen Code mehrmals (hier: 15 mal) ausführen lassen, und die 1x1-Tabellen zu einer langen Tabelle zusammen fügen lassen

map_dfr(1:15, ~{
  nhanes_men_only %>%
    sample_n(10) %>%
    summarise(sample_mean=mean(height))
})
# A tibble: 15 × 1
   sample_mean
         <dbl>
 1        173.
 2        172.
 3        176.
 4        171.
 5        174.
 6        174.
 7        176.
 8        176.
 9        177.
10        175.
11        176.
12        173.
13        170.
14        173.
15        176.

Notieren Sie jeweils, wie stark der Mittelwert über die 10 Männer jeweils abweicht vom Mittelwert aus der gesamten Tabelle mit allen Männern.

map_dfr(1:15, ~{
  nhanes_men_only %>%
    sample_n(10) %>%
    summarise(sample_mean=mean(height))
}) %>% 
  summarise( sd = sum((sample_mean - population_mean$avg_man_height)^2)/15)
# A tibble: 1 × 1
     sd
  <dbl>
1  7.05

Wiederholen Sie dies nun nochmals, aber ziehen Sie diesmal jeweils 100 statt nur 10 zufällig ausgewählte Zeilen.

map_dfr(1:15, ~{
  nhanes_men_only %>%
  sample_n(100) %>%
  summarise(sample_mean = mean(height))
}
)
# A tibble: 15 × 1
   sample_mean
         <dbl>
 1        175.
 2        173.
 3        175.
 4        174.
 5        174.
 6        173.
 7        175.
 8        174.
 9        174.
10        174.
11        176.
12        174.
13        174.
14        172.
15        173.

Notieren Sie wieder die Abweichungen der Stichproben-Mittelwerte vom gesamt-Mittelwert.

map_dfr(1:15, ~{
  nhanes_men_only %>%
  sample_n(100) %>%
  summarise(sample_mean = mean(height))
}) %>%
summarise(sd = sum((sample_mean - population_mean$avg_man_height)^2)/100)
# A tibble: 1 × 1
     sd
  <dbl>
1 0.122