Resampling Techniken

DS3 - Vom experimentellen Design zur
explorativen Datenanalyse & Data Mining

Saskia Otto

Universität Hamburg, IMF

Wintersemester 2024/2025

Lernziele

Am Ende dieser VL- und Übungseinheit werden Sie

  • das sog. Bootstrap-Verfahrens auf reale Datensätze anwenden können, um Konfidenzintervalle für Punktschätzungen zu erstellen.
  • wissen, wie das Bootstrap-Verfahren funktioniert, insbesondere wie Stichproben MIT Zurücklegen aus einer vorhandenen Stichprobe gezogen werden.
  • das Konzepts der Permutation von Beobachtungen kennen und wissen, wie Permutationstests auf verschiedene Arten von Daten anzuwenden sind, um Hypothesen über Populationsparameter zu testen.
  • in der Lage sein, die Ergebnisse von Permutationstests zu interpretieren und zu verstehen, wie statistische Signifikanz in diesem Kontext definiert ist.
  • Permutationstests mit anderen Hypothesentestverfahren zu vergleichen und die Vor- und Nachteile jeder Methode zu verstehen.

Unser Thema heute

Grafik von Nina Garman (Pixabay)

Schätzverfahren und Hypothesentest

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

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

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
  • → Wir sind uns zu 95% sicher, dass der wahre Mittelwert im Bereich 111.8 - 112.4 liegt (112.1 ± 0.3).
  • Unter der Annahme, dass die Verteilung einer symmetrischen Normalverteilung folgt!

Bootstrapping oder ‘Schnürsenkelmethode’

Nicht-parametrische KIs mittels Bootstrapping berechnen

  • Nützliche Methode, wenn Verteilung nicht einer theoretischen Verteilung entspricht (z.B. der z- und t-Verteilung).
  • 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 meist limitiert ist, werden Proben mit Rücklegen entnommen (somit unendlich groß).
    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.
  • Ergebnisse sind recht zuverlässig und mit Computern schnell umzusetzen.
  • Die KI sind nicht unbedingt symmetrisch.

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 = T) # (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








→ Wir sind uns nun zu 95% sicher, dass der wahre Mittelwert im Bereich 111.8 - 112.4 liegt.

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, das 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)
  y_sample <- df$y[indices]
  x_sample <- df$x[indices]
  model <- lm(y_sample ~ x_sample)
  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).

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_orig <- mean(Rennen$Geschw[Rennen$Art == "B"])- mean(
  Rennen$Geschw[Rennen$Art == "A"])
diff_orig
[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)
mean(Rennen_neu$Geschw[Rennen_neu$Art == "B"]) - mean(
  Rennen_neu$Geschw[Rennen_neu$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)
  mystats[i] <- mean(Rennen_neu$Geschw[Rennen_neu$Art == "B"]) -   
    mean(Rennen_neu$Geschw[Rennen_neu$Art == "A"])
}

Zwei-Stichproben-Vergleich | Verteilung

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

p-Werte für 1-/2-seitige Hypothese
# 1-seitig
(n1 <- sum(mystats >= diff_orig))
[1] 475
(p_val_1sided <- n1/10000)
[1] 0.0475
# 2-seitig (mit absoluten Werten)
(n2 <- sum(abs(mystats) >= 
    abs(diff_orig)))
[1] 940
(p_val_2sided <- n2/10000)
[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 # oder: 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_shuffled <- sample(df$sex)
  mod <- lm(body_mass_g ~ flipper_length_mm + sex_shuffled, 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.

Übungen

Übungswoche 5

1. Bootstrap-Konfidenzintervall bei Graugänsen

2. Permutationstest zur Alpakawolle

Vorbereitung @home

Mit Schleifen programmieren

  • Sehen Sie sich die letzte Vorlesung in Data Science 1 zum Thema Schleifen & Simulationen erneut an: 14-advanced-programming.html#/schleifen-simulationen
  • Bearbeiten Sie die Lektion 4 (L04-Schleifen und erste Simulationen) des swirl-Kurses ‘DSB-06-Fortgeschrittene_R_Programmierung’.
  • Beantworten Sie vor der fünften Übungsstunde die Fragen zu Schleifen im Moodle-Quiz.

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.

Kurswebseite: Data Science 3