Tag 4 Statistische Tests und Power

4.2 Der zentrale Grenzwertsatz (CLT)

Von Verteilungen zu Verteilungen von Mittelwerten

Was wäre, wenn Johanna als Master Studentin nicht ganz so fleißig gewesen wäre, und statt 300 je Zell Typ nur je 30 mal gezählt hätte?

Der Einfachheit nehmen wir an, dass die 300 Zellen komplett repräsentativ für alle WT bzw. KO Zellen sind, also die Gesamtpopulation darstellen.

  • Nun ziehen wir zufällig daraus je 30 Zellen.
  • Das machen wir nun 1000 mal.
  • Wie sieht nun die Verteilung der Mittelwerte aus?

4.2.1 Daten einlesen

Zunächst laden wir das Tidyverse und unsere Daten:

# load tidyverse
library(tidyverse)

# import data
data <- read_csv("data/03_inclusion_bodies.csv")

# attach data, handle with CARE!
attach(data)

attach macht die Spalten des Datensets (in unserem Fall ko and wt) global verfügbar. Daher können wir im Folgenden wt statt data$wt oder data["wt"] schreiben.

Verwende attach mit Vorsicht! Denn die Reihenfolge der beiden Vektoren ist nun unabhängig voneinander. Operationen, die Du mit den Vektoren wt und ko ausführst, beeinflussen nicht den jeweils anderen Vektor oder den zugrundeliegenden Data.frame data.

4.2.2 Der einfache Fall

Für den einfachen Fall ziehen wir 30 Punkte aus einem der beiden Vektoren und berechnen den Mittelwert.

# sample n = 30 points from a vector
subset <- sample(wt, 30)

# calculate mean
mean(subset)
## [1] 2.133333

4.2.3 Abstraktion für viele Fälle

Da wir diese Zeilen nicht 1000 mal tippen wollen schreiben wir eine Funktion, die genau das tut. Sie nimmt einen Vektor (x) und zieht n Punkte daraus.

get_subset_mean <- function(n, x) {
  subset <- sample(x, n)
  mean(subset)
}

Die Funktion funktioniert:

get_subset_mean(2, wt)
## [1] 6

Mit der map Funktion aus dem purrr-Package rufen wir diese Funktion nun 1000 mal auf uns speichern das Ergebnis in einem Vektor, bestehend aus Kommazahlen (daher map_dbl für map double).

N <- 1000
Ms <- map_dbl(1:N, ~ get_subset_mean(30, wt) )

4.2.4 Exkurs: Lambda-Funktionen

Das ~-Symbol (Tilde) erstellt eine sogenannte lambda Funktion, also eine Funktion, der wir keinen Namen geben, da wir sie nur einmal brauchen. Im Kontext der Familie von map Funktionen wird ~ übersetzt:

add_one <- function(x) x + 1 Ist equivalent zu ~ .x + 1. Hierbei heißt das Argument der Funktion immer .x.

Die lambda Funktion ~ get_subset_mean(30, wt) nimmt nacheinander die Zahlen in 1:N als Eingabe, nennt sie .x, aber verwendet das Argument .x überhaubt nicht. In unserem Fall bedeutet das also, dass die Funktion get_subset_mean bloß 1000 mal aufgerufen wird mit den Argumenten 30 und wt.

4.2.5 Populationsverteilung vs. Verteilung der Mittelwerte

Mit großem n wird die Verteilung der Mittelwerte symmetrisch, obwohl die Populationsverteilung unsymmetrisch ist.

hist(Ms)

hist(wt)

Nach dem Central Limit Theorem folgen die Mittelwerte einer Normalverteilung. Dies lässt sich zeigen:

# show qqnorm of Ms
qqnorm(Ms)
qqline(Ms, col = "red")

4.2.6 Exkurs: Quantiles

Am Beispiel einer Normalverteilung zeigt dnorm die Wahrscheinlichkeitsdichteverteilung (Im Englischen Probability Density Function, PDF).

curve(dnorm, -3, 3)
Die Wahrscheinlichkeit, einen Wert zwischen z.B. -1 und 1 zu ziehen, entspricht dem Integral der Funktion in diesem Bereich.

Abb. 4.1: Die Wahrscheinlichkeit, einen Wert zwischen z.B. -1 und 1 zu ziehen, entspricht dem Integral der Funktion in diesem Bereich.

Wenn wir nun das Integral von \(-\infty\) bis \(x\) berechnen erhalten wir die kumulative Wahrscheinlichkeitsfunktion pnorm.

curve(pnorm, -3, 3)
pnorm von x ist also die Wahrscheinlichkeit, einen Wert kleiner oder gleich x aus dnorm zu ziehen.

Abb. 4.2: pnorm von x ist also die Wahrscheinlichkeit, einen Wert kleiner oder gleich x aus dnorm zu ziehen.

qnomr, die Quantil-Funktion, ist nun die Inverse dieser Funkion, hat also die Achsen vertauscht.

curve(qnorm, 0, 1)

4.3 Die T-Verteilung

Der zentrale Grenzwertsatz gilt erst bei großen sample sizes. Darunter ist die Verteilung der Mittelwerte an den Rändern dicker, als eine Normalverteilung.

Daher verwenden wir (sicherheitshalber) für die meisten Daten und statistischen Tests die T-Verteilung, nicht die Normalverteilung.

Für 3 Freiheitsgrade:

curve(dnorm, -3, 3)
tdist <- function(x) dt(x, df = 3)
curve(tdist, -3, 3, add = TRUE, col = "red")
T-Verteilung in Rot

Abb. 4.3: T-Verteilung in Rot

Für 30 Freiheitsgrade

curve(dnorm, -3, 3)
tdist <- function(x) dt(x, df = 30)
curve(tdist, -3, 3, add = TRUE, col = "red")

4.3.1 Konfidenz Intervalle (CI)

Gleichermaßen können wir mit der t-Verteilung 95% Konfidenzintervalle berechnen.

Das 95% Konfidenzintervall eines Sample-Mittelwerts umschließt den wahren Mittelwert (den Mittelwert der Gesamtpopulation) in 95% der Fälle (wenn Du das Experiment unendlich oft wiederholen würdest).

# show hist of Ms
hist(Ms)

# add line for true mean of wt
abline(v = mean(wt), col = "red", lwd = 3)
Wahrer Mittelwert in Rot

Abb. 4.4: Wahrer Mittelwert in Rot 

# use t.test for conf.int
subset <- sample(wt, 30)
mean(subset)
## [1] 3.166667
lims <- t.test(subset)$conf.int

hist(subset, breaks = 20)
abline(v = lims[1], col = "red")
abline(v = lims[2], col = "red")

4.4 Signifikanztests

Signifikanztests beantworten die Frage:

Unter der Annahme, dass es zwischen zwei (oder mehr) Gruppen keinen Unterschied gibt (also sie aus der gleichen Verteilung stammen), wie wahrscheinlich ist es dann, einen so großen, oder größeren, Unterschied, wie beobachtet, zu finden?

4.4.1 Students T-Test

Für normalverteilte Daten verwenden wir den Students T-Test.

t.test(wt, ko)
## 
## 	Welch Two Sample t-test
## 
## data:  wt and ko
## t = -26.929, df = 476.21, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.19393 -11.39940
## sample estimates:
## mean of x mean of y 
##  3.043333 15.340000
subset_wt <- sample(wt, 3)
subset_ko <- sample(ko, 3)
t.test(subset_wt, subset_ko)
## 
## 	Welch Two Sample t-test
## 
## data:  subset_wt and subset_ko
## t = -6.8224, df = 2.6162, p-value = 0.009833
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -16.081856  -5.251477
## sample estimates:
## mean of x mean of y 
##   1.00000  11.66667

4.4.2 Wilcoxon Rank Sum Test

Sind unsere Daten nicht normalverteilt, wie es hier der Fall ist, ist ein sogenannter nicht-parametrischer Test angebracht. Diese Tests machen nicht die Annahme der Normalität.

Der Wilcoxon Rank Sum Test (Auch Mann-Whitney U Test genannt) wandelt die eigentlichen Werte der Datenpunkte zunächst in Ränge um. Also:

“Der wievieltniedrigste Punkt ist es?”

Beispiel:

x <- c(1, 2, 20, 1239132, 3, 5)
rank(x)
## [1] 1 2 5 6 3 4

Und an unseren Daten:

wilcox.test(wt, ko)
## 
## 	Wilcoxon rank sum test with continuity
## 	correction
## 
## data:  wt and ko
## W = 4707.5, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

4.5 Type I und Type II Errors

  • Type I: False Positives (rejection of a true null hypothesis)
  • Type II: False Negatives (non-rejection of false null hypothesis)

4.5.1 Type I: False Positives

Ein Typ I Fehler würde bedeuten, dass wir sagen, ein Unterschied zwischen den Gruppen existiere, obwohl alle Werte aus der gleichen Verteilung stammen und daher kein Unterschied besteht.

Bei einem P-Value von \(\leq\) 0.05 (= 5 %) wird typischerweise gesagt, dass ein “statistisch signifikanter” Unterschied besteht. Mit diesem typischen Cutoff von 0.05 akzeptieren wir also, allein durch unsere Definition, bereits (mindestens) 5 % falsch Positive.

Dieser Cutoff, das Signifikanzlevel, wird auch \(\alpha\) genannt.

\[\alpha=\text{Type I error rate}\]

curve(dnorm, -3 , 3)

# define n
n <- 3

# draw twice from the same normal distributions
draw1 <- rnorm(n)
draw2 <- rnorm(n)

# do the t.test
t.test(draw1, draw2)
## 
## 	Welch Two Sample t-test
## 
## data:  draw1 and draw2
## t = -0.31189, df = 3.0148, p-value = 0.7754
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.651341  2.997886
## sample estimates:
##  mean of x  mean of y 
## -0.1263381  0.2003894
# do the rank sum test
wilcox.test(draw1,  draw2)
## 
## 	Wilcoxon rank sum exact test
## 
## data:  draw1 and draw2
## W = 4, p-value = 1
## alternative hypothesis: true location shift is not equal to 0

Teste diesen Code mit unterschiedlichen Werten für n.

4.5.2 Type II errors

Wenn zwischen zwei Gruppen ein tatsächlicher Unterschied besteht, wir aber bei unserem Test ein nicht signifikantes Ergebnis erhalten, begehen wir einen Fehler vom Typ II.

Ein Typ II Fehler ist also ein falsch negatives Ergebnis.

\[\beta=\text{Type II error rate}\]

# define n
n <- 15

# draw from different normal distributions
draw1 <- rnorm(n, 0, 1)
draw2 <- rnorm(n, 1, 1)

# convert to tibble, pivot_longer, plot points
tibble(draw1, draw2) %>% 
  pivot_longer(c(1,2)) %>% 
  ggplot(aes(name, value, group = 1)) +
  geom_jitter(width = 0.2) +
  stat_summary(geom = "line", fun.y = mean, lwd = 1, lty = 2) +
  stat_summary(geom = "point", fun.y = mean, color = "red", size = 2)

# do the t.test
t.test(draw1, draw2)
## 
## 	Welch Two Sample t-test
## 
## data:  draw1 and draw2
## t = -2.5888, df = 27.319, p-value = 0.01525
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.9133144 -0.2219188
## sample estimates:
##  mean of x  mean of y 
## 0.06255136 1.13016800
# do the rank sum test
wilcox.test(draw1,  draw2)
## 
## 	Wilcoxon rank sum exact test
## 
## data:  draw1 and draw2
## W = 57, p-value = 0.0209
## alternative hypothesis: true location shift is not equal to 0

4.5.3 Statistical Power

Als statistische Power bezeichnen wir die Wahrscheinlichkeit eines Tests, einen wahren Unterschied zwischen Gruppen bei dem gewählten \(\alpha\) auch tatsächlich als statistisch signifikant zu kennzeichnen und damit ein wahr positives Ergebnis zu produzieren.

\[Power = 1-\beta\]

Für den T-Test können wir in R die folgende Funktion verwenden, die uns hier beispielsweise die Frage beantwortet: “Wie viele Proben muss ich pro Gruppe (mindestens) nehmen, um einen erwarteten Unterschied der Mittelwerte von 1 mit einer Standardabweichung von 1 in 80% der Fälle auch tatsächlich als solchen zu erkennen?”

# show power.t.test
power.t.test(delta = 1, sd = 1, power = 0.8)
## 
##      Two-sample t test power calculation 
## 
##               n = 16.71477
##           delta = 1
##              sd = 1
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

4.6 Probleme

4.6.1 The Jelly Bean Problem (Multiple Testing)

Die False Discovery Rate (FDA) kontrollieren:

  • Bonferroni Korrektur
# show p.adjust
p.adjust(c(0.05, 0.0001, 0.003323, 0.7), method = "bonferroni")
## [1] 0.200000 0.000400 0.013292 1.000000

Jeder P-Value wird mutipliziert mit der Zahl der Tests

  • Benjamini-Hochberg Prozedur
    • Ordne alle p-values in aufsteigender Reihenfolge
    • Wähle eine FDR (\(q\)) und nenne die Anzahl deiner Tests \(m\)
    • Finde den größten p-value für den gilt: \(p \leq iq/m\) mit dem Index des p-values \(i\).
# show p.adjust for BH
p.adjust(c(0.05, 0.0001, 0.003323, 0.7), method = "BH")
## [1] 0.06666667 0.00040000 0.00664600 0.70000000

4.6.2 The Base Rate Fallacy

Beispiel: Mammogramm

  • Sensitivity = Power = true positive rate
  • Specificity = true negative rate = \(1-\alpha\)
# calculate 
total <- 1000
positives <- 10
negatives <- total - positives
sensitivity <- 0.9
specificity <- 1 - 0.08
true_positives  <- sensitivity * positives
false_positives <- (1 - specificity) * negatives
p_positive <- true_positives / (true_positives + false_positives)
p_positive
## [1] 0.1020408

4.7 Übungen

4.7.1 Zeige die folgenden Statements anhand von Simulationen in R

4.7.1.1 Wenn Du die Sample Size n erhöhst, verringert sich der Standard Error of the Mean mit der Wurzel von n.

  • Ziehe 10 Zellen aus dem wt (oder ko) Vektor
  • Berechne den Mittelwert
  • Wiederhole das Ganze 1000 mal
  • Berechne die SD der 1000 Mittelwerte
  • Ziehe nun 40 statt 10 Zellen und wiederhole die Schritte
  • Schreibe eine Funktion, die n Zellen zieht, die Schritte ausführt und die SD ausgibt
  • Füttere die Zahlen von 1 bis 100 in die Funktion (mit map_dbl) und plotte das Ergebnis

4.7.1.2 Lösung

library(tidyverse)

read_csv("data/03_inclusion_bodies.csv") %>% attach()
draw <- sample(wt, 10)
mean(draw)
## [1] 3.5
get_sample_mean <- function(x, n) {
  draw <- sample(x, n, replace = TRUE)
  mean(draw)
}
get_sample_mean(wt, 10)
## [1] 3.9
N <- 1000
many_means <- map_dbl(1:N, ~ get_sample_mean(wt, 10) )
sd(many_means)
## [1] 1.245121
many_means <- map_dbl(1:N, ~get_sample_mean(wt, 40))
sd(many_means)
## [1] 0.6121342
# explain default arguments and scoping
get_sd_of_many_means <- function(n, x, N = 1000) {
  many_means <- map_dbl(1:N, ~get_sample_mean(x, n))
  sd(many_means)
}
n_max <- 100
sds_by_n <- map_dbl(1:n_max, get_sd_of_many_means, x = wt)
plot(x = 1:n_max,
     y = sds_by_n, type = "p")

curve(sd(wt)/sqrt(x), add = TRUE, col = "red")

# talk about difference between sample sd and population sd

4.7.1.3 Ein 95% Konfidenzintervall eines Samples enthält den Populationsmittelwert mit einer Wahrscheinlichkeit von 95%.

  • Ziehe 30 Zellen aus dem wt (oder ko) Vektor
  • Berechne die Limits des CI (Confidence Interval)
  • Schreibe eine Funktion, oben genanntes tut und führe sie 1000 mal aus.
  • Schreibe eine Funktion, die testet, ob ein Set an Limits den wahren Mittelwert einschließt
  • Wende sie auf die 1000 Sets der Limits an
  • Wie oft (prozentual) erhältst du TRUE?

4.7.1.4 Lösung

draw <- sample(wt, 30)
test_results <- t.test(draw)
test_results
## 
## 	One Sample t-test
## 
## data:  draw
## t = 3.9747, df = 29, p-value = 0.0004285
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.165051 3.634949
## sample estimates:
## mean of x 
##       2.4
summary(test_results)
##             Length Class  Mode     
## statistic   1      -none- numeric  
## parameter   1      -none- numeric  
## p.value     1      -none- numeric  
## conf.int    2      -none- numeric  
## estimate    1      -none- numeric  
## null.value  1      -none- numeric  
## stderr      1      -none- numeric  
## alternative 1      -none- character
## method      1      -none- character
## data.name   1      -none- character
str(test_results)
## List of 10
##  $ statistic  : Named num 3.97
##   ..- attr(*, "names")= chr "t"
##  $ parameter  : Named num 29
##   ..- attr(*, "names")= chr "df"
##  $ p.value    : num 0.000428
##  $ conf.int   : num [1:2] 1.17 3.63
##   ..- attr(*, "conf.level")= num 0.95
##  $ estimate   : Named num 2.4
##   ..- attr(*, "names")= chr "mean of x"
##  $ null.value : Named num 0
##   ..- attr(*, "names")= chr "mean"
##  $ stderr     : num 0.604
##  $ alternative: chr "two.sided"
##  $ method     : chr "One Sample t-test"
##  $ data.name  : chr "draw"
##  - attr(*, "class")= chr "htest"
test_results$conf.int
## [1] 1.165051 3.634949
## attr(,"conf.level")
## [1] 0.95
get_sample_ci <- function(x, n) {
  draw <- sample(x, n, replace = TRUE)
  t.test(draw)$conf.int
}
get_sample_ci(wt, 30)
## [1] 1.272275 3.727725
## attr(,"conf.level")
## [1] 0.95
CIs <- map(1:1000, ~get_sample_ci(wt, 30))
head(CIs)
## [[1]]
## [1] 1.319300 4.147366
## attr(,"conf.level")
## [1] 0.95
## 
## [[2]]
## [1] 1.55864 5.24136
## attr(,"conf.level")
## [1] 0.95
## 
## [[3]]
## [1] 1.564104 4.502562
## attr(,"conf.level")
## [1] 0.95
## 
## [[4]]
## [1] 1.402216 3.997784
## attr(,"conf.level")
## [1] 0.95
## 
## [[5]]
## [1] 1.211124 2.588876
## attr(,"conf.level")
## [1] 0.95
## 
## [[6]]
## [1] 2.414082 6.319252
## attr(,"conf.level")
## [1] 0.95
# explain list subsetting
CIs[1]
## [[1]]
## [1] 1.319300 4.147366
## attr(,"conf.level")
## [1] 0.95
CIs[[1]]
## [1] 1.319300 4.147366
## attr(,"conf.level")
## [1] 0.95
CIs[[1]][1]
## [1] 1.3193
test_ci <- function(limits, true_mean) {
  limits[1] < true_mean & limits[2] > true_mean
}
# sidenote
# note difference between & and &&
c(TRUE, TRUE) &  c(FALSE, TRUE)
## [1] FALSE  TRUE
c(TRUE, TRUE) && c(FALSE, TRUE)
## [1] FALSE
results <- map_lgl(CIs, test_ci, true_mean = mean(wt))
head(results)
## [1]  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
mean(results)
## [1] 0.936
# talk about organisation of code
# source("test.R")

4.7.1.5 Weitere Hinweise

# check out
?map
# for the subtle differences between:
map(1:100, ~ .x + 1)
## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 4
## 
## [[4]]
## [1] 5
## 
## [[5]]
## [1] 6
## 
## [[6]]
## [1] 7
## 
## [[7]]
## [1] 8
## 
## [[8]]
## [1] 9
## 
## [[9]]
## [1] 10
## 
## [[10]]
## [1] 11
## 
## [[11]]
## [1] 12
## 
## [[12]]
## [1] 13
## 
## [[13]]
## [1] 14
## 
## [[14]]
## [1] 15
## 
## [[15]]
## [1] 16
## 
## [[16]]
## [1] 17
## 
## [[17]]
## [1] 18
## 
## [[18]]
## [1] 19
## 
## [[19]]
## [1] 20
## 
## [[20]]
## [1] 21
## 
## [[21]]
## [1] 22
## 
## [[22]]
## [1] 23
## 
## [[23]]
## [1] 24
## 
## [[24]]
## [1] 25
## 
## [[25]]
## [1] 26
## 
## [[26]]
## [1] 27
## 
## [[27]]
## [1] 28
## 
## [[28]]
## [1] 29
## 
## [[29]]
## [1] 30
## 
## [[30]]
## [1] 31
## 
## [[31]]
## [1] 32
## 
## [[32]]
## [1] 33
## 
## [[33]]
## [1] 34
## 
## [[34]]
## [1] 35
## 
## [[35]]
## [1] 36
## 
## [[36]]
## [1] 37
## 
## [[37]]
## [1] 38
## 
## [[38]]
## [1] 39
## 
## [[39]]
## [1] 40
## 
## [[40]]
## [1] 41
## 
## [[41]]
## [1] 42
## 
## [[42]]
## [1] 43
## 
## [[43]]
## [1] 44
## 
## [[44]]
## [1] 45
## 
## [[45]]
## [1] 46
## 
## [[46]]
## [1] 47
## 
## [[47]]
## [1] 48
## 
## [[48]]
## [1] 49
## 
## [[49]]
## [1] 50
## 
## [[50]]
## [1] 51
## 
## [[51]]
## [1] 52
## 
## [[52]]
## [1] 53
## 
## [[53]]
## [1] 54
## 
## [[54]]
## [1] 55
## 
## [[55]]
## [1] 56
## 
## [[56]]
## [1] 57
## 
## [[57]]
## [1] 58
## 
## [[58]]
## [1] 59
## 
## [[59]]
## [1] 60
## 
## [[60]]
## [1] 61
## 
## [[61]]
## [1] 62
## 
## [[62]]
## [1] 63
## 
## [[63]]
## [1] 64
## 
## [[64]]
## [1] 65
## 
## [[65]]
## [1] 66
## 
## [[66]]
## [1] 67
## 
## [[67]]
## [1] 68
## 
## [[68]]
## [1] 69
## 
## [[69]]
## [1] 70
## 
## [[70]]
## [1] 71
## 
## [[71]]
## [1] 72
## 
## [[72]]
## [1] 73
## 
## [[73]]
## [1] 74
## 
## [[74]]
## [1] 75
## 
## [[75]]
## [1] 76
## 
## [[76]]
## [1] 77
## 
## [[77]]
## [1] 78
## 
## [[78]]
## [1] 79
## 
## [[79]]
## [1] 80
## 
## [[80]]
## [1] 81
## 
## [[81]]
## [1] 82
## 
## [[82]]
## [1] 83
## 
## [[83]]
## [1] 84
## 
## [[84]]
## [1] 85
## 
## [[85]]
## [1] 86
## 
## [[86]]
## [1] 87
## 
## [[87]]
## [1] 88
## 
## [[88]]
## [1] 89
## 
## [[89]]
## [1] 90
## 
## [[90]]
## [1] 91
## 
## [[91]]
## [1] 92
## 
## [[92]]
## [1] 93
## 
## [[93]]
## [1] 94
## 
## [[94]]
## [1] 95
## 
## [[95]]
## [1] 96
## 
## [[96]]
## [1] 97
## 
## [[97]]
## [1] 98
## 
## [[98]]
## [1] 99
## 
## [[99]]
## [1] 100
## 
## [[100]]
## [1] 101
# and
add_one <- function(x) x + 1
map(1:100, add_one)
## [[1]]
## [1] 2
## 
## [[2]]
## [1] 3
## 
## [[3]]
## [1] 4
## 
## [[4]]
## [1] 5
## 
## [[5]]
## [1] 6
## 
## [[6]]
## [1] 7
## 
## [[7]]
## [1] 8
## 
## [[8]]
## [1] 9
## 
## [[9]]
## [1] 10
## 
## [[10]]
## [1] 11
## 
## [[11]]
## [1] 12
## 
## [[12]]
## [1] 13
## 
## [[13]]
## [1] 14
## 
## [[14]]
## [1] 15
## 
## [[15]]
## [1] 16
## 
## [[16]]
## [1] 17
## 
## [[17]]
## [1] 18
## 
## [[18]]
## [1] 19
## 
## [[19]]
## [1] 20
## 
## [[20]]
## [1] 21
## 
## [[21]]
## [1] 22
## 
## [[22]]
## [1] 23
## 
## [[23]]
## [1] 24
## 
## [[24]]
## [1] 25
## 
## [[25]]
## [1] 26
## 
## [[26]]
## [1] 27
## 
## [[27]]
## [1] 28
## 
## [[28]]
## [1] 29
## 
## [[29]]
## [1] 30
## 
## [[30]]
## [1] 31
## 
## [[31]]
## [1] 32
## 
## [[32]]
## [1] 33
## 
## [[33]]
## [1] 34
## 
## [[34]]
## [1] 35
## 
## [[35]]
## [1] 36
## 
## [[36]]
## [1] 37
## 
## [[37]]
## [1] 38
## 
## [[38]]
## [1] 39
## 
## [[39]]
## [1] 40
## 
## [[40]]
## [1] 41
## 
## [[41]]
## [1] 42
## 
## [[42]]
## [1] 43
## 
## [[43]]
## [1] 44
## 
## [[44]]
## [1] 45
## 
## [[45]]
## [1] 46
## 
## [[46]]
## [1] 47
## 
## [[47]]
## [1] 48
## 
## [[48]]
## [1] 49
## 
## [[49]]
## [1] 50
## 
## [[50]]
## [1] 51
## 
## [[51]]
## [1] 52
## 
## [[52]]
## [1] 53
## 
## [[53]]
## [1] 54
## 
## [[54]]
## [1] 55
## 
## [[55]]
## [1] 56
## 
## [[56]]
## [1] 57
## 
## [[57]]
## [1] 58
## 
## [[58]]
## [1] 59
## 
## [[59]]
## [1] 60
## 
## [[60]]
## [1] 61
## 
## [[61]]
## [1] 62
## 
## [[62]]
## [1] 63
## 
## [[63]]
## [1] 64
## 
## [[64]]
## [1] 65
## 
## [[65]]
## [1] 66
## 
## [[66]]
## [1] 67
## 
## [[67]]
## [1] 68
## 
## [[68]]
## [1] 69
## 
## [[69]]
## [1] 70
## 
## [[70]]
## [1] 71
## 
## [[71]]
## [1] 72
## 
## [[72]]
## [1] 73
## 
## [[73]]
## [1] 74
## 
## [[74]]
## [1] 75
## 
## [[75]]
## [1] 76
## 
## [[76]]
## [1] 77
## 
## [[77]]
## [1] 78
## 
## [[78]]
## [1] 79
## 
## [[79]]
## [1] 80
## 
## [[80]]
## [1] 81
## 
## [[81]]
## [1] 82
## 
## [[82]]
## [1] 83
## 
## [[83]]
## [1] 84
## 
## [[84]]
## [1] 85
## 
## [[85]]
## [1] 86
## 
## [[86]]
## [1] 87
## 
## [[87]]
## [1] 88
## 
## [[88]]
## [1] 89
## 
## [[89]]
## [1] 90
## 
## [[90]]
## [1] 91
## 
## [[91]]
## [1] 92
## 
## [[92]]
## [1] 93
## 
## [[93]]
## [1] 94
## 
## [[94]]
## [1] 95
## 
## [[95]]
## [1] 96
## 
## [[96]]
## [1] 97
## 
## [[97]]
## [1] 98
## 
## [[98]]
## [1] 99
## 
## [[99]]
## [1] 100
## 
## [[100]]
## [1] 101
# and
map_dbl(1:100, add_one)
##   [1]   2   3   4   5   6   7   8   9  10  11  12  13
##  [13]  14  15  16  17  18  19  20  21  22  23  24  25
##  [25]  26  27  28  29  30  31  32  33  34  35  36  37
##  [37]  38  39  40  41  42  43  44  45  46  47  48  49
##  [49]  50  51  52  53  54  55  56  57  58  59  60  61
##  [61]  62  63  64  65  66  67  68  69  70  71  72  73
##  [73]  74  75  76  77  78  79  80  81  82  83  84  85
##  [85]  86  87  88  89  90  91  92  93  94  95  96  97
##  [97]  98  99 100 101
# especially the meaning of ~ (speak: lambda)
map(1:5, ~ print("hi"))
## [1] "hi"
## [1] "hi"
## [1] "hi"
## [1] "hi"
## [1] "hi"
## [[1]]
## [1] "hi"
## 
## [[2]]
## [1] "hi"
## 
## [[3]]
## [1] "hi"
## 
## [[4]]
## [1] "hi"
## 
## [[5]]
## [1] "hi"
# more examples
map(1:10, paste, "hi")
## [[1]]
## [1] "1 hi"
## 
## [[2]]
## [1] "2 hi"
## 
## [[3]]
## [1] "3 hi"
## 
## [[4]]
## [1] "4 hi"
## 
## [[5]]
## [1] "5 hi"
## 
## [[6]]
## [1] "6 hi"
## 
## [[7]]
## [1] "7 hi"
## 
## [[8]]
## [1] "8 hi"
## 
## [[9]]
## [1] "9 hi"
## 
## [[10]]
## [1] "10 hi"
is_even <- function(x) x %% 2 == 0
map_lgl(1:10, is_even)
##  [1] FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE
##  [9] FALSE  TRUE
# Note: Often, you don't need a map function
# Because many functions in R are vectorised by default:
x <- 1:10
y <- 1:10
# thus, just write
x + y
##  [1]  2  4  6  8 10 12 14 16 18 20
# instead of
map2_dbl(x,y, `+`)
##  [1]  2  4  6  8 10 12 14 16 18 20

4.7.2 Veranschauliche die folgenden Konzepte anhand von Simulationen auf den vorliegenden Daten

4.7.2.1 Sensitivity: Wie groß ist die Wahrscheinlichkeit abhängig von der Sample Size n mit einem Wilcoxon Rank Sum Test tatsächlich einen p-value <= 0.05 zu erhalten, wenn ein Unterschied existiert?

  • Ziehe 1000 mal ein Sample von 30 je aus wt und ko und teste auf Signifikanz
  • Wie viele der 1000 Versuche sind statistisch signifikant?
  • Wie verändert sich diese Zahl, wenn 10 statt 30 gezogen werden?

4.7.2.2 Lösung

test_samle_wilcox <- function(n) {
  draw_wt <- sample(wt, n, replace = TRUE)
  draw_ko <- sample(ko, n, replace = TRUE)
  wilcox.test(draw_wt, draw_ko, exact = FALSE)$p.value
}
many_p_values <- map_dbl(1:1000, ~test_samle_wilcox(30))
hist(many_p_values)

many_p_values <- map_dbl(1:1000, ~test_samle_wilcox(10))
hist(many_p_values)

many_p_values <- map_dbl(1:1000, ~test_samle_wilcox(4))
hist(many_p_values)

alpha = 0.05
mean(many_p_values <= alpha)
## [1] 0.642
# talk about power, effect size, difference between wilcox and t-test

4.7.2.3 Specificity: Unter der Voraussetzung, dass kein Unterschied zwischen den Bedingungen vorliegt, wie groß ist die Wahrscheinlichkeit, dennoch ein statistisch signifikantes Ergebnis zu erhalten?

  • Stell dir vor, alle Zellen seien wie wt-Zellen
  • Ziehe zwei mal 30 aus den wt-Zellen und lasse einen Wilcoxon Rank Sum Test laufen
  • Widerhole das das Prozedere 1000 mal
  • Wie oft ist das Ergebnis statistisch signifikant?

4.7.2.4 Lösung

test_same_wilcox <- function(n) {
  draw_wt <- sample(wt, n, replace = TRUE)
  draw_ko <- sample(wt, n, replace = TRUE)
  wilcox.test(draw_wt, draw_ko, exact = FALSE)$p.value
}
many_p_values <- map_dbl(1:1000, ~test_same_wilcox(30))
hist(many_p_values)

alpha = 0.05
mean(many_p_values <= alpha)
## [1] 0.065
# this is what we would expect from p-values!