Resampling Techniken

DS3 - Vom experimentellen Design zur
explorativen Datenanalyse & Data Mining

Saskia Otto & Monika Eberhard

Universität Hamburg, IMF

Wintersemester

Lernziele

Am Ende dieser VL- und Übungseinheit werden Sie

  • mit Schleifen in R wiederholte Zufallsziehungen programmieren können - die technische Grundlage der Resampling-Techniken.
  • das Bootstrap-Verfahren in R durchführen können, d. h. Stichproben mit Zurücklegen ziehen und daraus Konfidenzintervalle für Punktschätzungen berechnen.
  • einen Permutationstest durchführen können, d. h. Beobachtungen ohne Zurücklegen neu zuordnen, um Hypothesen über Populationsparameter zu testen und die statistische Signifikanz des Ergebnisses zu bestimmen.
  • die Ergebnisse beider Verfahren interpretieren und beurteilen können, wann ein Resampling-Ansatz gegenüber parametrischen Tests vorzuziehen ist.

Unser Thema heute

Grafik von Nina Garman (Pixabay)

Schätzverfahren und Hypothesentest | 1

Beide Verfahren der inferenziellen Statistik beruhen auf unterschiedlichen Stichprobenverteilung:

  1. Schätzverfahren
    • Verwendet die Stichprobenverteilung einer Schätzung.
    • Sie wird verwendet, um Standardfehler und Konfidenzintervalle zu erhalten.
    • Die meisten Methoden gehen davon aus, dass die Stichprobenverteilung annähernd normal ist.
  1. Hypothesentest
    • Verwendet die Null-Stichprobenverteilung (oder Nullverteilung): die Wahrscheinlichkeitsverteilung einer Teststatistik, wenn die Nullhypothese wahr ist.
    • Häufig werden die t-, F-, \chi^2- und Normalverteilungen zur Annäherung an die Nullverteilungen verwendet, aus denen die P-Werte berechnet werden.

Schätzverfahren und Hypothesentest | 2

Frage:

Was ist zu tun, wenn die Annahmen der besten verfügbaren Methode verletzt werden und wir nicht auf parametrische Tests und allgemeine lineare Modelle zurückgreifen können?

Antwort: Computerintensive Methoden

Ein Ansatz, bei dem die Leistungsfähigkeit des Computers genutzt wird, um eine Stichprobenverteilung zu erstellen:

  1. Schätzverfahren: Bootstrapping
  2. Hypothesentest: Permutationstests

→ 1 Datensatz, 10.000 Wiederholungen: Schleifen

Noch nie mit Schleifen gearbeitet? Kein Problem — die führe ich gleich ein.

Zur Erinnerung: Die 4 Grundelemente

Die 4 Grundelemente

von Objekt-orientierten Sprachen wie R, Python und Javascript

  • Variablen
  • Bedingungen
    • Eine Bedingung ist eine Abfrage; trifft diese zu, so wird das, was innerhalb der Bedingung ist, ausgeführt.
  • Schleifen
    • Eine Schleife enthält Code, der so lange ausgeführt wird, wie die Bedingung, die zur Schleife gehört, wahr ist.
  • Funktionen
    • Sammlung von Code, die jederzeit über den jeweiligen Funktionsnamen aufgerufen werden kann.

Bedingungen | 3 Typen in R

Wenn bei der Abfrage TRUE rauskommt, wird die (erste) Anweisung ausgeführt:

if.. Anweisung

if (Abfrage) {
   mache folgendes
} 

x <- 10L
if ( is.integer(x) ) {
  print(
    "Yes, x is an integer!")
}
[1] "Yes, x is an integer!"
if ( is.double(x) ) {
  print(
    "Yes, x is a double!")
}

→ hier passiert nichts weil die Abfrage FALSE ergibt!

if..else..

if (Abfrage) {
   mache dies
} else {
   mache das
} 

y <- 12.5
if (y < 6) {
  y <- 0
} else {
  y <- 1
}
y
[1] 1

ifelse() Funktion

ifelse (Abfrage, mache dies, mache das) 

z <- 1:3
ifelse(z == 2, z^2, 0)
[1] 0 4 0

Es können neben Operatoren auch Funktionen in der Abfrage verwendet werden, solange diese ein TRUE oder FALSE zurückgeben.

Bedingungen | Anwendungsbeispiel

Farbeinstellung bei Basisplots mit ifelse()

plot(x = iris$Sepal.Length, y = iris$Petal.Length, 
  col = ifelse(iris$Species == "setosa", "blue", "red"),
  pch = 20,  cex = 1, xlab = "Länge Kelchblatt", ylab = "Länge Bütenblatt", 
  main = "setosa vs. versicolor/virginica")

Schleifen

Zur Erinnerung

Eine Schleife enthält Code, der so lange ausgeführt wird, wie die Bedingung, die zur Schleife gehört, wahr ist. Der Code kann also z.B. zehn Mal ausgeführt werden, das hat den Vorteil, dass der Code nicht zehn Mal geschrieben werden muss.

  • Schleifen gehören neben Bedingungen zu den Kontrollstrukturen → sie kontrollieren den Programmablauf, indem Programmteile nur unter gewissen Bedingungen ausgeführt oder wiederholt werden.
  • Der Codeblock an Anweisungen innerhalb einer Schleife wird auch ‘Schleifenrumpf’ oder ‘Schleifenkörper’ (im Englischen ‘loop body’) genannt.
  • Schleifen, deren Schleifenbedingung immer zur Fortsetzung führt oder die keine Schleifenbedingung haben, sind Endlosschleifen. Sie sollten vermieden werden!
  • Schleifen können beliebig verschachtelt werden.

Schleifentypen in R

  • Zählschleifen (for): führen eine bestimmte Anzahl an Wiederholungen aus, die durch einen Zähler oder einen Index kontrolliert wird. Dieser erhöht sich mit jeder Iteration.
  • Kopfgesteuerte/vorprüfende (while) vs. fußgesteuerte/nachprüfende Schleifen (repeat): basieren auf ein Einsetzen und einer Verifikation durch eine logische Bedingung. Die Bedingung wird zu Beginn oder zum Ende des Schleifenkonstrukts getestet.

for Schleife


for Schleife


for Schleife

Ausführung in R
for (i in c(2,4,7) ) {
  print(i)
}
[1] 2
[1] 4
[1] 7
i
[1] 7

for Schleife - Stile

Bestes Vorgehen für die Ergebnisausgabe

  • Bevor eine Schleife gestartet wird, sollte immer genug Platz für das Ergebnis bereitgestellt werden.
  • Wenn die for Schleife mit jeder Wiederholung mit z.B. c() wächst, wird die for Schleife sehr langsam werden:
Objekt wachsen lassen
grow_obj <- function(x){
  y <- numeric()
  for (i in 1:x) { 
    y <- c(y, i)
  }
}
Indexierung
# Better: Indexing
index_obj <- function(x){
  y <- numeric(length(x))
  for (i in 1:x){ 
    y[i] <- i
  }
}

Hier wird ein leerer Vektor mit der Länge des Zählers erstellt, bevor die Schleife startet.

Bestes Vorgehen für die Ergebnisausgabe

Testen wir die Geschwindigkeit der beiden Funktionen
microbenchmark::microbenchmark(unit = "ms", times = 5,
    grow_obj(500), grow_obj(5000), index_obj(500), index_obj(5000) )
Unit: milliseconds
            expr       min        lq       mean    median        uq       max
   grow_obj(500)  0.229026  0.237103  0.2458032  0.243294  0.253052  0.266541
  grow_obj(5000) 28.640837 28.880523 30.4578422 30.028072 31.081567 33.658212
  index_obj(500)  0.061295  0.063468  0.0670842  0.063878  0.065354  0.081426
 index_obj(5000)  0.590113  0.599912  0.9261654  0.601101  0.616066  2.223635
 neval cld
     5  a 
     5   b
     5  a 
     5  a

In der Spalte mean ist zu sehen, dass die Variante mit der ‘Indexierung’ viel schneller ist als die mit wachsenden Objekten, besonders wenn viele Wiederholungen gemacht werden!

Your turn …

Quiz 1-4 | Bedingungen & Schleifen

Q1 - Bedingungen

Der Körper einer Bedingung (wie auch Schleife und Funktion) wird immer in geschweifte Klammern gesetzt, damit er über mehrere Zeilen gehen kann.

In der Abfrage muss ein Operator stehen, welcher die Abfrage wahr macht, denn der Codeblock in der Bedingung wird ja ausgeführt. Richtig ist daher:

x <- 25
if (x <= 25) {
  print('Yep, this is correct')
}
[1] "Yep, this is correct"

Q2 - Bedingungen

Anhand der zweimaligen geschweiften Klammern wird deutlich, dass es sich um eine ‘if..else..’ Bedingung handelt.

An der Ausgabe 0 ist zu erkennen, dass der Alternativbefehl (in ‘else’) ausgeführt wurde. In die zweite Lücke muss daher ein Operator, bei der die Abfrage FALSE ergibt:

i <- 1
if (i != 1) {
  i + 1
} else {
  i - 1
}
[1] 0

Q3 - Schleifen

Anhand der Abfrage ist erkennbar, dass es sich um eine ‘while..’ Schleife handeln muss. Auch hier gilt wieder: der Schleifenblock gehört in geschweifte Klammern!

Und da eine ‘while..’ Schleife nur mehrere Iterationen durchgeht, wenn der Input in der Abfrage sich ändert, muss sich hier a im Schleifenkörper um den Wert erhöhen, den der Output sich erhöht (weil dieser a darstellt.)

a <- 3
while (a < 12) {
  print(a)
  a <- a + 4
}
[1] 3
[1] 7
[1] 11

Q4 - Schleifen

Hier ist keine Abfrage jeweils zu sehen, sondern es wird ein Zählerindex definiert. Also muss es sich um die Zählschleife (‘for..’) handeln. Und zwar 2 ineinander verschachtelte. Wichtig ist hierbei, dass der Zähler unterschiedlich ist. Die erste Schleife hat meist den Index i, die nächste j, etc.).

Der innerste Schleifenkörper macht folgendes: es sollen die beiden Indizes zusammengefügt werden (mittels Funktion paste()) unter Verwendung des Trennzeichen ‘/’ . Damit dies auch wirklich in der Konsole sichtbar ist, braucht es aber die Funktion print(), welche als Input paste() enhält.

for (i in c(1,3)) {
  for (j in c(5,7)) {
    print(paste(i, j, sep = '/'))
  }
}
[1] "1/5"
[1] "1/7"
[1] "3/5"
[1] "3/7"

Bootstrapping

Punkt- und Intervallschätzung

Zur Erinnerung

  • Auf Basis einer Zufallsprobe wird der Bereich geschätzt, wo der gesuchte Populationsparameter liegen könnte.
  • Ausgangspunkt ist immer eine Punktschätzung → dann wird ein (symmetrisches) Intervall bestimmt, das Konfidenzintervall KI
  • Konfidenz wird als wiederholte Stichprobe interpretiert.
  • Kann für jeden Parameter (z.B. Mittelwert, Varianz, ..) berechnet werden.
  • Die Breite des Intervalls hängt ab
    • vom Stichprobenumfang, der Varianz der Stichprobe und
    • der festgelegten Wahrscheinlichkeit (= Konfidenzniveau) → üblich: 90%, 95%, 99%

KI vom Mittelwert wenn n > 30

KI_{95\%} = 1.96\cdot\sqrt{\frac{\sigma^2}{n}}

KI vom Mittelwert wenn n < 30

KI_{95\%} = t_{(\alpha/2,df)}\cdot\sqrt{\frac{s^2}{n}}

Intervallschätzung | 1

Zur Erinnerung ein Beispiel aus DS2

  • Wie präzise ist diese Schätzung von 117mm?
  • Kann es denn sein, dass der wahre Populationsmittelwert auch 100mm oder 140mm ist und wir einfach Pech mit der Probe hatten?
  • In welchem Bereich liegt der wahre Mittelwert höchstwahrscheinlich?

Intervallschätzung | 2

Zur Erinnerung ein Beispiel aus DS2

Die Population

mean(population)
[1] 112.0038

Die Stichprobe

x <- sample(population, 50)
mean(x)
[1] 112.1132
hist(x, main = "Stichprobe (N=50)")

Intervallschätzung | 3

Parametrisches KI basierend auf der Normalverteilung

(x_mean <- mean(x))
[1] 112.1132
(x_se <- sd(x)/sqrt(50))
[1] 0.1699562
(z_lower <- qnorm(p = 0.025, mean = 0, sd = 1))
[1] -1.959964
(z_upper <- qnorm(p = 0.975, mean = 0, sd = 1))
[1] 1.959964
(CI_lower <- z_lower*x_se)
[1] -0.3331081
(CI_upper <- z_upper*x_se)
[1] 0.3331081
  • → Der Populationsmittelwert liegt mit 95%-iger Konfidenz im Bereich 111.8 - 112.4 (112.1 ± 0.3) mm.
  • Unter der Annahme, dass die Verteilung einer symmetrischen Normalverteilung folgt!

Bootstrapping oder ‘Schnürsenkelmethode’ | 1

Nicht-parametrische KIs mittels Bootstrapping berechnen

Prinzip:

  1. Die Verteilung der Stichprobe dient als Grundlage einer ‘Pseudopopulation’.
  2. Aus dieser ‘Pseudopopulation’ entnehmen wir wiederholt Stichproben und berechnen jeweils die Kenngröße.
  3. Da der Stichprobenumfang begrenzt ist, wird mit Zurücklegen gezogen - so entstehen beliebig viele Wiederholungsstichproben gleicher Größe aus derselben Pseudopopulation.
  4. Aus der Grundgesamtheit der Kennwerte (z.B. der Mittelwerte) berechnen wir das untere und obere 2.5 % (oder 0.5 %,..) Quantil als unsere Konfidenzgrenze.

Bootstrapping oder ‘Schnürsenkelmethode’ | 2

Nicht-parametrische KIs mittels Bootstrapping berechnen

  • Nützliche Methode, wenn Verteilung nicht einer theoretischen Verteilung entspricht (z.B. der z- und t-Verteilung).
  • Die KI sind nicht unbedingt symmetrisch.
  • Ergebnisse sind recht zuverlässig und mit Computern schnell umzusetzen - allerdings nur, wenn die Stichprobe die Population einigermaßen gut repräsentiert.
    • Faustregel: n ≳ 20–30
    • Bei sehr kleinem n wird die Pseudopopulation zu grob und die KI-Grenzen können verzerrt sein.

Bootstrapping | Berechnung

Schritt 1-3: die Schleife
set.seed(123)
it <- 10000 # 10000 Iterationen
boot_mean <- numeric(it) 
for (i in 1:it){
  x_sample <- sample(x, replace = TRUE) # (N auch 20)
  boot_mean[i] <- mean(x_sample)
}
Schritt 4: die Quantilen = KI
# 2.5% und 97.5% Quantilen = 95% Konf.grenzen:
quantile(x = boot_mean, probs = c(0.025, 0.975))
    2.5%    97.5% 
111.7824 112.4487 
# KI:
quantile(boot_mean, c(0.025, 0.975)) - x_mean 
      2.5%      97.5% 
-0.3307953  0.3355195








→ Das 95%-Bootstrap-Konfidenzintervall für den Populationsmittelwert reicht von 111.8 bis 112.4.

Bootstrapping | Vergleich

Schätzung von KI für Regressionsparameter

Wenn Annahmen nicht erfüllt sind

Die Modellparameter
mod_lm <- lm(y ~ x, df)
summary(mod_lm)$coefficients
              Estimate Std. Error   t value     Pr(>|t|)
(Intercept) -55.960944  22.550127 -2.481624 1.761662e-02
x             8.135234   1.464779  5.553899 2.329342e-06
Parametrische Konfidenzintervalle
(ci_param <- c(
  lower_ci = 8.135234 - 1.96 * 1.464779,
  upper_ci = 8.135234 + 1.96 * 1.464779))
 lower_ci  upper_ci 
 5.264267 11.006201
  • Der Steigungsparameter b der Stichprobe als Punktschätzer für die Population beträgt 8.1.
  • Gegeben, dass die Daten normalverteilt und varianzhomogen sind → können wir uns zu 95 % sicher sein, dass bei einer Wiederholung der Datenerhebung die Regressionssteigung für diese neuen Daten zwischen 5.3 und 11.0 liegen würde.
Code 
df$res <- residuals(mod_lm)
df$fit <- fitted(mod_lm)
p1 <- ggplot(df, aes(y)) + 
  geom_histogram(bins = 10, fill = "grey50", colour = "grey10") +
  ggtitle("Histogramm von Y")
p2 <- ggplot(df, aes(x, y)) + geom_point() +
  geom_smooth(method = "lm", se = TRUE) + 
  ggtitle("Beziehung Y ~ X")
p3 <- ggplot(df, aes(res)) + 
  geom_histogram(bins = 10, fill = "grey50", colour = "grey10") +
  ggtitle("Histogramm der Residuen")
p4 <- ggplot(df, aes(fit, res)) + geom_point() + 
  geom_hline(yintercept = 0) +
  ggtitle("Residuen vs. fitted")
gridExtra::grid.arrange(p1, p2, p3, p4)

Schätzung von KI für Regressionsparameter

Bootstrapping

set.seed(32)
it <- 10000
boot_b <- numeric(it)
for (i in 1:it){
  n <- nrow(df)
  indices <- sample(1:n, replace = T)
  model <- lm(y ~ x, data = df[indices, ])
  boot_b[i] <- coef(model)[2]
}
quantile(boot_b,c(0.025,0.975))
     2.5%     97.5% 
 4.663306 11.757071

→ Der Konfidenzbereich liegt nun zwischen 4.7 und 11.8.

Permutationstest

Permutationstests | 1

  • Ähnlich wie Bootstrapping, mit dem Unterschied, dass
    • beim Permutationstest Stichproben OHNE Zurücklegen gezogen werden (d.h., wenn ein Wert ausgewählt wurde, kann er nicht erneut ausgewählt werden, so dass kein Wert ein Duplikat sein kann).
    • Dabei werden die Werte einfach neu gemischt.
  • Im Falle einer univariaten Statistik (z. B. Mittelwert) ändert sich dadurch nichts.
  • Bei zwei oder mehr Variablen bzw. Gruppen ändert sich durch die Umordnung einer Variablen allerdings die Teststatistik, z. B. die Korrelation oder Regression.

Permutationstests | 2

Prinzip der Permutation

Man generiert sich mit Hilfe des Computers seine eigene Zufalls- oder Nullverteilung und prüft dann, ob sich die gefundene Verteilung von dieser Zufallsverteilung unterscheidet.

  • Annahme: Die gemessenen Werte sind die Werte, die vorkommen können.
  • Frage: Mit welcher Häufigkeit tritt die gefundene Verteilung bei den oben generierten Zufallsverteilungen auf?
  • Signifikanz: Die Häufigkeit, mit der die gefundene Verteilung bei den durch zufällige Verteilung generierten Gruppen auftritt, entspricht der Irrtumswahrscheinlichkeit.

Demonstration aus dem Internet

Webseite: https://www.jwilber.me/permutationtest/

Was ist besser?

Nicht-parametrische Tests

  • Rang-basierte Tests, wie z. B. der Mann-Whitney-U-Test für zwei Stichproben, stellen Permutationstests dar:
    • Daten werden durch Ränge ersetzt und diese permutiert, um eine Nullverteilung zu erzeugen.
    • Die genaue Wahrscheinlichkeitsverteilung der U-Statistik ist bekannt.
    • Durch die Ersetzung der Daten durch ihre Ränge gehen jedoch Informationen verloren.

Permutationstest

  • Statt Daten durch die Ränge zu ersetzen können sie direkt permutiert werden.
  • Permutationstests haben bei kleinem Stichprobenumfang eine geringere Teststärke als parametrische Tests, aber sie sind aussagekräftiger als z.B. der Mann-Whitney-U-Test.
  • Bei großem Stichprobenumfang haben sie eine ähnliche Teststärke wie parametrische Tests.

Annahmen von Permutationstests

  • Zufällige Stichproben.
  • Bei Gruppenvergleichen der Mittelwert oder Mediane sollte die Verteilung der Variablen in jeder Population die gleiche Form haben.
  • Es gilt jedoch: Permutationstests sind robust gegenüber Abweichungen von der Gleichverteilungsannahme, wenn die Stichprobengröße groß ist (stärker als beim Mann-Whitney-U-Test).

2 Beispiele für Permutationstest

Beispiel 1: Zwei-Stichproben-Vergleich

Schneckenrennen

Unterscheiden sich zwei Schneckenarten in ihrer Laufgeschwindigkeit (in m/Stunde)?

Aus: https://www.globaltimes.cn/page/201807/1111861.shtml

Zwei-Stichproben-Vergleich | Daten

A <- c(4,2,7,9,2,8)
B <- c(7,10,9,5,8,11)
Rennen <- data.frame(
  Art = c(rep("A", 6), rep("B", 6)),
  Geschw = c(A,B)
)

Rennen |>
  group_by(Art) |>
  summarise(
    Mittelwert = mean(Geschw), 
    Standardabweichung = sd(Geschw)
  )
# A tibble: 2 × 3
  Art   Mittelwert Standardabweichung
  <chr>      <dbl>              <dbl>
1 A           5.33               3.08
2 B           8.33               2.16

Zwei-Stichproben-Vergleich | t-Test

t.test(Geschw ~ Art, Rennen) # 2-seitige Hypothese

    Welch Two Sample t-test

data:  Geschw by Art
t = -1.9547, df = 8.9658, p-value = 0.08247
alternative hypothesis: true difference in means between group A and group B is not equal to 0
95 percent confidence interval:
 -6.473935  0.473935
sample estimates:
mean in group A mean in group B 
       5.333333        8.333333 
t.test(Geschw ~ Art, Rennen, alternative = "less") # 1-seitige Hypothese

    Welch Two Sample t-test

data:  Geschw by Art
t = -1.9547, df = 8.9658, p-value = 0.04123
alternative hypothesis: true difference in means between group A and group B is less than 0
95 percent confidence interval:
       -Inf -0.1853508
sample estimates:
mean in group A mean in group B 
       5.333333        8.333333

Zwei-Stichproben-Vergleich | Permutation 1

Permutationstest - einmalig

Differenz beider Mittelwerte als Teststatistik
diff_obs <- with(Rennen, mean(Geschw[Art=="B"]) - mean(Geschw[Art=="A"]))
diff_obs
[1] 3
Teststatistik bei neuer Stichprobe
set.seed(321)
Rennen_neu <- Rennen
# Wir können die gruppierende Variable oder die Antwortvariable mischen
Rennen_neu$Geschw <- sample(Rennen_neu$Geschw, replace = FALSE) # (default)
with(Rennen_neu, mean(Geschw[Art=="B"]) - mean(Geschw[Art=="A"]))
[1] -0.3333333

Zwei-Stichproben-Vergleich | Permutation 2

Permutationstest - 10000fach

set.seed(321)
it <- 10000 # 10000 Iterationen
mystats <- numeric(it) 
for (i in 1:it){
  # Schleifenkörper: wir kopieren den Code von eben
}

Zwei-Stichproben-Vergleich | Permutation 2

Permutationstest - 10000fach

set.seed(321)
it <- 10000 # 10000 Iterationen
mystats <- numeric(it) 
for (i in 1:it){
  # Schleifenkörper: wir kopieren den Code von eben
  Rennen_neu <- Rennen
  Rennen_neu$Geschw <- sample(Rennen_neu$Geschw, replace = FALSE)
  mystats[i] <- with(Rennen_neu, 
    mean(Geschw[Art=="B"]) - mean(Geschw[Art=="A"]))
}

Zwei-Stichproben-Vergleich | Verteilung

Verteilung der 10000 Teststatistiken
hist(mystats)
abline(v = diff_obs, col = "red")

p-Werte für 1-/2-seitige Hypothese
### 1-seitig
(n1 <- sum(mystats >= diff_obs))
[1] 475
(p_val_1sided <- n1/it)
[1] 0.0475
### 2-seitig (mit absoluten Werten)
(n2 <- sum(abs(mystats) >= 
    abs(diff_obs)))
[1] 940
(p_val_2sided <- n2/it)
[1] 0.094
# direkte Berechnung von p mit mean()
mean(abs(mystats) >=  abs(diff_obs))
[1] 0.094

Zwei-Stichproben-Vergleich | Interpretation

  • Bei 475 von 10000 Permutationen sind die Differenzwerte zwischen den beiden Mittelwerten gleich 3 oder größer (also als die Teststatistik der Originaldaten).
  • Damit ist die Wahrscheinlichkeit, dass eine Differenz von ≥ 3 auftritt = 4.8 %.
  • Für die 2-seitige Hypothese vergleichen wir die absoluten Werte!
  • D.h., die Wahrscheinlichkeit eine solche Teststatistik wie die beobachtete zu erhalten - wenn die 1- bzw. 2-seitige H0 richtig ist - liegt bei 4.8 bzw. 9.4 %. = die Wahrscheinlichkeit, dass die beobachteten Daten unter der Annahme der H0 auftreten, liegt bei <5 bzw. < 10%.

Beispiel 2: ANCOVA

Schnabellänge bei Adélie Pinguinen

Das Modell
summary(mod_adelie)

Call:
lm(formula = body_mass_g ~ flipper_length_mm + sex, data = adelie)

Residuals:
    Min      1Q  Median      3Q     Max 
-662.85 -213.18   -6.27  207.42  743.73 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        305.087    755.580   0.404    0.687    
flipper_length_mm   16.314      4.019   4.059 8.08e-05 ***
sexmale            599.343     52.246  11.472  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 295.1 on 143 degrees of freedom
Multiple R-squared:  0.5918,    Adjusted R-squared:  0.5861 
F-statistic: 103.6 on 2 and 143 DF,  p-value: < 2.2e-16
Muster in der Varianz?
plot(residuals(mod_adelie) 
  ~ adelie$sex)

Beispiel 2: ANCOVA

Schnabellänge bei Adélie Pinguinen

Hier gibt es 2 Möglichkeiten der Permutation:

  1. Permutation des Faktors:
    • Testet den Effekt des Faktors auf die Antwortvariable, nachdem die Kovariate kontrolliert wurde. Fokus auf den Faktoreffekt.
  2. Permutation der Antwortvariable:
    • Testet die gesamte Modellanpassung, d. h. ob der Faktor und die Kovariate zusammen die Antwortvariable erklären. Fokus auf die Gesamtanpassung des Modells.

ANCOVA | Permutation des Faktors 1

Wir extrahieren zuerst die F-Statistik aus der drop1() Funktion.

Speichern der Original-Statistik
mod_drop1 <- drop1(mod_adelie, test = "F")
mod_drop1
Single term deletions

Model:
body_mass_g ~ flipper_length_mm + sex
                  Df Sum of Sq      RSS    AIC F value    Pr(>F)    
<none>                         12450275 1663.6                      
flipper_length_mm  1   1434486 13884760 1677.5  16.476 8.077e-05 ***
sex                1  11457596 23907871 1756.9 131.598 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
f_sex_orig <- 131.598 # besser: mod_drop1$`F value`[3]

ANCOVA | Permutation des Faktors 2

Wir schreiben die Schleife
set.seed(1234)
it <- 10000 # 10000 Iterationen
mystats_sex <- numeric(it)
for (i in 1:it){
  # Schleifenkörper
}

ANCOVA | Permutation des Faktors 2

Wir schreiben die Schleife
set.seed(1234)
it <- 10000 # 10000 Iterationen
mystats_sex <- numeric(it)
for (i in 1:it){
  # Schleifenkörper
  df <- adelie
  df$sex_perm <- sample(df$sex, replace = FALSE)
  mod <- lm(body_mass_g ~ flipper_length_mm + sex_perm, df)
  mystats_sex[i] <- drop1(mod, test = "F")$`F value`[3]
}

ANCOVA | Permutation des Faktors 3

Verteilung der Teststatistik

Code 
hist(mystats_sex, xlim = c(0,140))
abline(v = f_sex_orig, 
  col = "blue", lwd = 2)



p-Wert (2-seitige Hyp.)
(n_s <- sum(
  abs(mystats_sex) >= abs(f_sex_orig)))
[1] 0
(p_val_s <- n_s / it)
[1] 0

INTERPRETATION:

  • die Wahrscheinlichkeit, dass die beobachteten Daten unter der Annahme der H0 auftreten, liegt beim Faktor ‘sex’ unter 5 %.
  • Wir können daher die H0 ablehnen, es gibt einen signifikanten Effekt des Geschlechts.

Your turn …

Quiz 5+6 | Bootstrap-KI für den Mittelwert

Quiz 7 | Permutation 2-Stichproben-Vergleich

Durchführung

Quiz 8 | Permutation 2-Stichproben-Vergleich

Interpretation

1- und 2-seitiger p-Wert
p_onesided
[1] 0.0524
p_twosided
[1] 0.1004

Bootstrap & Permutation bei linearer Regression

Bei der folgenden einfachen linearen Regression sind nicht alle Modellannahmen erfüllt, wie Residuenplots (rechts) zeigen:

mod <- lm(counts ~ area, data = birds)
summary(mod)$coefficients
             Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 48.907219  10.776106 4.538487 1.373153e-05
area         5.403289   1.780359 3.034943 2.960263e-03
(ci_param <- c(
  lower_ci = 5.403289 - 1.96 * 1.780359,
  upper_ci = 5.403289 + 1.96 * 1.780359))
lower_ci upper_ci 
1.913785 8.892793

Wir wollen daher für den Steigungsparameter das Konfidenzintervall per Bootstrapverfahren und den p-Wert per Permutationstest berechnen.

Quiz 9 | Lineare Regression

Durchführung beider Verfahren

Quiz 10 | Lineare Regression

Interpretation

Vergleich der 95%-KIs
ci_param
lower_ci upper_ci 
1.913785 8.892793 
ci_boot
    2.5%    97.5% 
1.133575 8.856406 
Vergleich der p-Werte
p_param_twosided
[1] 1.373153e-05
p_perm_twosided
[1] 0.002

Übungen

Übungsaufgaben

Alle Aufgaben der Übungswoche finden Sie im Handbuch, gegliedert in drei Abschnitte:

  • Vorbereitung @home — Datenexploration oder Konzeptwiederholung, inkl. Moodle-Vorbereitungsquiz
  • Aufgaben @Übungsstunde — praktische R-Aufgaben in den 90 Minuten, inkl. Moodle-Übungsquiz und CodeRunner-Aufgaben im Anschluss
  • Challenges @Übungsstunde & @home — freiwillige Vertiefungsaufgaben

→ Die Vorbereitung bearbeiten Sie bitte vor der Übungsstunde.

  • Die Übungsnotebooks und Daten liegen in Ihrem R-Projektordner Data-Science-3.

Fragen..??


Total konfus?



Buchkapitel zum Nachlesen

Total gelangweilt?

Dann testen Sie doch Ihr Wissen in folgendem Abschlussquiz…

Abschlussquiz

Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de

Creative Commons License
Diese Arbeit is lizenziert unter einer Creative Commons Attribution-ShareAlike 4.0 International License mit Ausnahme der entliehenen und mit Quellenangabe versehenen Abbildungen.