81 Der Pre-Test oder Vortest

Letzte Änderung am 07. May 2025 um 20:37:43

“I struggle with some demons; They were middle class and tame.” — Leonard Cohen, You Want It Darker

In diesem Kapitel soll es um den Pre-Test oder auch Vortest gehen. Wir sind hier in einem experimentellen Design, welches verschiedene Gruppen beinhaltet oder aber wir wollen wissen, ob wir in einer linearen Regression die Normalverteilung unseres Messwertes $y$ vorliegen haben. Grundsätzlich geht es erstmal darum herauszufinden, ob die Annahmen an einen statistischen Test in deinen Daten erfüllt sind. Häufig wollen wir eine ANOVA für einen Gruppenvergleich rechnen und dann anschließend einen multiplen Test oder Post-hoc Test durchführen. In beiden Fällen wird es einfacher wenn wir eine Normalverteilung in unseren Messwert $y$ sowie eine Varianzhomogenität in unseren Behandlungsgruppen oder Faktoren $f$ vorliegen haben. Mit einfacher meine ich, dass du auch mit einer Abweichung von der Normalvertielung und auch Varianzhterogenität heutzutage umgehen kannst. Der Standard im statistischen Testen war aber immer die Normalverteilung und die Varianzhomogenität. Wenn beides nicht vorlag, dann wurde es manchmal etwas dunkel. Wir aber aber im 21. Jahrhundert ein paat mehr Pfeile im Köcher und können mit unterschiedlichsten Daten und deren Eigenschaften umgehen. Mehr dazu findest du im Teil zum statistsichen Modellieren und den nachfolgenden Kapiteln.

81.1 Allgemeiner Hintergrund

Wir werden uns in diesem Kapitel auf das faktorielle Experiment konzentrieren. Natürlich kannst du auch alle Funktionen in einem anderen Design anwenden. Wenn du wissen willst, ob eine Variable normalverteilt ist oder aber ein Gruppenfaktor homogen in den Varianzen, dann helfen dir hier auch die Funktionen weiter. Häufig werden aber die beiden Eigenschaften Normalverteilung und Varianzhomogenität in Gruppenvergleichen verwendet.

Eine Sache ist aber wichtig zu wissen. Wir untersuchen in unseren Experimenten ja immer nur eine Stichprobe der Grundgesamtheit und wollen dann von der Stichprobe einen Rückschluß auf die Grundgesamt machen. Wenn dich mehr dazu interessiert, dann schaue einmal in dem Kapitel zum Testen von Hypothesen rein. Es kann also sein, dass wir definitiv in der Grundgesamtheit einen normalverteilten Messwert vorliegen haben, wir aber noch zu wenige Beobachtungen in unsere Stichprobe vorliegen haben um diese Normalverteilung in einem Histogramm oder Densityplot zu sehen. Nehmen wir einmal die Körpergröße als ein normalverteilten Messwert $y$ an. Wir wissen, dass die Körpergröße einer Normalverteilung folgt. In der folgenden Abbildung Abbildung 81.1 siehst du einmal die Körpergrößen von unseren Gummibärchendaten. Insgesamt haben $432$ Männer und $418$ Frauen bei der Gummibärchenumfrage mitgemacht. Dennoch beobachten wir keine saubere Normalverteilung, wie wir sie erwarten würden. Wir haben noch zu wenige Beobachtungen gemacht.

Abbildung 81.1— Darstellung der Körpergröße in [cm] für die Geschlechter getrennt. Die Körpergröße ist normalverteilt. Die Farben repräsentieren die jeweiligen Geschlechter. Die Männer sind blau und die Frauen in lila dargestellt. **(A)** Histogramm. **(B)** Densityplot. *[Zum Vergrößern anklicken]*

Wir sehen also, nur weil etwas wie die Körpergröße wirklich normalverteilt ist, ist es noch etwas ganz anders diese Normalverteilung dann auch in den Messwert $y$ zu beobachten. Die Fallzahlen in der Grundgesamt und in der Stichprobe unterscheiden sich dann doch gewaltig und wir sind dann eben auch auf Annahmen angewiesen. Meistens passt es auch mit den Annahmen und wenn wir mal daneben liegen, kann es sein, dass es dann doch nicht so viel ausmacht, wenn der Effekt in unserer statistischen Auswertung groß genug ist.

Das Modell

Auch hier möchte ich einmal das statistische Modell besprechen was wir in dem Gruppenvergleich oder dem statistischen Modellieren benötigen. Im Folgenden findest du einmal ein faktorielles Modell mit einem Messwert $y$ und zwei Gruppenfaktoren $f_A$ und $f_B$. Die beiden Faktoren entsprechen zwei unterschiedlichen kategoriellen Variablen mit verschiedenen Gruppen. Wir wollen uns ja in diesem Kapitel auf die Normalverteilung und die Varianzhomogenität konzentrieren. Die beiden Gütekriterien können aber ganz klar dem Messwert $y$ und den experimentellen Faktoren zugeordnet werden.

\[ \underbrace{\;\mbox{Messwert}\; y\;}_{normalverteilt} \sim \overbrace{\;\mbox{Faktor}\; f_A + \mbox{Faktor}\; f_B\;}^{homogene\; Varianzen} \]

mit

$\mbox{Messwert}\; y$, gleich dem Messwert oder Outcome, wie die Sprungweite in [cm] als jump_length in unseren Beispieldaten.
$\mbox{Faktor}\; f_A$, gleich dem ersten Faktor $f_A$, wie die Tierart als animal mit unterschiedlichen Gruppen oder Leveln.
$\mbox{Faktor}\; f_B$, gleich dem zweiten Faktor $f_B$, wie der Messort als site mit unterschiedlichen Gruppen oder Leveln.

Damit haben wir uns erstmal für die Vortest für die Normalverteilung und die Varianzhomogenität geordnet. Wir wollen dann in den folgenden Abschnitten noch andere Gütekriterien eines Modells kurz anreißen, aber den Hauptteil findest du im Kapitel zur Modelgüte von linearen Modellen.

Gibt es noch mehr Vortests?

Jetzt könnte man meinen, dass mit der Normalverteilung und der Varianzhomogenität eigentlich die wichtigsten Gütekriterien vorgetestet werden. Es gibt aber für lineare Modelle, was ein Gruppenvergleich dann am Ende auch nur ist, noch andere Gütekriterien. Neben diesen beiden Eigenschaften können wir usn auch noch folgende weitere Gütekriterien anschauen. Ich verweise hier einmal auf die Hilfeseite des R Packetes {performance} für mehr Informationen und deren Referenzseite der Familie der check_*() Funktionen. Wie immer kommt es dann auf die Fragestellung und dann auf das enstprechende Modell sowie den verwendeten Algorithmus an. Je nachdem was du gemessen hast, also welche Werte dein $y$ annimmt, musst du einen anderes Modell wählen. Je nach Modell hast du dann auch andere Annahmen. Das würde hier aber das Kapitel sprengen. Gerne kannst du als Startpunkt einmal in das Teil zum statistsichen Modellieren reinschauen.

Betrachten wir also einmal im Folgenden die beiden wichtigsten Annahmen an ein faktorielles Design oder aber Gruppenvergleich. Wir fragen uns, haben wir eine Normalverteilung in den Messwerten $y$ und homogene Varianzen in den Faktoren oder Gruppen $f$ vorliegen? Dann können wir ganz normal eine ANOVA oder einen Tukey HSD Test rechnen.

81.1.1 Normalverteilung

Fangen wir also mit der Annahme der Normalverteilung an die Daten an. Hierbei ist wichtig, dass wir nicht die Daten insgesamt betrachten sondern uns fragen, ob der betrachtete Messwert $y$ im Modell oder statistischen Test normalverteilt ist. Wir haben uns den Zusammenhang ja schon oben einmal in dem statistischen Modell angeschaut. Häufig führt dies zu Verwirrungen, da verallgemeinert von den Daten gesprochen wird, die normalverteilt sein soll. Hier geht es dann wirklich nur um deinen Messwert $y$. Das nochmal als Erinnerung für die weiteren Betrachtungen. Was wären also beispielhaft normalverteilte Messwerte?

Beispiel: Frischgewicht, Trockengewicht, Chlorophyllgehalt, Pflanzenhöhe

Tabelle 81.1— Tabelle mit beispielhaften, normalverteilten Messwerten $y$.

freshmatter	drymatter	chlorophyll	height
8.23	1.21	45.88	24.19
2.61	0.87	43.91	18.51
4.81	0.34	37.44	21.74

Nach dem zentralen Grenzwertsatz können wir bei Merkmalen, die sich aus verschiedenen Einflussfaktoren zusammensetzen, allgemein von einer Normalverteilung ausgehen. Die Körpergröße oder das Körpergewicht ist normalverteilt, da wir hier es mit vielen Einflussgrößen zu tun haben, die das tatsächliche Körpergewicht einer Beobachtung ausmachen. Das Körpergewicht hängt eben von der täglichen Kalorienmenge, verschiedensten Genen, dem Muskelanteil, dem Aktivitätsgrad, der sozialen Stellung und vielen weiteren Einflusfaktoren ab. Alles zusammen addiert sich dann zum Körpergewicht wobei jeder Einflussfaktor nur einen kleinen Teil ausmacht.

Was heißt approximativ normalverteilt?: Wir sprechen von approximativ normalverteilt, wenn wir meinen, dass ein Messwert $y$ in unserer Stichprobe annähernd normalverteilt ist. Wir sind uns also nicht zu hundertprozent sicher, glauben aber das die Normalverteilungsannahme an unseren Messwert passen wird. Häufig sagen wir auch, dass gewisse Tranformationen approximativ normalverteilt sind. So haben wir nach einer log-Transformation log-normalverteilte Messwerte vorliegen. Wir sagen dann meistens, dass ein log-transformierter Messwert approximativ normalverteilt ist.

81.1.2 Varianzhomogenität

Kommen wir nun zur Varianzhomogenität oder Varianzheterogenität in den Gruppen des Behandlunsgfaktors. Je nachdem was du betrachtest, nennen wir es eben Varianzhomogenität oder Varianzheterogenität. Entweder sind die Varianzen gleich, dann haben wir Varianzhomogenität vorliegen oder die Varianzen in den Gruppen sind nicht gleich, dann hast du Varianzheterogenität in den Daten. Es gibt so ein paar Daumenregeln, die dir helfen abzuschätzen, ob in deinen Gruppen Varianzheterogenität vorliegt. Um es kurz zu machen, vermutlich hast du mindestens leichte Varianzheterogenität in den Daten vorliegen. Es ist bei kleinen Gruppengrößen nicht zu vermeiden, dass sich die Varianzen eben dann doch unterscheiden.

Du hat viele Behandlungsgruppen. Je mehr Gruppen du hast oder eben dann auch Faktorkombinationen, die du testen möchtest, desto wahrscheinlicher wird es, dass mindestens eine Gruppe eine unterschiedliche Varianz hat. Du hast Varianzheterogenität vorliegen.
Du misst deine Gruppen über die Zeit. Je größer, schwerer oder allgemein höher ein Messwert wird, desto größer wird auch die Varianz. Schaust du dir deine Messwerte über die Zeit an hast du meistens Varianzheterogenität vorliegen.
Deine Kontrolle ohne Behandlung verhält sich meistens nicht so, wie die Gruppen, die eine Behandlung erhalten haben. Wenn du nichts machst in deiner negativen Kontrolle, dann hast du meistens eine andere Streuung der Messwerte als unter einer Behandlung.
Deine Behandlungen sind stark unterschiedlich. Wenn deine Behandlungen sich biologisch oder chemisch in der Wirkung unterscheiden, denn werden vermutlich deine Messwerte auch anders streuen. Hier spielt auch die Anwendung der Behandlung und deren Bereitstellung eine Rolle. Wenn was nicht gleich ist, dann wird es vermutlich nicht gleiche Messwerte erzeugen.
Du hast wenig Fallzahl pro Gruppe oder Faktorkombination. Wenn du wenig Fallzahl in einer Gruppe hast, dann reicht schon eine (zufällige) Messabweichung und schon sind deine Varianzen heterogen.

In den folgenden Beispielen habe ich dir einige Experimente mit einem faktoriellen Design mitgebracht. Die Fotos stammen aus wissenschaftlichen Publikationen wie einer wissenschaftlichen Veröffentlichung oder aber wissenschaftlichen Postern hier auf dem Gelände der Hochschule Osnabrück.

Abbildung 81.2— Ein zweifaktorielles Experiment mit neuen Faktorkombinationen die alle miteinander paarweise vergleichen werden. Wir sehen sehr gut, dass die Kontrolle sehr kleine Werte hat und somit die Varianz in der Kontrolle sehr viel kleiner ist alles in den anderen Gruppen.

Abbildung 81.3— Ein zweifaktorielles Experiment mit sechs Faktorkombinationen. Wir sehen hier sehr gut, dass mit steigenden MIttelwerten, also höheren Barplots, auch die Varianz in den Gruppen zunimmt. Die Fehlerbalken werden immer länger.

Abbildung 81.4— Ein zweifaktorielles Experiment mit sehr vielen Faktorrkombinationen. Durch die unterschiedlichen miitelren Zählwerte ergeben sich sehr viele unterschiedlich große Mittelwerte. Darüber hinaus haben wir sehr viele Gruppen. Wir sehen hier sehr viel Varianzhterogenität.

Abbildung 81.5— Zwei einfaktorielle Experiemente in einer Abbildung dargestellt. Die linken Barplots und die rechten Barplots wurden getrennt voneinander ausgewertet. Auch hier sieht man sehr viel unterschiedliche Streuung in den Daten.

Jetzt haben wir uns einmal die wichtigsten Abbildungen angeschaut und haben so eine erste Idee was Varianzhomogenität sein könnte. Wir schauen uns dann in den folgenden Abschnitten noch mehr zu der Bestimmung an. Dann bleibt eigentlich noch eine abschließende Frage für den einführende Abschnitt.

Tut Varianzheterogenität anstatt Varianzhomogenität weh?: Nein. Meistens ist die Varianzheterogenität nicht so ausgeprägt, dass du nicht auch eine ANOVA oder anderen statistischen Test rechnen kannst. Über alle Gruppen hinweg wird dann zwar zum Beispiel in einer ANOVA die Varianz gemittelt und es kann dann zu weniger signifikanten Ergebnissen führen, aber so schlimm ist es nicht. Im Post-hoc Test solltest du aber die Varianzheterogenität berücksichtigen, da du ja immer nur zwei Gruppen gleichzeitig betrachtest. Aber auch hier gibt es dann die passenden Adjustierungen. Mehr dazu am Ende des Kapitels im Abschnitt zu den Auswegen.

81.2 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

R Code [zeigen / verbergen]

pacman::p_load(tidyverse, magrittr, olsrr,
               broom, car, performance, 
               see, scales, readxl, nlme,
               parameters, conflicted)
conflicts_prefer(dplyr::filter)
conflicts_prefer(dplyr::select)

An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.

81.3 Daten

Wir immer bringe ich hier ein paar Datensätze mit damit wir dann verstehen, was eigentlich in den folgenden Analysen in R und den entsprechenden R Paketen passiert. Ich zeige hier an den Daten nur die Anwendung in R. Deshalb fehlen dann hier auch die Mittelwerte und andere deskriptive Maßzahlen. Schauen wir jetzt also mal in unsere Beispieldaten für die einfaktorielle und zweifaktorielle Datenanalyse rein.

Einfaktorieller Datensatz

Beginnen wir mit einem einfaktoriellen Datensatz. Wir haben hier als Messwert die Sprungweite von Flöhen in [cm] vorliegen. Wissen wollen wir, ob sich die Sprungweite für drei verschiedene Floharten unterscheidet. Damit ist dann in unserem Modell der Faktor animal und die Sprungweite jump_length als Messwert. Ich lade einmal die Daten in das Objekt fac1_tbl. Hier haben wir dann ein simples Design vorliegen.

R Code [zeigen / verbergen]

fac1_tbl <- read_xlsx("data/flea_dog_cat_fox.xlsx") |>
  select(animal, jump_length) |> 
  mutate(animal = as_factor(animal))

Dann schauen wir uns die Daten einmal in der folgenden Tabelle als Auszug einmal an. Wichtig ist hier nochmal, dass du eben einen Faktor animal mit drei Leveln also Gruppen vorliegen hast. Wir wollen jetzt die drei Tierarten hinsichtlich ihrer Sprungweite in [cm] miteinander vergleichen.

Tabelle 81.2— Tabelle der Sprungweiten in [cm] als Messwert $y$ von Hunde-, Katzen- und Fuchsflöhen. Der Datensatz ist einfaktoriell, da wir nur einen Faktor vorliegen haben.

tinytable_bv8b72awxbre0rw1mobs

animal	jump_length
dog	5.7
dog	8.9
dog	11.8
...	...
fox	10.6
fox	8.6
fox	10.3

Und dann wollen wir uns noch einmal die Daten als einen einfachen Boxplot anschauen. Wir sehen, dass die Daten so gebaut sind, dass wir einen signifikanten Unterschied zwischend den Sprungweiten der Floharten erwarten. Die Boxen der Boxplots überlappen sich nicht und die Boxplots liegen auch nicht auf einer Ebene. Wir könnten hier von normalverteilten Daten und Varianzhomogenität ausgehen. Die Mediane liegen in der Mitte der Boxen und die Boxen sind ungefähr gleich groß.

Abbildung 81.6— Beispielhafter einfaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten.

Zweifaktorieller Datensatz

Neben dem einfaktoriellen Datensatz wollen wir uns noch den häufigeren Fall mit zwei Faktoren anschauen. Wir haben also nicht nur die drei Floharten vorliegen und wollen wissen ob diese unterschiedlich weit springen. Darüber hinaus haben wir noch einen zweiten Faktor gewählt. Wir haben die Sprungweiten der Hunde-, Katzen- und Fuchsflöhe nämlich für die beiden Geschlechter getrennt gemessen. Dadurch haben wir jetzt den Faktor animal und den Faktor sex vorliegen. Wiederum fragen wir uns, ob sich die Sprungweite in [cm] der drei Floharten in den beiden Geschlechtern unterscheidet. Darüber hinaus haben wir neben der Sprungweite noch die Schlupfzeiten in [m] gemessen. Im Folgenden lade ich einmal den Datensatz in das Objekt fac2_tbl.

R Code [zeigen / verbergen]

fac2_tbl <- read_xlsx("data/flea_dog_cat_length_weight.xlsx") |> 
  select(animal, sex, jump_length, hatch_time) |> 
  mutate(animal = as_factor(animal),
         sex = as_factor(sex))

Betrachten wir als erstes einen Auszug aus der Datentabelle. Wir haben hier als Messwert oder Outcome $y$ die Sprungweite jump_length sowie die Schlupfzeiten hatch_time vorliegen. Als ersten Faktor die Variable animal und als zweiten Faktor die Variable sex festgelegt.

Tabelle 81.3— Tabelle der Sprungweiten in [cm] und Schlupfzeiten [m] als Messwert $y$ von Hunde-, Katzen- und Fuchsflöhen der beiden Geschlechter. Der Datensatz ist zweifaktoriell, da wir einen Behandlungsfaktor mit animal und einen zweiten Faktor mit sex vorliegen haben.

tinytable_hpvtfkb5hhaa62ra1k4i

animal	sex	jump_length	hatch_time
cat	male	15.79	483.6
cat	male	18.33	82.56
cat	male	17.58	296.73
...	...	...	...
fox	female	24.35	182.68
fox	female	24.36	104.89
fox	female	22.13	62.99

Auch hier schauen wir uns einmal die Daten in einem Boxplot und einem Densityplot an. Wir wollen ja sehen, ob sich zum einen die Gruppen unterscheiden und zum anderen wie unsere Messwerte der Sprungweiten und der Schlupfzeiten verteilt sind. Wir erkennen in den Boxplots und auch in den Densityplots, dass wir vermutlich eine approximative Normalverteilung in den Sprungweiten vorliegen haben, aber auf keinen Fall eine Normalverteilung in den Schlupfzeiten. Du siehst hier nochmal in den beiden Abbildungen die Schiefe in der Verteilung der Schlupfzeiten.

Boxplot
Densityplot

Abbildung 81.7— Zweifaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Geschlechtern.

Abbildung 81.8— Zweifaktorieller Boxplot für die Schlupfzeiten in [m] gruppiert nach den Floharten und den beiden Geschlechtern.

Abbildung 81.9— Densityplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Geschlechtern.

Abbildung 81.10— Densityplot für die Schlupfzeiten in [m] gruppiert nach den Floharten und den beiden Geschlechtern.

81.4 Visuelle Überprüfung

text

Hier einmal die theoretischen Abbilungen von drei Gruppen mit Varianzhomogenität und Varianzheterogenität. Hier wird schnell klar, dass du mit mehr Behandlungsgruppen kaum die Annahme für Varianzhomogenität aufrecht erhalten kannst. Insbesondere Kontrollen oder Standardbehandlungen haben große oder kleinere Varianzen als der Rest der Behandlungen. Manchaml hat man aber Glück und die Varianzen sind eher homogen.

Abbildung 81.11— Darstellung der Varianzhomogenität und Varianzheterogenität in einem Barplot. *[Zum Vergrößern anklicken]*

Abbildung 81.12— Darstellung der Varianzhomogenität und Varianzheterogenität in einem Boxplot. *[Zum Vergrößern anklicken]*

Wir bauen wir usn jetzt die Barplots oder Säulendiagramme? Entweder direkt aus {emmeans}, wie ich dir das weiter unten zum Compact letter display zeige oder eben als Wiederholung des Kapitels zur Visualisierung von Daten. Wir bauen uns also einmal den einfaktoriellen Barplot wie auch den zweifaktoriellen Barplot.

81.4.1 Einfaktorieller Barplot

Für den einfaktoriellen Barplot brauchen wir einmal den Mittelwert und die Standardabweichung der einzelne Floharten.

R Code [zeigen / verbergen]

stat_fac1_tbl <- fac1_tbl |> 
  group_by(animal) |> 
  summarise(mean = mean(jump_length),
            sd = sd(jump_length))
stat_fac1_tbl

# A tibble: 3 × 3
  animal  mean    sd
  <fct>  <dbl> <dbl>
1 dog     8.13  2.14
2 cat     4.74  1.90
3 fox     9.16  1.10

Dann können wir uns auch schon den Barplot mit {ggplot} erstellen. Später willst du dann noch das Compact letter display über den Balken als Mittelwerte ergänzen.

R Code [zeigen / verbergen]

ggplot(data = stat_fac1_tbl, 
       aes(x = animal, y = mean, fill = animal)) +
  theme_minimal() +
  geom_bar(stat = "identity") + 
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), 
                width = 0.2) + 
  labs(x = "Flohart", y = "Sprungweite in [cm]") +
  theme(legend.position = "none") + 
  scale_fill_okabeito()

Abbildung 81.13— Beispielhafter einfaktorieller Barplot für die Sprungweiten in [cm] gruppiert nach den Floharten.

81.4.2 Zweifaktorieller Barplot

Für den zweifaktoriellen Barplot brauchen wir einmal den Mittelwert und die Standardabweichung der einzelne Floharten getrennt für die beiden Messorte. Das geht natürlich auch umgedreht, also die Messorte für die Floharten. Das kommt dann auf deine Fragestellung an.

R Code [zeigen / verbergen]

stat_fac2_tbl <- fac2_tbl |> 
  group_by(animal, sex) |> 
  summarise(mean = mean(jump_length),
            sd = sd(jump_length))
stat_fac2_tbl

# A tibble: 6 × 4
# Groups:   animal [3]
  animal sex     mean    sd
  <fct>  <fct>  <dbl> <dbl>
1 cat    male    15.4  1.85
2 cat    female  20.5  1.78
3 dog    male    18.1  1.90
4 dog    female  22.9  2.06
5 fox    male    20.7  1.95
6 fox    female  25.5  2.06

Und dann können wir auch schon den zweifaktoriellen Barplot in {ggplot} erstellen. Du musst schauen, was du auf die x-Achse legst und was du dann auf die Legende und daher mit fill gruppierst. Damit die Positionen passen, spiele ich hier noch mit der Funktion position_dodge() rum.

R Code [zeigen / verbergen]

ggplot(data = stat_fac2_tbl, 
       aes(x = sex, y = mean, fill = animal)) +
  theme_minimal() +
  geom_bar(stat = "identity", width = 0.9, 
           position = position_dodge(0.9)) + 
  geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), 
                width = 0.2, 
                position = position_dodge(0.9)) + 
  labs(x = "Flohart", y = "Sprungweite in [cm]", fill = "Tierart") +
  scale_fill_okabeito()

Abbildung 81.14— Beispielhafter zweifaktorieller Barplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Messorten.

81.4.3 Einfaktorieller Boxplot

Das praktische bei den Boxplots ist, dass wir hier nichts mehr vorrechnen müssen, sondern direkt die Boxplots in {ggplot} erstellen können. Ich finde man sieht immer in einem Boxplot besser, ob die Streuung um den Median eher homogen oder eher heterogen ist. Gerne ergänze ich noch den Mittelwert mit der Funktion stat_summary().

R Code [zeigen / verbergen]

ggplot(data = fac1_tbl, 
       aes(x = animal, y = jump_length, fill = animal)) +
  theme_minimal() +
  geom_boxplot() +
  stat_summary(fun.y = mean, geom = "point", 
               shape=23, size = 5, fill = "gray50") +
  labs(x = "Flohart", y = "Sprungweite in [cm]") +
  theme(legend.position = "none") + 
  scale_fill_okabeito()

Abbildung 81.15— Beispielhafter einfaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten.

81.4.4 Zweifaktorieller Boxplot

Den zweifaktoriellen Boxplot erstellen wir für die einzelnen Floharten getrennt für die beiden Messorte. Du musst schauen, was du auf die x-Achse legst und was du dann auf die Legende und daher mit fill gruppierst. Gerne ergänze ich noch den Mittelwert mit der Funktion stat_summary(), muss hier aber schauen, dass ich nach dem Faktor animal gruppiere und dann noch mit der Funktion position_dodge() die richtige Position finde.

R Code [zeigen / verbergen]

ggplot(data = fac2_tbl, 
       aes(x = sex, y = jump_length, fill = animal)) +
  theme_minimal() +
  geom_boxplot() + 
  stat_summary(fun.y = mean, geom = "point", aes(group = animal), 
               shape=23, size = 5, fill = "gray50",
               position = position_dodge(0.75)) +
  labs(x = "Flohart", y = "Sprungweite in [cm]", fill = "Tierart") +
  scale_fill_okabeito()

Abbildung 81.16— Beispielhafter zweifaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Messorten.

81.5 Statistische Überprüfung

81.5.1 Normalverteilung

How to Calculate Skewness & Kurtosis in R

Die Schiefe (eng. skewness) ist ein Maß für die Asymmetrie einer Verteilung. Dieser Wert kann positiv oder negativ sein.

Eine negative Schiefe deutet darauf hin, dass sich der Schwanz auf der linken Seite der Verteilung befindet, die sich in Richtung negativer Werte erstreckt. Eine positive Schiefe zeigt an, dass sich der Schwanz auf der rechten Seite der Verteilung befindet, die sich in Richtung positiver Werte erstreckt. Ein Wert von Null bedeutet, dass die Verteilung überhaupt nicht schief ist, d. h. die Verteilung ist vollkommen symmetrisch.

Die Kurtosis (eng. kurtosis) ist ein Maß dafür, ob eine Verteilung im Vergleich zu einer Normalverteilung ein starkes oder schwaches Schwanzende aufweist.

Die Kurtosis einer Normalverteilung beträgt 3. Wenn eine gegebene Verteilung eine Kurtosis von weniger als 3 aufweist, wird sie als playkurtisch bezeichnet, was bedeutet, dass sie dazu neigt, weniger und weniger extreme Ausreißer zu produzieren als die Normalverteilung. Wenn eine gegebene Verteilung eine Kurtosis größer als 3 hat, wird sie als leptokurtisch bezeichnet, was bedeutet, dass sie dazu neigt, mehr Ausreißer als die Normalverteilung zu produzieren.

In einigen Formeln (Definition von Fisher) wird von der Kurtosis 3 abgezogen, um den Vergleich mit der Normalverteilung zu erleichtern. Nach dieser Definition hätte eine Verteilung eine größere Kurtosis als eine Normalverteilung, wenn der Kurtosis-Wert größer als 0 wäre.

Abbildung 81.17— Die t-Verteilung für drei beispielhafte Freiheitsgrade. Je größer die Freiheitsgrade und damit die Fallzahl, desto näher kommt die t-Verteilung einer Normalverteilung nahe. *[Zum Vergrößern anklicken]*

R Code [zeigen / verbergen]

pacman::p_load(moments, broom, scales, report)

#calculate skewness
skewness(fac2_tbl$hatch_time)

[1] 21.30159

R Code [zeigen / verbergen]

#calculate kurtosis
kurtosis(fac2_tbl$hatch_time)

[1] 496.1632

R Code [zeigen / verbergen]

fac2_tbl |> 
  group_by(animal, sex) |> 
  summarise(kurtosis = kurtosis(hatch_time),
            skewness = skewness(hatch_time))

`summarise()` has grouped output by 'animal'. You can override using the
`.groups` argument.

# A tibble: 6 × 4
# Groups:   animal [3]
  animal sex    kurtosis skewness
  <fct>  <fct>     <dbl>    <dbl>
1 cat    male      31.3      4.59
2 cat    female     9.23     2.46
3 dog    male      94.6      9.59
4 dog    female     7.22     1.87
5 fox    male       8.74     2.09
6 fox    female     9.13     2.27

Die Momente-Bibliothek bietet auch die Funktion jarque.test(), die einen Anpassungsgütetest durchführt, der feststellt, ob die Stichprobendaten eine Schiefe und eine Wölbung aufweisen, die einer Normalverteilung entsprechen oder nicht. Die Null- und Alternativhypothesen dieses Tests lauten wie folgt:

Nullhypothese: Der Datensatz weist eine Schiefe und Wölbung auf, die einer Normalverteilung entspricht.

Alternativhypothese: Der Datensatz weist eine Schiefe und eine Kurtosis auf, die nicht mit einer Normalverteilung übereinstimmen.

R Code [zeigen / verbergen]

jarque.test(fac2_tbl$hatch_time)


    Jarque-Bera Normality Test

data:  fac2_tbl$hatch_time
JB = 6125624, p-value < 2.2e-16
alternative hypothesis: greater

R Code [zeigen / verbergen]

fac2_tbl |> 
  split(~ animal + sex) |> 
  map(~jarque.test(.x$hatch_time)) |> 
  map(tidy) |> 
  bind_rows(.id = "test") |>
  select(test, p.value) |> 
  mutate(decision = ifelse(p.value <= 0.05, "reject normal", "normal"),
         p.value = pvalue(p.value, accuracy = 0.001))

# A tibble: 6 × 3
  test       p.value decision     
  <chr>      <chr>   <chr>        
1 cat.male   <0.001  reject normal
2 dog.male   <0.001  reject normal
3 fox.male   <0.001  reject normal
4 cat.female <0.001  reject normal
5 dog.female <0.001  reject normal
6 fox.female <0.001  reject normal

R Code [zeigen / verbergen]

fac2_tbl |>
  select(jump_length, hatch_time) |> 
  report()

The data contains 600 observations of the following 2 variables:

  - jump_length: n = 600, Mean = 20.51, SD = 3.77, Median = 20.51, MAD = 3.93,
range: [11.42, 29.75], Skewness = 0.02, Kurtosis = -0.52, 0% missing
  - hatch_time: n = 600, Mean = 251.61, SD = 793.13, Median = 145.00, MAD =
123.83, range: [3.17, 18761.35], Skewness = 21.36, Kurtosis = 497.31, 0%
missing

text

shapiro.test()

Wenn du überprüfen willst, ob dein Messwert $y$ einer Normalverteilung folgt, dann kannst du die Funktion check_normality() nutzt. Die Funktion rechnet dann den Shapiro-Wilk-Test um auf eine Abweichung von der Normalverteilung zu testen. Hierzu ist anzumerken, dass der Test relativ empfindlich bei Abweichungen in den Verteilungsschwänzen ist. Darüber hinaus braucht der Shapiro-Wilk-Test auch etwas Fallzahl, damit er auf die Normalverteilung testen kann. Im Folgenden schauen wir uns den Code für ein einfaktorielles, zweifaktorielles und dreifaktorielles Modell einmal an. Am Ende des Abschnitts gehe ich nochmal darauf ein, was du machen kannst, wenn du keine Normalverteilung in deinem Messwert $y$ vorliegen hast.

Einfaktoriell

Beginnen wir wieder mit einem einfaktoriellen Modell. Wir wollen wissen, ob unsere Sprungweite in [cm] über unsere verschiedenen Floharten normalverteilt ist. Wir bauen also erstmal das Modell und schicken es dann in die Funktion check_normality() aus dem R Paket {performance}.

R Code [zeigen / verbergen]

lm(jump_length ~ animal, data = fac1_tbl) |> 
  check_normality()

OK: residuals appear as normally distributed (p = 0.514).

Die Funktion sagt, dass wir eine Normalverteilung in unseren Daten vorliegen haben. Wir können uns auch einen Diagnoseplot wiedergeben lassen. Dafür müssen wir die Funktion nur an die Funktion plot() weiterleiten. Das Schöne ist, dass die Abbildung uns auch gleich sagt, was wir zu erwarten haben um eine Normalverteilung anzunehmen.

R Code [zeigen / verbergen]

lm(jump_length ~ animal, data = fac1_tbl) |> 
  check_normality() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.18— Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem einfaktoriellen Modell.

Zweifaktoriell

Im zweifaktoriellen Fall ändert sich jetzt nur das Model. Wir haben eben zwei Faktoren vorliegen und diese müssen wir dann mit ins Modell nehmen. Ich habe hier auch gleich den Interaktionsterm mit ergänzt, ich teste gerne das Modell, was ich dann später auch auswerten möchte.

Sprungweite jump_length
Schlupfzeit hatch_time

R Code [zeigen / verbergen]

lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_normality()

OK: residuals appear as normally distributed (p = 0.135).

R Code [zeigen / verbergen]

lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_normality()

Warning: Non-normality of residuals detected (p < .001).

Auch hier haben wir eine Normalverteilung in den Daten vorliegen. Gerne schaue ich mir auch die Abbildung der Residuen einmal an und das geht dann flott über die Funktion plot(). Da musst du nur die Ausgabe der Funktion check_normality() weiterleiten. Die leichten Bögen in den Punkten kommen von den unterschiedlichen Faktoren und deren Effekten auf die Sprungweite.

Sprungweite jump_length
Schlupfzeit hatch_time

R Code [zeigen / verbergen]

lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_normality() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.19— Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem zweifaktoriellen Modell.

R Code [zeigen / verbergen]

lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_normality() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.20— Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem zweifaktoriellen Modell.

81.5.2 Varianzhomogenität

bartlett.test()

leveneTest()

check_homogeneity()

Wir können auch die Funktion check_homogeneity() aus dem Paket {performance} nutzen. Wir erhalten hier auch gleich eine Entscheidung in englischer Sprache ausgegeben. Die Funktion check_homogeneity() nutzt den Bartlett-Test. Wir können in Funktion auch andere Methoden mit method = c("bartlett", "fligner", "levene", "auto") wählen.

Wir sehen, dass sich die Implementierung des Bartlett-Tests in check_homogeneity() nicht von der Funktion bartlett.test() unterscheidet, aber die Entscheidung gegen die Varianzhomogenität zu einem Signifikanzniveau von 5% gefällt wird. Nicht immer hilft einem der Entscheidungtext einer Funktion.

Nachdem wir die Normalverteilung des Messwertes $y$ getestet haben, wollen wir jetzt noch schauen, ob unsere Faktoren einer Varianzhomogenität über die Gruppen genügt. Die Funktion check_homogeneity() nutzt den Bartlett-Test um auf eine Abweichung von der Varianzhomogenität zu testen. Es gibt noch verschiedene andere Tests, aber ich würde hier empfehlen bei eiem Test zu bleiben und dann mit dem Ergebnis zu leben. Wild verschiedene Tests auf Varianzhomogenität durchzuführen führt dann auch nur zu einem wilden Fruchtsalat an Ergebnissen mit denen keiner was anfangen kann. Im Folgenden schauen wir uns den Code für ein einfaktorielles, zweifaktorielles und dreifaktorielles Modell einmal an. Am Ende des Abschnitts gehe ich nochmal darauf ein, was du machen kannst, wenn du keine Varianzhomogenität in deinen Faktoren vorliegen hast.

Einfaktoriell

Beginnen wir wieder mit dem einfaktoriellen Modell. Wir stecken das Modell dann einfach in die Funktion check_homogeneity() und erhalten die Information über die Varianzhomogenität wiedergegeben.

R Code [zeigen / verbergen]

lm(jump_length ~ animal, data = fac1_tbl) |> 
  check_homogeneity()

OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.297).

Wunderbar, wir haben keine Abweichung von der Varianzhomongenität. Wir können uns auch die Daten nochmal anschauen. Hier sehen wir aber schon, dass die Daten etwas heterogen aussehen der Test aber die Homogenität nicht ablehnt. Das ist immer schwierig bei kleinen Fallzahlen. Aber an irgendwas muss man eben glauben.

R Code [zeigen / verbergen]

lm(jump_length ~ animal, data = fac1_tbl) |> 
  check_homogeneity() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.21— Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem einfaktoriellen Modell.

Zweifaktoriell

Jetzt können wir mal schauen, was passiert wenn wir die Anzahl an möglichen Faktorkombinationen erhöhen indem wir ein zweifaktorielles Modell nutzen. Hier haben wir dann ja sechs Faktorkombinationen oder Gruppen die dann alle homogen in den Varianzen sein müssen.

Sprungweite jump_length
Schlupfzeit hatch_time

R Code [zeigen / verbergen]

lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_homogeneity()

OK: There is not clear evidence for different variances across groups (Bartlett Test, p = 0.651).

R Code [zeigen / verbergen]

lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_normality()

Warning: Non-normality of residuals detected (p < .001).

Wir haben auch hier Varianzhomogenität über alle Gruppen der Faktoren vorliegen. Wenn du dir jetzt die Abbildung zu dem Test anschaust, dann siehst du auch hier, dass die Violinplots eben dann doch alle etas anders aussehen. Wir haben aber hier auch das gleiche Problem wie bei dem einfaktoriellen Fall, wir haben eben dann doch recht wenig Fallzahl in unseren Daten.

Sprungweite jump_length
Schlupfzeit hatch_time

R Code [zeigen / verbergen]

lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_homogeneity() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.22— Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem zweifaktoriellen Modell.

R Code [zeigen / verbergen]

lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> 
  check_homogeneity() |> 
  plot() +
  scale_fill_okabeito()

Abbildung 81.23— Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem zweifaktoriellen Modell.

81.6 Auswege

Flowchart

text

flowchart TD
    A("Normalverteiltes Outcome
       in jeder Versuchsgruppe"):::factor --- B(((ja))) --> B1 
    A("Normalverteiltes Outcome
       in jeder Versuchsgruppe"):::factor --- F(((nein))) --> B2
    subgraph B1["Mittelwertsvergleiche"]
    C("Varianzhomogenität
       über alle Gruppen"):::factor --- D((("ja"))) --> E("lm()"):::factor
    C("Varianzhomogenität
       über alle Gruppen"):::factor --- J(((nein))) --> K("gls()"):::factor 
    end
    subgraph B2["Medianvergleiche oder Anteile"]
    G("Varianzhomogenität
       über alle Gruppen"):::factor --- H(((ja))) --> I("rq()"):::factor  
    G("Varianzhomogenität
       über alle Gruppen"):::factor --- L(((nein))) 
    end
    classDef factor fill:#56B4E9,stroke:#333,stroke-width:0.75px

Abbildung 81.24— Flowchart für die Entscheidung welches statistische Modell gerechnet werden kann: lm(), lineare Regression, gls(), generalized least squares Regression, rq(), Quantilesregression. Die Funktionen finden sich teilweise in eigenen R Paketen. Bei einem nicht-normalverteilten Outcome mit Varianzheterogenität über die Gruppen müssen wir nochmal gemeinsam in die Daten schauen.

Du findest die vorgeschlagenen Funktionen dann in den entsprechenden Kapiteln zur ANOVA und den Posthoc-Tests. Du kannst dir dann dort den Code zusammensuchen. Je nach deiner Datenlage musst du dann nochmal etwas an dem R Code programmieren. Beachte, dass die Funktionen sich teilweise in eigenen R Paketen finden lassen. So ist die Funktion gls() im R Paket {nlme} und die Funktion rq() im R Paket {quantreg} zu finden. Du kannst auch bei Varianzheterogenität das R Paket {sandwich} nutzen und einen entsprechend angepassten Varianzschätzer. Mehr findest du dazu bei den Posthoc-Test in dem Abschnitt zu dem Gruppenvergleich unter Varianzheterogenität oder gleich im ersten Zerforschenbeispiel zum einfaktoriellen Barplot.

81.6.1 Transformation

R Code [zeigen / verbergen]

log_tbl <- fac2_tbl |> 
  mutate(log_hatch_time = log(hatch_time))

Abbildung 81.25— Histogramm der nicht transfomierten und transformierten Daten. **(A)** Nicht transformierte, rohe Daten. **(B)** $log$-transformierte Daten. *[Zum Vergrößern anklicken]*

R Code [zeigen / verbergen]

rank_tbl <- fac2_tbl |> 
  mutate(rank_hatch_time = rank(hatch_time))

Abbildung 81.26— Histogramm der nicht transfomierten und transformierten Daten. **(A)** Nicht transformierte, rohe Daten. **(B)** $log$-transformierte Daten. *[Zum Vergrößern anklicken]*

R Code [zeigen / verbergen]

jump_mod <- lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl)

R Code [zeigen / verbergen]

bc_obj <- boxcox(jump_mod, plotit = FALSE)
lambda <- bc_obj$x[which.max(bc_obj$y)]
lambda

[1] 0

R Code [zeigen / verbergen]

boxcox_tbl <- fac2_tbl |> 
  mutate(jump_boxcox = ((jump_length ^ lambda - 1)/lambda))

R Code [zeigen / verbergen]

ggplot(data = boxcox_tbl) +
  theme_minimal() +
  geom_histogram(aes(x = jump_boxcox), alpha = 0.9,
                 fill = cbbPalette[3], color = "black") +
  labs(x = 'Box-Cox transformierte Sprungweiten', y = 'Anzahl')

Warning: Removed 600 rows containing non-finite outside the scale range
(`stat_bin()`).

Abbildung 81.27— Histogramm der Box-Cox transfomierten Daten in `ggplot` dargestellt.

R Paket {MASS}

81.6.2 Modellierung

R Paket `{performance}`

Wenn du noch etwas weiter gehen möchtest, dann kannst du dir noch die Hilfeseite von dem R Paket {performance} Robust Estimation of Standard Errors, Confidence Intervals, and p-values anschauen. Die Idee ist hier, dass wir die Varianz/Kovarianz robuster daher mit der Berücksichtigung von Varianzheterogenität (eng. heteroskedasticity) schätzen.

R Code [zeigen / verbergen]

fit_fac2 <- lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl)

R Code [zeigen / verbergen]

fit_fac2 |> 
  model_parameters()

Parameter                   | Coefficient |   SE |         95% CI | t(594) |      p
-----------------------------------------------------------------------------------
(Intercept)                 |       15.38 | 0.19 | [15.00, 15.76] |  79.46 | < .001
animal [dog]                |        2.73 | 0.27 | [ 2.19,  3.27] |   9.98 | < .001
animal [fox]                |        5.33 | 0.27 | [ 4.79,  5.87] |  19.47 | < .001
sex [female]                |        5.08 | 0.27 | [ 4.54,  5.61] |  18.55 | < .001
animal [dog] × sex [female] |       -0.24 | 0.39 | [-1.00,  0.52] |  -0.62 | 0.535 
animal [fox] × sex [female] |       -0.28 | 0.39 | [-1.04,  0.48] |  -0.71 | 0.475


Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

R Code [zeigen / verbergen]

fit_fac2 |> 
  model_parameters(vcov = "HC3")

Parameter                   | Coefficient |   SE |         95% CI | t(594) |      p
-----------------------------------------------------------------------------------
(Intercept)                 |       15.38 | 0.19 | [15.01, 15.74] |  82.64 | < .001
animal [dog]                |        2.73 | 0.27 | [ 2.21,  3.26] |  10.23 | < .001
animal [fox]                |        5.33 | 0.27 | [ 4.80,  5.86] |  19.72 | < .001
sex [female]                |        5.08 | 0.26 | [ 4.57,  5.58] |  19.69 | < .001
animal [dog] × sex [female] |       -0.24 | 0.38 | [-0.99,  0.51] |  -0.63 | 0.530 
animal [fox] × sex [female] |       -0.28 | 0.38 | [-1.03,  0.48] |  -0.72 | 0.472


Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

R Paket `{emmeans}`

R Code [zeigen / verbergen]

fit_fac2 |>
  emmeans(~ animal * sex, vcov. = sandwich::vcovHAC)

 animal sex    emmean    SE  df lower.CL upper.CL
 cat    male     15.4 0.188 594     15.0     15.7
 dog    male     18.1 0.176 594     17.8     18.5
 fox    male     20.7 0.197 594     20.3     21.1
 cat    female   20.5 0.181 594     20.1     20.8
 dog    female   22.9 0.189 594     22.6     23.3
 fox    female   25.5 0.211 594     25.1     25.9

Confidence level used: 0.95

R Paket `{nlme}`

R Code [zeigen / verbergen]

fit_fac2 <- gls(jump_length ~ animal + sex + animal:sex, data = fac2_tbl,
                weights = varIdent(form = ~ 1 | animal*sex))

R Code [zeigen / verbergen]

fit_fac2 |> 
  model_parameters()

# Fixed Effects

Parameter                   | Coefficient |   SE |         95% CI | t(594) |      p
-----------------------------------------------------------------------------------
(Intercept)                 |       15.38 | 0.19 | [15.01, 15.74] |  83.05 | < .001
animal [dog]                |        2.73 | 0.27 | [ 2.21,  3.25] |  10.29 | < .001
animal [fox]                |        5.33 | 0.27 | [ 4.80,  5.86] |  19.82 | < .001
sex [female]                |        5.08 | 0.26 | [ 4.57,  5.58] |  19.78 | < .001
animal [dog] × sex [female] |       -0.24 | 0.38 | [-0.99,  0.51] |  -0.63 | 0.528 
animal [fox] × sex [female] |       -0.28 | 0.38 | [-1.03,  0.47] |  -0.72 | 0.470


Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

Referenzen

```{r echo = FALSE} #| message: false #| warning: false pacman::p_load(tidyverse, readxl, knitr, kableExtra, Hmisc, quantreg, multcomp, emmeans, ggpubr, multcompView, nlme, tinytable, see, patchwork, conflicted) conflicts_prefer(dplyr::filter) conflicts_prefer(dplyr::select) ## gg_template <- ggplot() + theme_minimal() + theme(#axis.text.x = element_blank(), #axis.ticks.x = element_blank(), # axis.text.y = element_text(color = "#CC79A7", face = "bold", size = 14), axis.text.y = element_text(), axis.ticks.y = element_blank(), axis.text = element_text(size = 12), #axis.title = element_text(size = 14, face = "bold"), plot.title = element_text(size = 16, face = "bold"), plot.subtitle = element_text(size = 12, face = "italic"), plot.caption = element_text(face = "italic"), legend.position = "none") + scale_color_okabeito() + scale_fill_okabeito() ## cbbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7") ## gummi_tbl <- read_excel("data/gummibears.xlsx") |> select(gender, height) |> na.omit() |> mutate(gender = as_factor(gender)) ## stat_gummi_tbl <- gummi_tbl |> group_by(gender) |> summarise(n = n()) ``` # Der Pre-Test oder Vortest {#sec-pretest} *Letzte Änderung am `r format(fs::file_info("stat-tests-pretest2.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`* > *"I struggle with some demons; They were middle class and tame." --- Leonard Cohen, You Want It Darker* In diesem Kapitel soll es um den Pre-Test oder auch Vortest gehen. Wir sind hier in einem experimentellen Design, welches verschiedene Gruppen beinhaltet oder aber wir wollen wissen, ob wir in einer linearen Regression die Normalverteilung unseres Messwertes $y$ vorliegen haben. Grundsätzlich geht es erstmal darum herauszufinden, ob die Annahmen an einen statistischen Test in deinen Daten erfüllt sind. Häufig wollen wir eine [ANOVA für einen Gruppenvergleich](#sec-anova) rechnen und dann anschließend einen [multiplen Test oder Post-hoc Test](#sec-posthoc) durchführen. In beiden Fällen wird es einfacher wenn wir eine Normalverteilung in unseren Messwert $y$ sowie eine Varianzhomogenität in unseren Behandlungsgruppen oder Faktoren $f$ vorliegen haben. Mit einfacher meine ich, dass du auch mit einer Abweichung von der Normalvertielung und auch Varianzhterogenität heutzutage umgehen kannst. Der Standard im statistischen Testen war aber immer die Normalverteilung und die Varianzhomogenität. Wenn beides nicht vorlag, dann wurde es manchmal etwas dunkel. Wir aber aber im 21. Jahrhundert ein paat mehr Pfeile im Köcher und können mit unterschiedlichsten Daten und deren Eigenschaften umgehen. Mehr dazu findest du im [Teil zum statistsichen Modellieren](https://jkruppa.github.io/stat-modeling-preface.html) und den nachfolgenden Kapiteln. ## Allgemeiner Hintergrund Wir werden uns in diesem Kapitel auf das faktorielle Experiment konzentrieren. Natürlich kannst du auch alle Funktionen in einem anderen Design anwenden. Wenn du wissen willst, ob eine Variable normalverteilt ist oder aber ein Gruppenfaktor homogen in den Varianzen, dann helfen dir hier auch die Funktionen weiter. Häufig werden aber die beiden Eigenschaften Normalverteilung und Varianzhomogenität in Gruppenvergleichen verwendet. Eine Sache ist aber wichtig zu wissen. Wir untersuchen in unseren Experimenten ja immer nur eine Stichprobe der Grundgesamtheit und wollen dann von der Stichprobe einen Rückschluß auf die Grundgesamt machen. Wenn dich mehr dazu interessiert, dann schaue einmal in dem [Kapitel zum Testen von Hypothesen](https://jkruppa.github.io/stat-tests-preface-theory.html) rein. Es kann also sein, dass wir definitiv in der Grundgesamtheit einen normalverteilten Messwert vorliegen haben, wir aber noch zu wenige Beobachtungen in unsere Stichprobe vorliegen haben um diese Normalverteilung in einem Histogramm oder Densityplot zu sehen. Nehmen wir einmal die Körpergröße als ein normalverteilten Messwert $y$ an. Wir wissen, dass die Körpergröße einer Normalverteilung folgt. In der folgenden Abbildung @fig-pretest-gummi-01 siehst du einmal die Körpergrößen von unseren Gummibärchendaten. Insgesamt haben $`r pull(filter(stat_gummi_tbl, gender == "m"), n)`$ Männer und $`r pull(filter(stat_gummi_tbl, gender == "w"), n)`$ Frauen bei der Gummibärchenumfrage mitgemacht. Dennoch beobachten wir keine saubere Normalverteilung, wie wir sie erwarten würden. Wir haben noch zu wenige Beobachtungen gemacht. ```{r} #| echo: false #| message: false #| warning: false #| fig-align: center #| fig-height: 5 #| fig-width: 10 #| label: fig-pretest-gummi-01 #| fig-cap: "Darstellung der Körpergröße in [cm] für die Geschlechter getrennt. Die Körpergröße ist normalverteilt. Die Farben repräsentieren die jeweiligen Geschlechter. Die Männer sind blau und die Frauen in lila dargestellt. **(A)** Histogramm. **(B)** Densityplot. *[Zum Vergrößern anklicken]*" mean_tbl <- gummi_tbl |> group_by(gender) |> summarise(mean = round(mean(height), 1)) y_max <- 55 p1 <- ggplot(data = gummi_tbl, aes(x = height, fill = gender)) + theme_minimal() + geom_histogram(position = "dodge", color = "black") + labs(x = "", y = "Anzahl", fill = "Geschlecht") + scale_x_continuous(breaks = seq(150, 210, by = 5)) + scale_y_continuous(breaks = seq(0, y_max, 5), limits = c(0, y_max)) + scale_fill_okabeito(order = c(2, 7), labels = c("männlich", "weiblich")) + theme(legend.position = "none") #c(0.85, 0.8)) p2 <- ggplot(data = gummi_tbl, aes(x = height, fill = gender)) + theme_minimal() + geom_density(alpha = 0.75) + labs(x = "", y = "", fill = "Geschlecht") + scale_fill_okabeito(order = c(2, 7), labels = c("männlich", "weiblich")) + scale_x_continuous(breaks = seq(145, 210, by = 5), limits = c(145, 210)) + scale_y_continuous(breaks = seq(0, 0.055, 0.01), limits = c(0, 0.055)) + geom_vline(xintercept = mean_tbl$mean, color = cbbPalette[c(6, 8)], size = 1) + annotate("label", x = mean_tbl$mean, y = 0.055, label = mean_tbl$mean) + theme(legend.position = "none") #c(0.85, 0.8)) p1 + p2 + plot_layout(ncol = 2) + plot_annotation(tag_levels = 'A') + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` Wir sehen also, nur weil etwas wie die Körpergröße wirklich normalverteilt *ist*, ist es noch etwas ganz anders diese Normalverteilung dann auch in den Messwert $y$ zu *beobachten*. Die Fallzahlen in der Grundgesamt und in der Stichprobe unterscheiden sich dann doch gewaltig und wir sind dann eben auch auf Annahmen angewiesen. Meistens passt es auch mit den Annahmen und wenn wir mal daneben liegen, kann es sein, dass es dann doch nicht so viel ausmacht, wenn der Effekt in unserer statistischen Auswertung groß genug ist. #### Das Modell {.unnumbered .unlisted} Auch hier möchte ich einmal das statistische Modell besprechen was wir in dem Gruppenvergleich oder dem statistischen Modellieren benötigen. Im Folgenden findest du einmal ein faktorielles Modell mit einem Messwert $y$ und zwei Gruppenfaktoren $f_A$ und $f_B$. Die beiden Faktoren entsprechen zwei unterschiedlichen kategoriellen Variablen mit verschiedenen Gruppen. Wir wollen uns ja in diesem Kapitel auf die Normalverteilung und die Varianzhomogenität konzentrieren. Die beiden Gütekriterien können aber ganz klar dem Messwert $y$ und den experimentellen Faktoren zugeordnet werden. $$ \underbrace{\;\mbox{Messwert}\; y\;}_{normalverteilt} \sim \overbrace{\;\mbox{Faktor}\; f_A + \mbox{Faktor}\; f_B\;}^{homogene\; Varianzen} $$ mit - $\mbox{Messwert}\; y$, gleich dem Messwert oder Outcome, wie die Sprungweite in \[cm\] als `jump_length` in unseren Beispieldaten. - $\mbox{Faktor}\; f_A$, gleich dem ersten Faktor $f_A$, wie die Tierart als `animal` mit unterschiedlichen Gruppen oder Leveln. - $\mbox{Faktor}\; f_B$, gleich dem zweiten Faktor $f_B$, wie der Messort als `site` mit unterschiedlichen Gruppen oder Leveln. Damit haben wir uns erstmal für die Vortest für die Normalverteilung und die Varianzhomogenität geordnet. Wir wollen dann in den folgenden Abschnitten noch andere Gütekriterien eines Modells kurz anreißen, aber den Hauptteil findest du im [Kapitel zur Modelgüte](#sec-lin-reg-quality) von linearen Modellen. #### Gibt es noch mehr Vortests? {.unnumbered .unlisted} Jetzt könnte man meinen, dass mit der Normalverteilung und der Varianzhomogenität eigentlich die wichtigsten Gütekriterien vorgetestet werden. Es gibt aber für lineare Modelle, was ein Gruppenvergleich dann am Ende auch nur ist, noch andere Gütekriterien. Neben diesen beiden Eigenschaften können wir usn auch noch folgende weitere Gütekriterien anschauen. Ich verweise hier einmal auf die [Hilfeseite des R Packetes `{performance}`](https://easystats.github.io/performance/articles/check_model.html) für mehr Informationen und deren [Referenzseite der Familie der `check_*()` Funktionen](https://easystats.github.io/performance/reference/index.html#check-model-assumptions-or-data-properties). Wie immer kommt es dann auf die Fragestellung und dann auf das enstprechende Modell sowie den verwendeten Algorithmus an. Je nachdem was du gemessen hast, also welche Werte dein $y$ annimmt, musst du einen anderes Modell wählen. Je nach Modell hast du dann auch andere Annahmen. Das würde hier aber das Kapitel sprengen. Gerne kannst du als Startpunkt einmal in das [Teil zum statistsichen Modellieren](https://jkruppa.github.io/stat-modeling-preface.html) reinschauen. Betrachten wir also einmal im Folgenden die beiden wichtigsten Annahmen an ein faktorielles Design oder aber Gruppenvergleich. Wir fragen uns, haben wir eine Normalverteilung in den Messwerten $y$ und homogene Varianzen in den Faktoren oder Gruppen $f$ vorliegen? Dann können wir ganz normal eine ANOVA oder einen Tukey HSD Test rechnen. ### Normalverteilung Fangen wir also mit der Annahme der Normalverteilung an die Daten an. Hierbei ist wichtig, dass wir nicht die Daten insgesamt betrachten sondern uns fragen, ob der betrachtete Messwert $y$ im Modell oder statistischen Test normalverteilt ist. Wir haben uns den Zusammenhang ja schon oben einmal in dem statistischen Modell angeschaut. Häufig führt dies zu Verwirrungen, da verallgemeinert von den Daten gesprochen wird, die normalverteilt sein soll. Hier geht es dann wirklich *nur* um deinen Messwert $y$. Das nochmal als Erinnerung für die weiteren Betrachtungen. Was wären also beispielhaft normalverteilte Messwerte? *Beispiel: Frischgewicht, Trockengewicht, Chlorophyllgehalt, Pflanzenhöhe* | freshmatter | drymatter | chlorophyll | height | |:-----------:|:---------:|:-----------:|--------| | 8.23 | 1.21 | 45.88 | 24.19 | | 2.61 | 0.87 | 43.91 | 18.51 | | 4.81 | 0.34 | 37.44 | 21.74 | : Tabelle mit beispielhaften, normalverteilten Messwerten $y$. {#tbl-bsp-normal} Nach dem [zentralen Grenzwertsatz](https://de.wikipedia.org/wiki/Zentraler_Grenzwertsatz) können wir bei Merkmalen, die sich aus verschiedenen Einflussfaktoren zusammensetzen, allgemein von einer Normalverteilung ausgehen. Die Körpergröße oder das Körpergewicht ist normalverteilt, da wir hier es mit vielen Einflussgrößen zu tun haben, die das tatsächliche Körpergewicht einer Beobachtung ausmachen. Das Körpergewicht hängt eben von der täglichen Kalorienmenge, verschiedensten Genen, dem Muskelanteil, dem Aktivitätsgrad, der sozialen Stellung und vielen weiteren Einflusfaktoren ab. Alles zusammen addiert sich dann zum Körpergewicht wobei jeder Einflussfaktor nur einen kleinen Teil ausmacht. Was heißt *approximativ normalverteilt*? : Wir sprechen von approximativ normalverteilt, wenn wir meinen, dass ein Messwert $y$ in unserer Stichprobe annähernd normalverteilt ist. Wir sind uns also nicht zu hundertprozent sicher, glauben aber das die Normalverteilungsannahme an unseren Messwert passen wird. Häufig sagen wir auch, dass gewisse Tranformationen approximativ normalverteilt sind. So haben wir nach einer log-Transformation log-normalverteilte Messwerte vorliegen. Wir sagen dann meistens, dass ein log-transformierter Messwert approximativ normalverteilt ist. ### Varianzhomogenität Kommen wir nun zur Varianzhomogenität oder Varianzheterogenität in den Gruppen des Behandlunsgfaktors. Je nachdem was du betrachtest, nennen wir es eben Varianzhomogenität oder Varianzheterogenität. Entweder sind die Varianzen gleich, dann haben wir Varianzhomogenität vorliegen oder die Varianzen in den Gruppen sind nicht gleich, dann hast du Varianzheterogenität in den Daten. Es gibt so ein paar Daumenregeln, die dir helfen abzuschätzen, ob in deinen Gruppen Varianzheterogenität vorliegt. Um es kurz zu machen, vermutlich hast du mindestens leichte Varianzheterogenität in den Daten vorliegen. Es ist bei kleinen Gruppengrößen nicht zu vermeiden, dass sich die Varianzen eben dann doch unterscheiden. 1) Du hat viele Behandlungsgruppen. Je mehr Gruppen du hast oder eben dann auch Faktorkombinationen, die du testen möchtest, desto wahrscheinlicher wird es, dass mindestens eine Gruppe eine unterschiedliche Varianz hat. Du hast Varianzheterogenität vorliegen. 2) Du misst deine Gruppen über die Zeit. Je größer, schwerer oder allgemein höher ein Messwert wird, desto größer wird auch die Varianz. Schaust du dir deine Messwerte über die Zeit an hast du meistens Varianzheterogenität vorliegen. 3) Deine Kontrolle ohne Behandlung verhält sich meistens nicht so, wie die Gruppen, die eine Behandlung erhalten haben. Wenn du nichts machst in deiner negativen Kontrolle, dann hast du meistens eine andere Streuung der Messwerte als unter einer Behandlung. 4) Deine Behandlungen sind stark unterschiedlich. Wenn deine Behandlungen sich biologisch oder chemisch in der Wirkung unterscheiden, denn werden vermutlich deine Messwerte auch anders streuen. Hier spielt auch die Anwendung der Behandlung und deren Bereitstellung eine Rolle. Wenn was nicht gleich ist, dann wird es vermutlich nicht gleiche Messwerte erzeugen. 5) Du hast wenig Fallzahl pro Gruppe oder Faktorkombination. Wenn du wenig Fallzahl in einer Gruppe hast, dann reicht schon eine (zufällige) Messabweichung und schon sind deine Varianzen heterogen. In den folgenden Beispielen habe ich dir einige Experimente mit einem faktoriellen Design mitgebracht. Die Fotos stammen aus wissenschaftlichen Publikationen wie einer wissenschaftlichen Veröffentlichung oder aber wissenschaftlichen Postern hier auf dem Gelände der Hochschule Osnabrück. ::: panel-tabset ## Beispiel 1 ![Ein zweifaktorielles Experiment mit neuen Faktorkombinationen die alle miteinander paarweise vergleichen werden. Wir sehen sehr gut, dass die Kontrolle sehr kleine Werte hat und somit die Varianz in der Kontrolle sehr viel kleiner ist alles in den anderen Gruppen.](images/eda/zerforschen_barplot_2fac.png){#fig-pretest-barplot-01 fig-align="center" width="100%"} ## Beispiel 2 ![Ein zweifaktorielles Experiment mit sechs Faktorkombinationen. Wir sehen hier sehr gut, dass mit steigenden MIttelwerten, also höheren Barplots, auch die Varianz in den Gruppen zunimmt. Die Fehlerbalken werden immer länger.](images/eda/zerforschen_barplot_selen.jpg){#fig-pretest-barplot-02 fig-align="center" width="100%"} ## Beispiel 3 ![Ein zweifaktorielles Experiment mit sehr vielen Faktorrkombinationen. Durch die unterschiedlichen miitelren Zählwerte ergeben sich sehr viele unterschiedlich große Mittelwerte. Darüber hinaus haben wir sehr viele Gruppen. Wir sehen hier sehr viel Varianzhterogenität.](images/eda/zerforschen_barplot_pest_count.jpg){#fig-pretest-barplot-02 fig-align="center" width="100%"} ## Beispiel 4 ![Zwei einfaktorielle Experiemente in einer Abbildung dargestellt. Die linken Barplots und die rechten Barplots wurden getrennt voneinander ausgewertet. Auch hier sieht man sehr viel unterschiedliche Streuung in den Daten.](images/eda/zerforschen_barplot_root.jpg){#fig-pretest-barplot-02 fig-align="center" width="100%"} ::: Jetzt haben wir uns einmal die wichtigsten Abbildungen angeschaut und haben so eine erste Idee was Varianzhomogenität sein könnte. Wir schauen uns dann in den folgenden Abschnitten noch mehr zu der Bestimmung an. Dann bleibt eigentlich noch eine abschließende Frage für den einführende Abschnitt. Tut Varianzheterogenität anstatt Varianzhomogenität weh? : Nein. Meistens ist die Varianzheterogenität nicht so ausgeprägt, dass du nicht auch eine ANOVA oder anderen statistischen Test rechnen kannst. Über alle Gruppen hinweg wird dann zwar zum Beispiel in einer ANOVA die Varianz gemittelt und es kann dann zu weniger signifikanten Ergebnissen führen, aber so schlimm ist es nicht. Im Post-hoc Test solltest du aber die Varianzheterogenität berücksichtigen, da du ja immer nur zwei Gruppen gleichzeitig betrachtest. Aber auch hier gibt es dann die passenden Adjustierungen. Mehr dazu am Ende des Kapitels im Abschnitt zu den Auswegen. ## Genutzte R Pakete Wir wollen folgende R Pakete in diesem Kapitel nutzen. ```{r echo = TRUE} #| message: false pacman::p_load(tidyverse, magrittr, olsrr, broom, car, performance, see, scales, readxl, nlme, parameters, conflicted) conflicts_prefer(dplyr::filter) conflicts_prefer(dplyr::select) ``` An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst. ## Daten Wir immer bringe ich hier ein paar Datensätze mit damit wir dann verstehen, was eigentlich in den folgenden Analysen in R und den entsprechenden R Paketen passiert. Ich zeige hier an den Daten nur die Anwendung in R. Deshalb fehlen dann hier auch die Mittelwerte und andere deskriptive Maßzahlen. Schauen wir jetzt also mal in unsere Beispieldaten für die einfaktorielle und zweifaktorielle Datenanalyse rein. #### Einfaktorieller Datensatz {.unnumbered .unlisted} Beginnen wir mit einem einfaktoriellen Datensatz. Wir haben hier als Messwert die Sprungweite von Flöhen in \[cm\] vorliegen. Wissen wollen wir, ob sich die Sprungweite für drei verschiedene Floharten unterscheidet. Damit ist dann in unserem Modell der Faktor `animal` und die Sprungweite `jump_length` als Messwert. Ich lade einmal die Daten in das Objekt `fac1_tbl`. Hier haben wir dann ein simples Design vorliegen. ```{r} #| message: false fac1_tbl <- read_xlsx("data/flea_dog_cat_fox.xlsx") |> select(animal, jump_length) |> mutate(animal = as_factor(animal)) ``` Dann schauen wir uns die Daten einmal in der folgenden Tabelle als Auszug einmal an. Wichtig ist hier nochmal, dass du eben einen Faktor `animal` mit drei Leveln also Gruppen vorliegen hast. Wir wollen jetzt die drei Tierarten hinsichtlich ihrer Sprungweite in \[cm\] miteinander vergleichen. ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-1fac-table #| tbl-cap: "Tabelle der Sprungweiten in [cm] als Messwert $y$ von Hunde-, Katzen- und Fuchsflöhen. Der Datensatz ist einfaktoriell, da wir nur einen Faktor vorliegen haben." fac1_raw_tbl <- read_xlsx("data/flea_dog_cat_fox.xlsx") |> select(animal, jump_length) rbind(head(fac1_raw_tbl, n = 3), rep("...", times = ncol(fac1_raw_tbl)), tail(fac1_raw_tbl, n = 3)) |> tt(width = 2/3, align = "c", theme = "striped") ``` Und dann wollen wir uns noch einmal die Daten als einen einfachen Boxplot anschauen. Wir sehen, dass die Daten so gebaut sind, dass wir einen signifikanten Unterschied zwischend den Sprungweiten der Floharten erwarten. Die Boxen der Boxplots überlappen sich nicht und die Boxplots liegen auch nicht auf einer Ebene. Wir könnten hier von normalverteilten Daten und Varianzhomogenität ausgehen. Die Mediane liegen in der Mitte der Boxen und die Boxen sind ungefähr gleich groß. ```{r} #| message: false #| echo: false #| warning: false #| label: fig-ggplot-anova-boxplot-1fac #| fig-align: center #| fig-height: 4 #| fig-width: 4 #| fig-cap: "Beispielhafter einfaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten." ggplot(data = fac1_tbl, aes(x = animal, y = jump_length, fill = animal)) + theme_minimal() + geom_boxplot() + stat_summary(fun.y = mean, geom = "point", shape=23, size = 3, fill = "gray50") + labs(x = "Flohart", y = "Sprungweite in [cm]") + theme(legend.position = "none") + scale_fill_okabeito() ``` #### Zweifaktorieller Datensatz {.unnumbered .unlisted} Neben dem einfaktoriellen Datensatz wollen wir uns noch den häufigeren Fall mit zwei Faktoren anschauen. Wir haben also nicht nur die drei Floharten vorliegen und wollen wissen ob diese unterschiedlich weit springen. Darüber hinaus haben wir noch einen zweiten Faktor gewählt. Wir haben die Sprungweiten der Hunde-, Katzen- und Fuchsflöhe nämlich für die beiden Geschlechter getrennt gemessen. Dadurch haben wir jetzt den Faktor `animal` und den Faktor `sex` vorliegen. Wiederum fragen wir uns, ob sich die Sprungweite in \[cm\] der drei Floharten in den beiden Geschlechtern unterscheidet. Darüber hinaus haben wir neben der Sprungweite noch die Schlupfzeiten in \[m\] gemessen. Im Folgenden lade ich einmal den Datensatz in das Objekt `fac2_tbl`. ```{r} #| message: false fac2_tbl <- read_xlsx("data/flea_dog_cat_length_weight.xlsx") |> select(animal, sex, jump_length, hatch_time) |> mutate(animal = as_factor(animal), sex = as_factor(sex)) ``` Betrachten wir als erstes einen Auszug aus der Datentabelle. Wir haben hier als Messwert oder Outcome $y$ die Sprungweite `jump_length` sowie die Schlupfzeiten `hatch_time` vorliegen. Als ersten Faktor die Variable `animal` und als zweiten Faktor die Variable `sex` festgelegt. ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-2fac-table #| tbl-cap: "Tabelle der Sprungweiten in [cm] und Schlupfzeiten [m] als Messwert $y$ von Hunde-, Katzen- und Fuchsflöhen der beiden Geschlechter. Der Datensatz ist zweifaktoriell, da wir einen Behandlungsfaktor mit `animal` und einen zweiten Faktor mit `sex` vorliegen haben." fac2_raw_tbl <- read_xlsx("data/flea_dog_cat_length_weight.xlsx") |> select(animal, sex, jump_length, hatch_time) rbind(head(fac2_raw_tbl, n = 3), rep("...", times = ncol(fac2_raw_tbl)), tail(fac2_raw_tbl, n = 3)) |> tt(width = 1, align = "c", theme = "striped") ``` Auch hier schauen wir uns einmal die Daten in einem Boxplot und einem Densityplot an. Wir wollen ja sehen, ob sich zum einen die Gruppen unterscheiden und zum anderen wie unsere Messwerte der Sprungweiten und der Schlupfzeiten verteilt sind. Wir erkennen in den Boxplots und auch in den Densityplots, dass wir vermutlich eine approximative Normalverteilung in den Sprungweiten vorliegen haben, aber auf keinen Fall eine Normalverteilung in den Schlupfzeiten. Du siehst hier nochmal in den beiden Abbildungen die Schiefe in der Verteilung der Schlupfzeiten. ::: panel-tabset ## Boxplot ```{r} #| message: false #| echo: false #| warning: false #| label: fig-ggplot-anova-boxplot-2fac-1 #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Zweifaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Geschlechtern." ggplot(data = fac2_tbl, aes(x = animal, y = jump_length, fill = sex)) + theme_minimal() + geom_boxplot() + labs(x = "Flohart", y = "Sprungweite in [cm]", fill = "Geschlecht") + scale_fill_okabeito() ``` ```{r} #| message: false #| echo: false #| warning: false #| label: fig-ggplot-anova-boxplot-2fac-2 #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Zweifaktorieller Boxplot für die Schlupfzeiten in [m] gruppiert nach den Floharten und den beiden Geschlechtern." ggplot(data = fac2_tbl, aes(x = animal, y = hatch_time, fill = sex)) + theme_minimal() + geom_boxplot() + labs(x = "Flohart", y = "Schlupfzeiten in [m]", fill = "Geschlecht") + scale_fill_okabeito() + ylim(0, 2000) ``` ## Densityplot ```{r} #| message: false #| echo: false #| warning: false #| label: fig-ggplot-pretest-dens-2fac-1 #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Densityplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Geschlechtern." ggplot(data = fac2_tbl, aes(x = jump_length, fill = interaction(animal, sex))) + theme_minimal() + geom_density(alpha = 0.25) + labs(x = "Sprungweite in [cm]", y = "", fill = "Tierart & Geschlecht") + scale_fill_okabeito() + xlim(9, 32) ``` ```{r} #| message: false #| echo: false #| warning: false #| label: fig-ggplot-pretest-dens-2fac-2 #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Densityplot für die Schlupfzeiten in [m] gruppiert nach den Floharten und den beiden Geschlechtern." ggplot(data = fac2_tbl, aes(x = hatch_time, fill = interaction(animal, sex))) + theme_minimal() + geom_density(alpha = 0.25) + labs(x = "Schlupfzeiten in [m]", y = "", fill = "Tierart & Geschlecht") + scale_fill_okabeito() + xlim(0, 2000) + geom_segment(aes(x = 0, y = 0, xend = 0, yend = 0.0017)) ``` ::: ## Visuelle Überprüfung text ::: panel-tabset ## Theoretisch Hier einmal die theoretischen Abbilungen von drei Gruppen mit Varianzhomogenität und Varianzheterogenität. Hier wird schnell klar, dass du mit mehr Behandlungsgruppen kaum die Annahme für Varianzhomogenität aufrecht erhalten kannst. Insbesondere Kontrollen oder Standardbehandlungen haben große oder kleinere Varianzen als der Rest der Behandlungen. Manchaml hat man aber Glück und die Varianzen sind eher homogen. ```{r} #| echo: false #| eval: true #| message: false #| warning: false #| label: fig-posthoc-var-homo-bar #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Darstellung der Varianzhomogenität und Varianzheterogenität in einem Barplot. *[Zum Vergrößern anklicken]*" p1 <- gg_template %+% tibble(x_fct = c("A", "B", "C"), mean_normal = c(1.5, 2.5, 5), sd_normal = c(0.5, 0.5, 0.5)) + aes(x_fct, mean_normal, fill = x_fct) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = mean_normal-(sd_normal), ymax = mean_normal+(sd_normal)), width = 0.2, linewidth = 1) + labs(x = "Behandlungsgruppen (X)", y = "Mittler Messwert (Y)", title = "Varianzhomogenität", subtitle = "Varianzen sind gleich") + ylim(0, 8) p2 <- gg_template %+% tibble(x_fct = c("A", "B", "C"), mean_normal = c(1.5, 2.5, 5), sd_normal = c(0.25, 1.25, 2.5)) + aes(x_fct, mean_normal, fill = x_fct) + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = mean_normal-(sd_normal), ymax = mean_normal+(sd_normal)), width = 0.2, linewidth = 1) + labs(x = "Behandlungsgruppen (X)", y = "Mittler Messwert (Y)", title = "Varianzheterogenität", subtitle = "Varianzen sind unterschiedlich") + ylim(0, 8) p1 + p2 + plot_layout(ncol = 2) + plot_annotation(tag_levels = 'A') + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` ```{r} #| echo: false #| eval: true #| message: false #| warning: false #| label: fig-posthoc-var-homo-box #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Darstellung der Varianzhomogenität und Varianzheterogenität in einem Boxplot. *[Zum Vergrößern anklicken]*" p1 <- gg_template %+% tibble(x_fct = gl(3, 5, labels = c("A", "B", "C")), y_val = c(0.25, 1, 1.5, 2, 2.75, 1.25, 2, 2.5, 3, 3.75, 3.75, 4.5, 5, 5.5, 6.25)) + aes(x_fct, y_val, fill = x_fct) + geom_boxplot(linewidth = 1) + labs(x = "Behandlungsgruppen (X)", y = "Mittler Messwert (Y)", title = "Varianzhomogenität", subtitle = "Varianzen sind gleich") + ylim(0, 8) p2 <- gg_template %+% tibble(x_fct = gl(3, 5, labels = c("A", "B", "C")), y_val = c(0.5, 1, 1.5, 2, 2.5, 0, 1, 2.5, 4, 6.5, 2, 3, 5, 6, 8)) + aes(x_fct, y_val, fill = x_fct) + geom_boxplot(linewidth = 1) + labs(x = "Behandlungsgruppen (X)", y = "Mittler Messwert (Y)", title = "Varianzheterogenität", subtitle = "Varianzen sind unterschiedlich") + ylim(0, 8) p1 + p2 + plot_layout(ncol = 2) + plot_annotation(tag_levels = 'A') + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` ## `geom_bar()` Wir bauen wir usn jetzt die Barplots oder Säulendiagramme? Entweder direkt aus `{emmeans}`, wie ich dir das weiter unten zum *Compact letter display* zeige oder eben als Wiederholung des [Kapitels zur Visualisierung von Daten](#sec-eda-ggplot). Wir bauen uns also einmal den einfaktoriellen Barplot wie auch den zweifaktoriellen Barplot. ### Einfaktorieller Barplot Für den einfaktoriellen Barplot brauchen wir einmal den Mittelwert und die Standardabweichung der einzelne Floharten. ```{r} #| message: false #| warning: false stat_fac1_tbl <- fac1_tbl |> group_by(animal) |> summarise(mean = mean(jump_length), sd = sd(jump_length)) stat_fac1_tbl ``` Dann können wir uns auch schon den Barplot mit `{ggplot}` erstellen. Später willst du dann noch das *Compact letter display* über den Balken als Mittelwerte ergänzen. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-posthoc-barplot-1fac #| fig-align: center #| fig-height: 4 #| fig-width: 4 #| fig-cap: "Beispielhafter einfaktorieller Barplot für die Sprungweiten in [cm] gruppiert nach den Floharten." ggplot(data = stat_fac1_tbl, aes(x = animal, y = mean, fill = animal)) + theme_minimal() + geom_bar(stat = "identity") + geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.2) + labs(x = "Flohart", y = "Sprungweite in [cm]") + theme(legend.position = "none") + scale_fill_okabeito() ``` ### Zweifaktorieller Barplot Für den zweifaktoriellen Barplot brauchen wir einmal den Mittelwert und die Standardabweichung der einzelne Floharten getrennt für die beiden Messorte. Das geht natürlich auch umgedreht, also die Messorte für die Floharten. Das kommt dann auf deine Fragestellung an. ```{r} #| message: false #| warning: false stat_fac2_tbl <- fac2_tbl |> group_by(animal, sex) |> summarise(mean = mean(jump_length), sd = sd(jump_length)) stat_fac2_tbl ``` Und dann können wir auch schon den zweifaktoriellen Barplot in `{ggplot}` erstellen. Du musst schauen, was du auf die x-Achse legst und was du dann auf die Legende und daher mit `fill` gruppierst. Damit die Positionen passen, spiele ich hier noch mit der Funktion `position_dodge()` rum. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-posthoc-barplot-2fac #| fig-align: center #| fig-height: 4 #| fig-width: 4 #| fig-cap: "Beispielhafter zweifaktorieller Barplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Messorten." ggplot(data = stat_fac2_tbl, aes(x = sex, y = mean, fill = animal)) + theme_minimal() + geom_bar(stat = "identity", width = 0.9, position = position_dodge(0.9)) + geom_errorbar(aes(ymin = mean-sd, ymax = mean+sd), width = 0.2, position = position_dodge(0.9)) + labs(x = "Flohart", y = "Sprungweite in [cm]", fill = "Tierart") + scale_fill_okabeito() ``` ## `geom_boxplot()` ### Einfaktorieller Boxplot Das praktische bei den Boxplots ist, dass wir hier nichts mehr vorrechnen müssen, sondern direkt die Boxplots in `{ggplot}` erstellen können. Ich finde man sieht immer in einem Boxplot besser, ob die Streuung um den Median eher homogen oder eher heterogen ist. Gerne ergänze ich noch den Mittelwert mit der Funktion `stat_summary()`. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-posthoc-boxplot-1fac #| fig-align: center #| fig-height: 4 #| fig-width: 4 #| fig-cap: "Beispielhafter einfaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten." ggplot(data = fac1_tbl, aes(x = animal, y = jump_length, fill = animal)) + theme_minimal() + geom_boxplot() + stat_summary(fun.y = mean, geom = "point", shape=23, size = 5, fill = "gray50") + labs(x = "Flohart", y = "Sprungweite in [cm]") + theme(legend.position = "none") + scale_fill_okabeito() ``` ### Zweifaktorieller Boxplot Den zweifaktoriellen Boxplot erstellen wir für die einzelnen Floharten getrennt für die beiden Messorte. Du musst schauen, was du auf die x-Achse legst und was du dann auf die Legende und daher mit `fill` gruppierst. Gerne ergänze ich noch den Mittelwert mit der Funktion `stat_summary()`, muss hier aber schauen, dass ich nach dem Faktor `animal` gruppiere und dann noch mit der Funktion `position_dodge()` die richtige Position finde. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-posthoc-boxplot-2fac #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Beispielhafter zweifaktorieller Boxplot für die Sprungweiten in [cm] gruppiert nach den Floharten und den beiden Messorten." ggplot(data = fac2_tbl, aes(x = sex, y = jump_length, fill = animal)) + theme_minimal() + geom_boxplot() + stat_summary(fun.y = mean, geom = "point", aes(group = animal), shape=23, size = 5, fill = "gray50", position = position_dodge(0.75)) + labs(x = "Flohart", y = "Sprungweite in [cm]", fill = "Tierart") + scale_fill_okabeito() ``` ::: ## Statistische Überprüfung ### Normalverteilung ::: panel-tabset ## Theoretisch [How to Calculate Skewness & Kurtosis in R](https://www.statology.org/skewness-kurtosis-in-r/) Die Schiefe (eng. *skewness*) ist ein Maß für die Asymmetrie einer Verteilung. Dieser Wert kann positiv oder negativ sein. Eine negative Schiefe deutet darauf hin, dass sich der Schwanz auf der linken Seite der Verteilung befindet, die sich in Richtung negativer Werte erstreckt. Eine positive Schiefe zeigt an, dass sich der Schwanz auf der rechten Seite der Verteilung befindet, die sich in Richtung positiver Werte erstreckt. Ein Wert von Null bedeutet, dass die Verteilung überhaupt nicht schief ist, d. h. die Verteilung ist vollkommen symmetrisch. Die Kurtosis (eng. *kurtosis*) ist ein Maß dafür, ob eine Verteilung im Vergleich zu einer Normalverteilung ein starkes oder schwaches Schwanzende aufweist. Die Kurtosis einer Normalverteilung beträgt 3. Wenn eine gegebene Verteilung eine Kurtosis von weniger als 3 aufweist, wird sie als playkurtisch bezeichnet, was bedeutet, dass sie dazu neigt, weniger und weniger extreme Ausreißer zu produzieren als die Normalverteilung. Wenn eine gegebene Verteilung eine Kurtosis größer als 3 hat, wird sie als leptokurtisch bezeichnet, was bedeutet, dass sie dazu neigt, mehr Ausreißer als die Normalverteilung zu produzieren. In einigen Formeln (Definition von Fisher) wird von der Kurtosis 3 abgezogen, um den Vergleich mit der Normalverteilung zu erleichtern. Nach dieser Definition hätte eine Verteilung eine größere Kurtosis als eine Normalverteilung, wenn der Kurtosis-Wert größer als 0 wäre. ```{r} #| echo: false #| message: false #| warning: false #| label: fig-pretest-rlnorm-01 #| fig-align: center #| fig-height: 4 #| fig-width: 8 #| fig-cap: "Die t-Verteilung für drei beispielhafte Freiheitsgrade. Je größer die Freiheitsgrade und damit die Fallzahl, desto näher kommt die t-Verteilung einer Normalverteilung nahe. *[Zum Vergrößern anklicken]*" p1 <- ggplot(data.frame(x = c(-3.25, 3.25)), aes(x)) + theme_minimal() + geom_hline(yintercept = 0, gray = "gray50") + geom_segment(aes(x = 0, xend = 0, y = 0, yend = 0.4)) + stat_function(fun = dnorm, linewidth = 1, args = list(mean = 0, sd = 1), xlim = c(-4, 4), color = "#E89F00") + scale_x_continuous(breaks = 0, labels = "Mittelwert = Median = Modus", guide = guide_axis(n.dodge = 2)) + labs(x = "", y = "") + theme(legend.text = element_text(size = 14), legend.title = element_text(size = 14, face = 2), axis.text.y = element_blank(), axis.title.x = element_text(size = 16, face = 2), axis.text.x = element_text(size = 12), plot.title = element_text(size = 17), plot.subtitle = element_text(size = 12, face = "italic"), plot.caption = element_text(size = 12), panel.grid.major.y = element_blank(), panel.grid.major.x = element_blank(), panel.grid.minor.y = element_blank(), legend.position = "top") p2 <- ggplot(data.frame(x = c(-3.25, 3.25)), aes(x)) + theme_minimal() + geom_hline(yintercept = 0, gray = "gray50") + geom_segment(aes(x = 0.8, xend = 0.8, y = 0, yend = 2.46)) + geom_segment(aes(x = 0.625, xend = 0.625, y = 0, yend = 1.71), linetype = 11) + geom_segment(aes(x = 0.45, xend = 0.45, y = 0, yend = 0.68)) + stat_function(fun = dbeta, linewidth = 1, args = list(shape1 = 5, shape2 = 2), xlim = c(0.15, 1), color = "#E89F00") + labs(x = "", y = "", title = "Linksschief", subtitle = "Negative Schiefe") + scale_x_continuous(breaks = c(0.45, 0.625, 0.8), labels = c("Mittelwert", "Median", "Modus"), guide = guide_axis(n.dodge = 2)) + theme(legend.text = element_text(size = 14), legend.title = element_text(size = 14, face = 2), axis.text.y = element_blank(), axis.title.x = element_text(size = 16, face = 2), axis.text.x = element_text(size = 12), plot.title = element_text(size = 17), plot.subtitle = element_text(size = 12, face = "italic"), plot.caption = element_text(size = 12), panel.grid.major.y = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), panel.grid.minor.y = element_blank(), legend.position = "top") p3 <- ggplot(data.frame(x = c(-3.25, 3.25)), aes(x)) + theme_minimal() + geom_hline(yintercept = 0, gray = "gray50") + geom_segment(aes(x = 2, xend = 2, y = 0, yend = 0.183)) + geom_segment(aes(x = 5, xend = 5, y = 0, yend = 0.102), linetype = 11) + geom_segment(aes(x = 10, xend = 10, y = 0, yend = 0.017)) + stat_function(fun = dchisq, linewidth = 1, args = list(df = 4), xlim = c(0, 16), color = "#E89F00") + scale_x_continuous(breaks = c(2, 5, 10), labels = c("Modus", "Median", "Mittelwert"), guide = guide_axis(n.dodge = 2)) + labs(x = "", y = "", title = "Rechtsschief", subtitle = "Positive Schiefe") + theme(legend.text = element_text(size = 14), legend.title = element_text(size = 14, face = 2), axis.text.y = element_blank(), axis.title.x = element_text(size = 16, face = 2), axis.text.x = element_text(size = 12), plot.title = element_text(size = 17), plot.subtitle = element_text(size = 12, face = "italic"), plot.caption = element_text(size = 12), panel.grid.major.y = element_blank(), panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank(), panel.grid.minor.y = element_blank(), legend.position = "top") p2 + p1 + p3 + plot_layout(ncol = 3) + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` ## `{moments}` ```{r} pacman::p_load(moments, broom, scales, report) #calculate skewness skewness(fac2_tbl$hatch_time) ``` ```{r} #calculate kurtosis kurtosis(fac2_tbl$hatch_time) ``` ```{r} fac2_tbl |> group_by(animal, sex) |> summarise(kurtosis = kurtosis(hatch_time), skewness = skewness(hatch_time)) ``` Die Momente-Bibliothek bietet auch die Funktion `jarque.test()`, die einen Anpassungsgütetest durchführt, der feststellt, ob die Stichprobendaten eine Schiefe und eine Wölbung aufweisen, die einer Normalverteilung entsprechen oder nicht. Die Null- und Alternativhypothesen dieses Tests lauten wie folgt: Nullhypothese: Der Datensatz weist eine Schiefe und Wölbung auf, die einer Normalverteilung entspricht. Alternativhypothese: Der Datensatz weist eine Schiefe und eine Kurtosis auf, die nicht mit einer Normalverteilung übereinstimmen. ```{r} jarque.test(fac2_tbl$hatch_time) ``` ```{r} fac2_tbl |> split(~ animal + sex) |> map(~jarque.test(.x$hatch_time)) |> map(tidy) |> bind_rows(.id = "test") |> select(test, p.value) |> mutate(decision = ifelse(p.value <= 0.05, "reject normal", "normal"), p.value = pvalue(p.value, accuracy = 0.001)) ``` ## `{report}` ```{r} fac2_tbl |> select(jump_length, hatch_time) |> report() ``` ::: text ::::: panel-tabset ## Theoretisch ## `{stats}` `shapiro.test()` ## `{performance}` Wenn du überprüfen willst, ob dein Messwert $y$ einer Normalverteilung folgt, dann kannst du die Funktion `check_normality()` nutzt. Die Funktion rechnet dann den [Shapiro-Wilk-Test](https://de.wikipedia.org/wiki/Shapiro-Wilk-Test) um auf eine Abweichung von der Normalverteilung zu testen. Hierzu ist anzumerken, dass der Test relativ empfindlich bei Abweichungen in den Verteilungsschwänzen ist. Darüber hinaus braucht der Shapiro-Wilk-Test auch etwas Fallzahl, damit er auf die Normalverteilung testen kann. Im Folgenden schauen wir uns den Code für ein einfaktorielles, zweifaktorielles und dreifaktorielles Modell einmal an. Am Ende des Abschnitts gehe ich nochmal darauf ein, was du machen kannst, wenn du keine Normalverteilung in deinem Messwert $y$ vorliegen hast. #### Einfaktoriell {.unnumbered .unlisted} Beginnen wir wieder mit einem einfaktoriellen Modell. Wir wollen wissen, ob unsere Sprungweite in \[cm\] über unsere verschiedenen Floharten normalverteilt ist. Wir bauen also erstmal das Modell und schicken es dann in die Funktion `check_normality()` aus dem R Paket `{performance}`. ```{r} lm(jump_length ~ animal, data = fac1_tbl) |> check_normality() ``` Die Funktion sagt, dass wir eine Normalverteilung in unseren Daten vorliegen haben. Wir können uns auch einen Diagnoseplot wiedergeben lassen. Dafür müssen wir die Funktion nur an die Funktion `plot()` weiterleiten. Das Schöne ist, dass die Abbildung uns auch gleich sagt, was wir zu erwarten haben um eine Normalverteilung anzunehmen. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-normal-f1 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem einfaktoriellen Modell." lm(jump_length ~ animal, data = fac1_tbl) |> check_normality() |> plot() + scale_fill_okabeito() ``` #### Zweifaktoriell {.unnumbered .unlisted} Im zweifaktoriellen Fall ändert sich jetzt nur das Model. Wir haben eben zwei Faktoren vorliegen und diese müssen wir dann mit ins Modell nehmen. Ich habe hier auch gleich den Interaktionsterm mit ergänzt, ich teste gerne das Modell, was ich dann später auch auswerten möchte. ::: panel-tabset ## Sprungweite `jump_length` ```{r} lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> check_normality() ``` ## Schlupfzeit `hatch_time` ```{r} lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> check_normality() ``` ::: Auch hier haben wir eine Normalverteilung in den Daten vorliegen. Gerne schaue ich mir auch die Abbildung der Residuen einmal an und das geht dann flott über die Funktion `plot()`. Da musst du nur die Ausgabe der Funktion `check_normality()` weiterleiten. Die leichten Bögen in den Punkten kommen von den unterschiedlichen Faktoren und deren Effekten auf die Sprungweite. ::: panel-tabset ## Sprungweite `jump_length` ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-normal-f2-1 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem zweifaktoriellen Modell." lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> check_normality() |> plot() + scale_fill_okabeito() ``` ## Schlupfzeit `hatch_time` ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-normal-f2-2 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_normality()` zur Überprüfung der Normalverteilung des Messwerts in einem zweifaktoriellen Modell." lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> check_normality() |> plot() + scale_fill_okabeito() ``` ::: ::::: ### Varianzhomogenität ::::: panel-tabset ## Theoretisch ## `{stats}` `bartlett.test()` ## `{car}` `leveneTest()` ## `{performance}` `check_homogeneity()` Wir können auch die Funktion `check_homogeneity()` aus dem Paket `{performance}` nutzen. Wir erhalten hier auch gleich eine Entscheidung in englischer Sprache ausgegeben. Die Funktion `check_homogeneity()` nutzt den Bartlett-Test. Wir können in Funktion auch andere Methoden mit `method = c("bartlett", "fligner", "levene", "auto")` wählen. Wir sehen, dass sich die Implementierung des Bartlett-Tests in `check_homogeneity()` nicht von der Funktion `bartlett.test()` unterscheidet, aber die Entscheidung gegen die Varianzhomogenität zu einem Signifikanzniveau von 5% gefällt wird. Nicht immer hilft einem der Entscheidungtext einer Funktion. Nachdem wir die Normalverteilung des Messwertes $y$ getestet haben, wollen wir jetzt noch schauen, ob unsere Faktoren einer Varianzhomogenität über die Gruppen genügt. Die Funktion `check_homogeneity()` nutzt den [Bartlett-Test](https://de.wikipedia.org/wiki/Bartlett-Test) um auf eine Abweichung von der Varianzhomogenität zu testen. Es gibt noch verschiedene andere Tests, aber ich würde hier empfehlen bei eiem Test zu bleiben und dann mit dem Ergebnis zu leben. Wild verschiedene Tests auf Varianzhomogenität durchzuführen führt dann auch nur zu einem wilden Fruchtsalat an Ergebnissen mit denen keiner was anfangen kann. Im Folgenden schauen wir uns den Code für ein einfaktorielles, zweifaktorielles und dreifaktorielles Modell einmal an. Am Ende des Abschnitts gehe ich nochmal darauf ein, was du machen kannst, wenn du keine Varianzhomogenität in deinen Faktoren vorliegen hast. #### Einfaktoriell {.unnumbered .unlisted} Beginnen wir wieder mit dem einfaktoriellen Modell. Wir stecken das Modell dann einfach in die Funktion `check_homogeneity()` und erhalten die Information über die Varianzhomogenität wiedergegeben. ```{r} lm(jump_length ~ animal, data = fac1_tbl) |> check_homogeneity() ``` Wunderbar, wir haben keine Abweichung von der Varianzhomongenität. Wir können uns auch die Daten nochmal anschauen. Hier sehen wir aber schon, dass die Daten etwas heterogen *aussehen* der Test aber die Homogenität nicht ablehnt. Das ist immer schwierig bei kleinen Fallzahlen. Aber an irgendwas muss man eben glauben. ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-variance-f1 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem einfaktoriellen Modell." lm(jump_length ~ animal, data = fac1_tbl) |> check_homogeneity() |> plot() + scale_fill_okabeito() ``` #### Zweifaktoriell {.unnumbered .unlisted} Jetzt können wir mal schauen, was passiert wenn wir die Anzahl an möglichen Faktorkombinationen erhöhen indem wir ein zweifaktorielles Modell nutzen. Hier haben wir dann ja sechs Faktorkombinationen oder Gruppen die dann alle homogen in den Varianzen sein müssen. ::: panel-tabset ## Sprungweite `jump_length` ```{r} lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> check_homogeneity() ``` ## Schlupfzeit `hatch_time` ```{r} lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> check_normality() ``` ::: Wir haben auch hier Varianzhomogenität über alle Gruppen der Faktoren vorliegen. Wenn du dir jetzt die Abbildung zu dem Test anschaust, dann siehst du auch hier, dass die Violinplots eben dann doch alle etas anders aussehen. Wir haben aber hier auch das gleiche Problem wie bei dem einfaktoriellen Fall, wir haben eben dann doch recht wenig Fallzahl in unseren Daten. ::: panel-tabset ## Sprungweite `jump_length` ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-variance-f2-1 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem zweifaktoriellen Modell." lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) |> check_homogeneity() |> plot() + scale_fill_okabeito() ``` ## Schlupfzeit `hatch_time` ```{r} #| message: false #| echo: true #| warning: false #| label: fig-ggplot-check-anova-variance-f2-2 #| fig-align: center #| fig-height: 4 #| fig-width: 6 #| fig-cap: "Schnelle Abbildung der Residuen aus `check_homogeneity()` zur Überprüfung der Varianzhomogenität der Faktoren in einem zweifaktoriellen Modell." lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) |> check_homogeneity() |> plot() + scale_fill_okabeito() ``` ::: ::::: ## Auswege #### Flowchart {.unnumbered .unlisted} text ```{mermaid} %%| label: fig-mermaid-pretest-01 %%| fig-width: 6 %%| fig-cap: "Flowchart für die Entscheidung welches statistische Modell gerechnet werden kann: `lm()`, lineare Regression, `gls()`, generalized least squares Regression, `rq()`, Quantilesregression. Die Funktionen finden sich teilweise in eigenen R Paketen. Bei einem nicht-normalverteilten Outcome mit Varianzheterogenität über die Gruppen müssen wir nochmal gemeinsam in die Daten schauen." flowchart TD A("Normalverteiltes Outcome in jeder Versuchsgruppe"):::factor --- B(((ja))) --> B1 A("Normalverteiltes Outcome in jeder Versuchsgruppe"):::factor --- F(((nein))) --> B2 subgraph B1["Mittelwertsvergleiche"] C("Varianzhomogenität über alle Gruppen"):::factor --- D((("ja"))) --> E("lm()"):::factor C("Varianzhomogenität über alle Gruppen"):::factor --- J(((nein))) --> K("gls()"):::factor end subgraph B2["Medianvergleiche oder Anteile"] G("Varianzhomogenität über alle Gruppen"):::factor --- H(((ja))) --> I("rq()"):::factor G("Varianzhomogenität über alle Gruppen"):::factor --- L(((nein))) end classDef factor fill:#56B4E9,stroke:#333,stroke-width:0.75px ``` Du findest die vorgeschlagenen Funktionen dann in den entsprechenden Kapiteln zur ANOVA und den Posthoc-Tests. Du kannst dir dann dort den Code zusammensuchen. Je nach deiner Datenlage musst du dann nochmal etwas an dem R Code programmieren. Beachte, dass die Funktionen sich teilweise in eigenen R Paketen finden lassen. So ist die Funktion `gls()` im R Paket `{nlme}` und die Funktion `rq()` im R Paket `{quantreg}` zu finden. Du kannst auch bei Varianzheterogenität das R Paket `{sandwich}` nutzen und einen entsprechend angepassten Varianzschätzer. Mehr findest du dazu bei den Posthoc-Test in dem Abschnitt zu dem [Gruppenvergleich unter Varianzheterogenität](#sec-posthoc-var-heterogen) oder gleich im ersten Zerforschenbeispiel zum einfaktoriellen Barplot. ### Transformation ::: panel-tabset ## `log()` ```{r} log_tbl <- fac2_tbl |> mutate(log_hatch_time = log(hatch_time)) ``` ```{r} #| message: false #| echo: false #| label: fig-log-pretest-1 #| fig-align: center #| fig-height: 3.5 #| fig-width: 8 #| fig-cap: "Histogramm der nicht transfomierten und transformierten Daten. **(A)** Nicht transformierte, rohe Daten. **(B)** $log$-transformierte Daten. *[Zum Vergrößern anklicken]*" #| layout-nrow: 1 p1 <- ggplot(log_tbl, aes(hatch_time)) + theme_minimal() + geom_density(fill = cbbPalette[2], color = "black") + labs(x = "Zeit bis zum Schlüpfen in [h]", y = "") + xlim(-100, NA) p2 <- ggplot(log_tbl, aes(log_hatch_time)) + theme_minimal() + geom_density(fill = cbbPalette[3], color = "black") + labs(x = "Zeit bis zum Schlüpfen in log(h)", y = "") p1 + p2 + theme(panel.grid.minor.x = element_blank()) + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` ## `rank()` ```{r} rank_tbl <- fac2_tbl |> mutate(rank_hatch_time = rank(hatch_time)) ``` ```{r} #| message: false #| echo: false #| label: fig-log-pretest-2 #| fig-align: center #| fig-height: 3.5 #| fig-width: 8 #| fig-cap: "Histogramm der nicht transfomierten und transformierten Daten. **(A)** Nicht transformierte, rohe Daten. **(B)** $log$-transformierte Daten. *[Zum Vergrößern anklicken]*" #| layout-nrow: 1 p1 <- ggplot(rank_tbl, aes(hatch_time)) + theme_minimal() + geom_density(fill = cbbPalette[2], color = "black") + labs(x = "Zeit bis zum Schlüpfen in [h]", y = "") + xlim(-100, NA) p2 <- ggplot(rank_tbl, aes(rank_hatch_time)) + theme_minimal() + geom_density(fill = cbbPalette[3], color = "black") + labs(x = "Zeit bis zum Schlüpfen in log(h)", y = "") p1 + p2 + theme(panel.grid.minor.x = element_blank()) + plot_annotation(tag_levels = 'A', tag_prefix = '(', tag_suffix = ')') & theme(plot.tag = element_text(size = 16, face = "bold")) ``` ## `boxcox()` ```{r} jump_mod <- lm(hatch_time ~ animal + sex + animal:sex, data = fac2_tbl) ``` ```{r} bc_obj <- boxcox(jump_mod, plotit = FALSE) lambda <- bc_obj$x[which.max(bc_obj$y)] lambda ``` ```{r} boxcox_tbl <- fac2_tbl |> mutate(jump_boxcox = ((jump_length ^ lambda - 1)/lambda)) ``` ```{r} #| message: false #| echo: true #| label: fig-boxcox-pretest #| fig-align: center #| fig-height: 5 #| fig-width: 5 #| fig-cap: "Histogramm der Box-Cox transfomierten Daten in `ggplot` dargestellt." ggplot(data = boxcox_tbl) + theme_minimal() + geom_histogram(aes(x = jump_boxcox), alpha = 0.9, fill = cbbPalette[3], color = "black") + labs(x = 'Box-Cox transformierte Sprungweiten', y = 'Anzahl') ``` R Paket `{MASS}` ::: ### Modellierung #### R Paket `{performance}` {.unnumbered .unlisted} Wenn du noch etwas weiter gehen möchtest, dann kannst du dir noch die Hilfeseite von dem R Paket `{performance}` [Robust Estimation of Standard Errors, Confidence Intervals, and p-values](https://easystats.github.io/parameters/articles/model_parameters_robust.html?q=Heteroskedasticity#robust-covariance-matrix-estimation-from-model-parameters) anschauen. Die Idee ist hier, dass wir die Varianz/Kovarianz robuster daher mit der Berücksichtigung von Varianzheterogenität (eng. *heteroskedasticity*) schätzen. ```{r} fit_fac2 <- lm(jump_length ~ animal + sex + animal:sex, data = fac2_tbl) ``` ```{r} fit_fac2 |> model_parameters() ``` ```{r} fit_fac2 |> model_parameters(vcov = "HC3") ``` #### R Paket `{emmeans}` {.unnumbered .unlisted} ```{r} fit_fac2 |> emmeans(~ animal * sex, vcov. = sandwich::vcovHAC) ``` #### R Paket `{nlme}` {.unnumbered .unlisted} ```{r} fit_fac2 <- gls(jump_length ~ animal + sex + animal:sex, data = fac2_tbl, weights = varIdent(form = ~ 1 | animal*sex)) ``` ```{r} fit_fac2 |> model_parameters() ``` ## Referenzen {.unnumbered}

81.1 Allgemeiner Hintergrund

Das Modell

Gibt es noch mehr Vortests?

81.1.1 Normalverteilung

81.1.2 Varianzhomogenität

81.2 Genutzte R Pakete

81.3 Daten

Einfaktorieller Datensatz

Zweifaktorieller Datensatz

81.4 Visuelle Überprüfung

81.4.1 Einfaktorieller Barplot

81.4.2 Zweifaktorieller Barplot

81.4.3 Einfaktorieller Boxplot

81.4.4 Zweifaktorieller Boxplot

81.5 Statistische Überprüfung

81.5.1 Normalverteilung

Einfaktoriell

Zweifaktoriell

81.5.2 Varianzhomogenität

Einfaktoriell

Zweifaktoriell

81.6 Auswege

Flowchart

81.6.1 Transformation

81.6.2 Modellierung

R Paket {performance}

R Paket {emmeans}

R Paket {nlme}

Referenzen

R Paket `{performance}`

R Paket `{emmeans}`

R Paket `{nlme}`