Statistische Verfahren sind dazu da, Antworten auf Fragen zu finden. Ist Medikamt X besser als Medikamt Y? Oder, ist Diät P besser als Diät S? Solche Fragen kennt jeder, sie trifft man häufig in Alltagsgesprächen an und sie bewegen sich noch nicht in Bereichen, bei denen Statistiklaien am liebsten sofort den Raum verlassen möchten, sobald der Vortragende mit Begriffen wie statistische Kennwerte, Stichprobenverteilung oder Standardabweichung um sich wirft.

Um solche Unterschiedsfragen (ist X besser als Y) zu beantworten braucht es nichtsdestotrotz statistische Verfahren. Diese Verfahren müssen auch gar nicht kompliziert sein. Stelle dir vor, du möchtest Wissen, ob Diät 1, Diät 2 oder Diät 3 besser ist, um Gewicht abzunehmen. Du teilst 60 Leute in je 3 Gruppen (20 Menschen pro Gruppe). Gruppe 1 muss Diät 1 durchziehen, Gruppe 2 muss Diät 2 durchziehen und Gruppe 3 muss Diät 3 durchziehen. Jede Gruppe für exakt 6 Wochen. Bevor der Versuch beginnt, wird nochmal jeder gewogen. Am Ende des Versuchs müssen wieder alle auf die Waage. Welche Diät ist nun am besten? Die einfachste Antwort wäre diejenige Gruppe, die am meisten abgenommen hat. Lass uns das prüfen.

R, R-Studio und das Tidyverse

Ich werde alle Analysen gleich mit R durchspielen. R ist eine Programmiersprache, die für statistische Auswertungen geschrieben wurde. Wenn du mitmachen möchtest, brauchst du daher R und R-Studio. Den Datensatz für dieses Beispiel nehmen wir uns von dieser Webseite (lade dir die csv-Datei zum Thema Diät herunter). Solltest du es noch nicht gemacht, haben installiere zunächst das Paket tidyverse und das Paket gghighlight:

install.packages("tidyverse")
install.packags("gghighlight")

Anschließend können wir die Pakete laden und den Datensatz einlesen:

library(tidyverse)
library(gghighlight)
diet_dataset <- read_csv("stcp-Rdataset-Diet.csv")

Wenn du noch nicht mit R vertraut bist, stelle sicher, dass du das richtige Arbeitsverzeichnis bestimmst, bevor du die Daten einliest. Am einfachsten gibst du folgenden Befehl in die Konsole ein und wählst denjenigen Ordner aus, in der sich die stcp-Rdataset-Diet.csv Datei befindet.

setwd()

Die Daten müssten nun geladen sein und du kannst dir mit glimpse die Daten ansehen:

> glimpse(diet_dataset)
Observations: 78
Variables: 7
$ Person       <int> 25, 26, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 27, 28, 29, 30, 31, 32, 33, 34,...
$ gender       <int> NA, NA, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
$ Age          <int> 41, 32, 22, 46, 55, 33, 50, 50, 37, 28, 28, 45, 60, 48, 41, 37, 44, 37, 41, 43, 20, 51...
$ Height       <int> 171, 174, 159, 192, 170, 171, 170, 201, 174, 176, 165, 165, 173, 156, 163, 167, 174, 1...
$ pre.weight   <int> 60, 103, 58, 60, 64, 64, 65, 66, 67, 69, 70, 70, 72, 72, 72, 82, 58, 58, 59, 61, 62, 6...
$ Diet         <int> 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...
$ weight6weeks <dbl> 60.0, 103.0, 54.2, 54.0, 63.3, 61.1, 62.2, 64.0, 65.0, 60.5, 68.1, 66.9, 70.5, 69.0, 6...

Wir beschäftigen uns jetzt nur noch mit zwei dieser Variablen: weight6weeks und Diet. Diet ist nichts anderes als eine Variable, die die Zahlen 1, 2 und 3 enthält. 1 steht für Diät 1 und so weiter. weight6weeks kennzeichnet das Gewicht der Menschen nach 6 Wochen.

Wir sind allerdings an dem Gewichtsunterschied nach 6 Wochen interessiert. Dazu erstellen wir eine neue Variable:

diet_dataset <- diet_dataset %>%
  mutate(
    weight_diff = weight6weeks - pre.weight
  )

mutate ist eine Funktion, mit der wir neue Variablen erstellen können.

Jetzt sind wir bereit heraus zu finden, ob die Diäten nun unterschiedlich erfolgreich sind.

Macht die Diät einen Unterschied? Eine einfache Methode, die Frage zu beantworten

Was wäre die einfachste Methode, zu prüfen, ob sich die Gruppen unterscheiden? Wir vergleichen die Mittelwerte.

x=i=1nxin\overline{x} = \frac{\displaystyle\sum_{i=1}^{n}x_i}{n}

Oder vereinfacht gesagt, wir summieren alle Werte und teilen diese Summe durch die Anzahl der Werte. Lass uns zunächst den Mittelwert der drei Gruppen bestimmen.

(mean_between_groups <- diet_dataset %>%
  group_by(Diet) %>%
  summarise(
    mean = mean(weight_diff)
  ))
# A tibble: 3 x 2
   Diet  mean
  <int> <dbl>
1     1 -3.3 
2     2 -3.03
3     3 -5.15

Mmmmh, ok Gruppe 3 scheint am meisten abgenommen zu haben. Ist also Diät 3 am besten? Oder anders gefragt, ab welchem Unterschied ist eine Diät besser als die andere? Oder, kann es nicht einfach sein, dass es die Unterschiede zufälligen Schwankungen unterliegen? Überlge dir das mal selber. Wenn jetzt der Mittelwert der Diät 3 bei -4 läge, würdest du dann sagen, dass die Diät besser ist als die andere? Was wäre bei einem Mittelwert von -3.5? Du merkst, die Entscheidung ist ein wenig willkürlich. Wir brauchen in der Regel einen Cut-Off-Wert, anhand wir entscheiden, ob eine Diät einen Unterschied ausmacht. Diesen Cut-Off-Wert bestimmen wir durch andere Methoden.

Eine andere Methode Unterschiede, zwischen Gruppen zu bestimmen

Anstatt den Mittelwert zu nehmen, um heraus zu finden, ob sich die Gruppen unterscheiden, könnten wir auch fragen, ob diese Unterschiede zufälligen statistischen Schwankungen unterliegen, oder ob es tatsächlich unwahrscheinlich ist, dass sich 3 Gruppen so stark voneinander unterscheiden. Um zu verstehen, was ich damit meine, ein kleines Beispiel.

Ein Würfelbeispiel

Nehmen wir an, du hast einen 6-seitigen Würfel. Du möchtest wissen, ob dieser Würfel auch wirklich jede Zahl mit gleicher Wahrscheinlichkeit anzeigt. Das bedeutet, jede Zahl sollte zu (1/6 * 100)-prozentiger, also 16.67%-iger Wahrscheinlichkeit erscheinen. Wenn wir allerdings 6 mal würfeln ist es sehr unwahrscheinlich, dass jede Zahl genau einmal vorkommt. Das Vorkommen der Zahlen unterliegt bestimmten statistischen Schwankungen. Sollten wir allerdings 10000 mal würfen, können wir annehmen, dass jede Zahl zu etwa 16.67% auftritt.

Ein nicht-manipulierter Würfel

Zunächst würfeln wir mit einen nicht-manupulierten Würfen 10000 mal:

set.seed(100)
my_sample <- sample(1:6, 10000, TRUE)

Mit set.seet(100) bekommst du die gleichen Ergebnisse wie ich auch. TRUE bedeutet, dass wir jedes mal erneut die Chance haben, Zahlen von 1 bis 6 zu würfen.

Wie oft sollte jede Zahl vorkommen? Genau 1667 mal in etwa. Da jede Zahl die gleiche Chance hat, aufzutreten, teilen wir 10000 durch 6 und erhalten 1667 (gerundet). Ähnliche Werte bekommen wir von R zurück:

table(my_sample)
my_sample
   1    2    3    4    5    6 
1657 1712 1667 1673 1656 1635 

Wir sehen, dass nicht jede Zahl exakt 1667 mal gewürfelt wurde, aber diese statistischen Schwankungen sind zu erwarten. Würden wir dieses Experiment nicht 10000, sondern eine Millionen mal durchführen, würden wir fast exakt für jeden Würfel die gleiche Wahrscheinlichkeit erhalten. Reale Würfel hingegen sind nicht perfekt. Schon eine abgebrochene Kante könnte dazu führen, dass manche Zahlen wahrscheinlicher gewürfelt werden.

Ein manipulierter Würfel

Was passiert nun, wenn der Würfel manipuliert wurde? Wie können wir feststellen, ob die Unterschiede im Auftreten der Zahlen nicht zufälligen statistischen Schwankungen (wie beim nicht-manipulierten Würfen) unterliegen, sondern systematisch sind? Um diese Frage zu beantworten wiederholen wir das Experiment 1000 mal. 1000 mal Würfeln wir den Würfel 100 mal und schauen, wie oft die Zahl 4 darin vorkommt. Statistisch sollten wir annehmen, dass die Zahl vier am häufigsten in etwa 1667 mal auftritt.

sample_distribution <- 1:1000 %>% 
  map_dbl(~ sample(1:6, 100, TRUE) %>% 
            table(.) %>% .[names(.) == 4] %>% 
            as.vector) %>% table

Diese Funktion ist komplizierter. Stück für Stück:

  • 1:1000: Wir erstellen einen Vector mit den Zahlen 1 bis 1000
  • map_dbl: Für jedes dieser Zahlen von 1 bis 1000 lassen wir eine Funktion drüber laufen
  • sample(1:6, 1000, TRUE): Wir würfen den Würfel 100 mal
  • table(.): Wir lassen uns die Häufigkeiten der Zahlen bei diesen 100 Würfen anzeigen
  • .[names(.) == 4]: Wir zählen, wie oft die Zahl 4 in diesen Häufigkeiten vor kommt
  • table: Aus den 1000 Experimenten zählen wir, wie oft die 4 in jedem der Experimente vorgekommen ist.

Wir können nun diese Verteilung visuell darstellen:

ggplot(as.data.frame(sample_distribution), aes(x = . , y = Freq)) + 
  geom_bar(stat = "identity") +
  theme_bw()

Sampling distribution

Ok, anscheinend ist es sehr wahrscheinlich bei 100 Würfen die Zahl 4 16 oder 17 mal zu erhalten. Die Zahl 4 aber 27 oder gar 6 mal zu erhalten ist schon äußerst unwahrscheinlich, wenn wir annehmen, dass der Würfel nicht manipuliert wurde.

Was du hier siehst ist eine Wahrscheinlichkeitsverteilung. Anhand der Verteilung können wir bestimmen, wie wahrscheinlich ein Ereignis (das Auftreten der Zahl 4 bei 100 Würfen) auftritt. Auf Grundlage der von uns simulierten Grafik können wir sagen, dass ein Würfel vermutlich manipuliert wurde, wenn die 4 27 mal auftritt.

Binomialverteilung

Die eben erstellte Wahrscheinlicheitsverteilung können wir heran ziehen, um zu bestimmen, ob der Würfel manipuliert wurde. Würden wir beispielsweise die Zahl 4 nur 6 mal bei 100 Würfen erhalten, wäre es sehr unwahrscheinlich, dass der Würfel manipuliert wurde. Und genau so testen wir in der Regel wissenschaftliche Fragestellungen. Wir fragen uns nicht, ob unsere wissenschaftliche Frage korrekt ist, sondern ob unser Ereignis (die Unterschiede der Gewichtsreduzierung der Diäten) unwahrscheinlich ist, wenn wir keine Unterschiede annehmen.

Diese Wahrscheinlichkeitsverteilungen müssen wir aber nicht jedes Mal neu simulieren. Sie wurden bereits berechnet. Für das Würfelbeispiel ist es sinnvoll, eine Binomialverteilung zur verwenden. Nichts anderes ist im Prinzip unsere Wahrscheinlichkeitsverteilung. Mit einer Binomialverteilung testen wir die Wahrscheinlichkeit für n Ereignisse mit einer bestimmten Wahrscheinlichkeit. Zum Beispiel: Wie wahrscheinlich ist es, dass wir bei 5 Würfen 3 mal die Zahl 4 würfeln?

dbinom(3, 5, 1/6)
# 0.03215021 -> bei etwa 3%

Unsere Wahrscheinlichkeitsverteilung können wir nun auch als Binomialverteilung darstellen:

dist <- tibble(x = 4:32, y = dbinom(4:32, 100, 1/6))

ggplot(dist, aes(x, y)) + 
  geom_bar(stat = "identity")

Ab wann glauben wir also, dass der Würfel vermutlich kein normaler Würfel ist (sondern manipuliert wurde)? Vorhin hatten wir einen Cut-Off willkürlich bestimmt. Wissenschaftler machen es ganz ähnlich.

Wenn die Wahrscheinlichkeit für ein Ereignis unter 5% auf Grundlage einer Wahrscheinlichkeitsverteilung ist, ist ein Ereignis unwahrscheinlich.

In unserer Verteilung sind dies folgende Ereignisse:

Wenn wir also die Zahl 4 22 mal oder öfter würfeln, müssen wir davon ausgehen, dass es kein normaler Würfel ist, sondern, dass der Würfel vermutlich manipuliert wurde. Wir sprechen dann von einem signifikanten Ereignis.

Die F-Verteilung

Mit der Binomialverteilung testen wir, ob das Auftreten einer bestimmten Anzahl von Ereignissen wahrscheinlich oder unwahrscheinlich ist. Mit der F-Verteilung testen wir, ob die Unterschiede zweier Varianzen statistischen Schwankungen unterliegen oder systematisch sind. Um das zu verstehen, müssen wir erstmal verstehen, was mit Varianz gemeint ist.

Varianz

Varianz ist ein statistisches Maß, mit dem wir anzeigen können, wie stark eine Variable variiert. Menschen beispielsweise variieren in ihrem Gewicht. Manche Menschen sind schwer, andere leicht. Berechnet wir die Formel folgendermaßen:

σ2=i=1n(xix)2n1\sigma^2 = \frac{\displaystyle\sum_{i=1}^{n}(x_i - \overline{x})^2}{n - 1}
  • x\overline{x}: Der Mittelwert der Variable: mean(diet_dataset$weight6weeks, na.rm = TRUE)
  • xix_i: Jeder einzelne Gewichtswert
  • n1n - 1: Die Anzahl der Menschen, deren Variable Gewicht minus 1.

In anderen Worten gesprochen, summieren wir die quadrierten Abweichungen der Einzelgewichtswerte mit dem Mittelwert der Stichprobe und teilen diese Summe durch n - 1. Es ist immer ganz praktisch, sich solche mathematische Ideen auch grafisch vorzustellen.

sample <- diet_dataset[c(1:5), ]

ggplot(sample, aes(Person, weight6weeks)) +
  geom_point() +
  geom_hline(yintercept = mean(sample$weight6weeks)) +
  geom_rect(aes(xmin = Person, 
                xmax = Person + (abs(weight6weeks - mean(weight6weeks))),   
                ymin = weight6weeks, 
                ymax = mean(weight6weeks), 
                fill = "red",
                alpha = .1), 
            data = sample) +
  guides(fill = FALSE, alpha = FALSE) +
  theme_bw() +
  coord_fixed()

Die horizontale Linie stellt den Mittelwert unser Variable sample dar. Jeder Punkt ist das Gewicht einer jeden Person nach 6 Wochen.

Varianz ist in diesem Beispiel nichts anderes als die durchschnittliche Fläche dieser Quadrate geteilt durch die Anzahl der Quadrate - 1. Wir haben bisher offen gelassen, weshalb wir eigentlich diese Quadrate nicht einfach durch die Anzahl der Quadrate teilen, sondern durch die Anzahl der Quadrate minus 1. Dies liegt daran, dass wir die Varianz einer Population in der Regel unterschätzen, wenn wir nur ein paar Personen aus einer Stichprobe ziehen (hier die Anzahl der Personen, die an unserem Experiment teilnehmen). Man nennt diese Korrektur auch Bessel’s Correction.

Varianz für unsere Variable weight6weeks sieht grafisch also folgendermaßen aus:

ggplot(diet_dataset, aes(Person, weight6weeks)) +
  geom_point() +
  geom_hline(yintercept = mean(diet_dataset$weight6weeks)) +
  geom_rect(aes(xmin = Person, 
                xmax = Person + (abs(weight6weeks - mean(weight6weeks))),   
                ymin = weight6weeks, 
                ymax = mean(weight6weeks), 
                fill = Person,
                alpha = .1), 
            data = diet_dataset) +
  guides(fill = FALSE, alpha = FALSE) +
  theme_bw() +
  coord_fixed()

Erneut müssen wir lediglich die Summe der Fläche dieser Quadrate durch die Anzahl der Quadrate minus 1 teilen und erhalten die Varianz.

Je kleiner die Fläche dieser Quadrate ist, desto geringer ist die Varianz der Variable. In anderen Worten, je geringer die Einzelwerte vom Mittelwert der Variable abweichen, desto geringer ist die Varianz.

In R können wir diese Varianz folgendermaßen berechnen:

((diet_dataset$weight6weeks - mean(diet_dataset$weight6weeks))**2 %>% sum) / (nrow(diet_dataset) - 1)

oder

var(diet_dataset$weight6weeks)

Beide Befehle ergeben eine Varianz von 79.64677.

Varianz-Quotienten

Wir können Varianzen miteinander vergleichen, indem wir einen Quotienten zweier Varianzen bilden. Beispielsweise könnten wir die Varianz der ersten Diätgruppe mit der Varianz der zweiten Diätgruppe vergleichen.

diet_dataset %>%
  group_by(Diet) %>%
  summarise(
    variance = var(weight_diff)
  )
# A tibble: 3 x 2
   Diet variance
  <int>    <dbl>
1     1     5.02
2     2     6.37
3     3     5.74
5.02 / 6.37
# [1] 0.7880691

Bei einem Wert über 1 wäre die Varianz der ersten Gruppe größer als die Varianz der zweiten Gruppe. Bei einem Wert unter 1 ist die Varianz der ersten Gruppe kleiner als die Varianz der zweiten Gruppe.

Die F-Wert ergibt sich aus einem solchen Quotienten. Aber anstatt willkürlich Varianzen mit einem Quotienten zu vergleichen, vergleichen wir in der Regel eine systematische gegen eine nicht-systematische Varianz. Eine systematische Varianz tritt in der Regel durch Manipulationen auf, die wir selber herbei führen. Beispielsweise haben wir die 60 Menschen 3 Diäten zugeordnet. Diese Zuordnung war systematisch. Nicht-systematische Varianz ist Varianz, die durch die Daten gegeben ist. Beispielsweise die Varianz aller Daten um den Mittelwert der Variable.

F-Wert

Der F-Wert ist der Quotient der Varianz zwischen unseren 3 Diätgruppen (SSB) und der unsystematischen Varianz in unseren Daten (SSW).

Sum of Squares Between (SSB)

Die Varianz zwischen den Gruppen berechnen wir, indem wir den Mittelwert der Gruppen vom Gesamtmittelwert der Variable Gewichtsreduzierung abziehen, quadrieren und mal der Anzahl der Fälle berechnen.

diet_dataset %>%
  group_by(Diet) %>%
  summarise(
    variance = var(weight_diff),
    mean = mean(weight_diff),
    n = n()
  )

ss_between <- tibble(group = c(rep(-3.3, 24), rep(-3.03, 27), 
                               rep(-5.15, 27)),
                     person = seq(1:nrow(diet_dataset)),
                     diet = as.factor(c(rep("Diet One", 24), 
                                        rep("Diet Two", 27), 
                                        rep("Diet Three", 27))))

ggplot(ss_between, aes(person, group)) +
  geom_point() +
  geom_hline(yintercept = mean(diet_dataset$weight_diff)) +
  geom_rect(aes(xmin = person,
                xmax = person + abs(group - mean(diet_dataset$weight_diff)),
                ymin = group,
                ymax = mean(diet_dataset$weight_diff),
                fill = diet,
                alpha = .01,),
            data = ss_between) +
  guides(fill = FALSE, alpha = FALSE) +
  theme_bw() +
  coord_fixed()

Die Quadrate sind nur schlecht zu erkennen, da die Skalierung so unterschiedlich ist. Der Strich in der Mitte kennzeichnet den Mittelwert der Variable, die Punkte sind die Mittelwerte der Einzelgruppen. Es gibt so viele Punkte, wie es Menschen in jeder Gruppe gibt. Wenn wir die Quadrate zerren, sieht dies so aus:

Diese Quadrate müssen wir erneut aufsummieren.

Sum of Squares Within (SSW)

Die Varianz innerhalb der Gruppe ist nichts anderes als die Summe der Varianz innerhalb der Gruppen.

diet_dataset <- diet_dataset %>%
  arrange(Diet) %>%
  mutate(
    mean_group = c(rep(-3.3, 24), rep(-3.03, 27), rep(-5.15, 27))
  )

ggplot(diet_dataset, aes(Person, weight_diff)) +
  geom_point() +
  geom_segment(aes(x = 0, xend = 24, y = -3.3, yend = -3.3)) +
  geom_segment(aes(x = 25, xend = 51, y = -3.03, yend = -3.03)) +
  geom_segment(aes(x = 52, xend = 78, y = -5.15, yend = -5.15)) +
  geom_rect(aes(xmin = Person,
                xmax = Person + abs(weight_diff - mean_group),
                ymin = weight_diff,
                ymax = mean_group,
                fill = as.factor(Diet),
                alpha = .01,),
            data = diet_dataset) +
  theme_bw() +
  guides(fill = FALSE, alpha = FALSE) +
  coord_fixed()    

F-Wert

Wir haben bisher nur die Quadrate aufaddiert, aber noch keine Varianz berechnet. Dazu müssen wir SSB und SSW durch einen Nenner teilen:

  • Mean Squares Between (MSB): SSB / (k - 1). K sind die Anzahl der Gruppen, hier 2.
  • Mean Squares Within (MSW): SSW / (n - k). N ist die Anzahl der Personen in allen Gruppen, hier 78

Der F-Wert bildet sich nun aus dem Quotienten zwischen MSB und MSW.