pacman::p_load(tidyverse, magrittr, broom, quantreg,
see, performance, emmeans, multcomp, janitor,
parameters, conflicted)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")30 Die ANCOVA
Letzte Änderung am 02. April 2024 um 09:52:43
“It’s better to solve the right problem approximately than to solve the wrong problem exactly.” — John Tukey
Eine Kovariate ist eine Variable, die mit berücksichtigt wird, um mögliche verzerrende Einflüsse auf die Analyseergebnisse (Konfundierung) abzuschätzen oder zu verringern.
Eigentlich hat sich die Analysis of Covariance (ANCOVA) etwas überlebt. Wir können mit dem statistischen Modellieren eigentlich alles was die ANCOVA kann plus wir erhalten auch noch Effektschätzer für die Kovariaten und die Faktoren. Dennoch hat die ANCOVA ihren Platz in der Auswertung von Daten. Wenn du ein oder zwei Faktoren hast plus eine numerische Variable, wie das Startgewicht, für die du die Analyse adjustieren möchtest, dann ist die ANCOVA für dich gemacht.
Also kurz gesprochen adjustiert die Analysis of Covariance (ANCOVA) die Faktoren einer ANOVA um eine kontinuierliche Covariate. Adjustiert bedeutet in dem Fall, dass die Effekte des unterschiedlichen Startgewichts von Pflanzen durch das Einbringen der Kovariate mit in der statistischen Analyse berücksichtigt werden. Wir werden hier auch nur über die Nutzung in R sprechen und auf die theoretische Herleitung verzichten.
Wir können die einfaktorielle ANCOVA in folgender Form schreiben. Wir haben haben einen Faktor \(x_1\) und eine Kovariate oder aber ein numerisches \(x_2\). Damit sähe die ANCOVA wie folgt aus.
\[ y \sim x_1 + x_2 \]
Damit ist die ANCOVA aber sehr abstrakt beschrieben. Der eine Faktor kommt damit gar nicht zur Geltung. Deshalb schreiben wir die ANCOVA wie folgt mit einem \(f_1\) für den Faktor und einem \(c_1\) für eine numerische Kovariate. Damit haben wir einen bessere Übersicht.
\[ y \sim f_1 + c_1 \]
Somit erklärt sich die zweifaktorielle ANCOVA schon fast von alleine. Wir erweitern einfach das Modell um einen zweiten Faktor \(f_2\) und haben somit eine zweifaktorielle ANCOVA.
\[ y \sim f_1 + f_2 + c_1 \]
Im Folgenden schauen wir uns einmal die Daten und die Hypothesen zu einer möglichen Fragestellung an.
30.1 Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
30.2 Daten
Für unser Beispiel nutzen wir die Daten der Sprungweite in [cm] von Flöhen auf Hunde-, Katzen- und Füchsen. Damit haben wir den ersten Faktor animal mit drei Leveln. Als Kovariate schauen wir uns das Gewicht als numerische Variable an. Schlussendlich brauchen wir noch das Outcome jump_length als \(y\). Für die zweifaktorielle ANCOVA nehmen wir noch den Faktor sex mit zwei Leveln hinzu.
ancova_tbl <- read_csv2("data/flea_dog_cat_length_weight.csv") |>
select(animal, sex, jump_length, weight) |>
mutate(animal = as_factor(animal))In der Tabelle 56.1 ist der Datensatz ancova_tbl nochmal dargestellt.
| animal | sex | jump_length | weight |
|---|---|---|---|
| cat | male | 15.79 | 6.02 |
| cat | male | 18.33 | 5.99 |
| cat | male | 17.58 | 8.05 |
| cat | male | 14.09 | 6.71 |
| cat | male | 18.22 | 6.19 |
| cat | male | 13.49 | 8.18 |
| … | … | … | … |
| fox | female | 27.81 | 8.04 |
| fox | female | 24.02 | 9.03 |
| fox | female | 24.53 | 7.42 |
| fox | female | 24.35 | 9.26 |
| fox | female | 24.36 | 8.85 |
| fox | female | 22.13 | 7.89 |
Unser zweiter Datensatz ist ein Anwendungsdatensatz aus dem Gemüsebau. Wir schauen uns das Wachstum von drei Gurkensorten über siebzehn Wochen an. Die Gurkensorten sind hier unsere Versuchsgruppen. Da wir es hier mit echten Daten zu tun haben, müssen wir uns etwas strecken damit die Daten dann auch passen. Wir wollen das Wachstum der drei Gurkensorten über die Zeit betrachten - also faktisch den Verlauf des Wachstums.
gurke_raw_tbl <- read_excel("data/wachstum_gurke.xlsx") |>
clean_names() |>
select(-pfl, -erntegewicht) |>
mutate(versuchsgruppe = as_factor(versuchsgruppe)) In der Tabelle 56.2 sehen wir einmal die rohen Daten dargestellt.
| versuchsgruppe | t1 | t2 | t3 | t4 | t5 | t6 | t7 | t8 | t9 | t10 | t11 | t12 | t13 | t14 | t15 | t16 | t17 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Proloog L | 5.5 | 6.1 | 7.4 | 8.9 | 9.9 | 12 | 14.4 | 17 | 19.8 | 21.2 | 23.2 | 24 | 29.7 | 32.8 | NA | NA | NA |
| Proloog L | 4.6 | 5.1 | 6.4 | 5.7 | 5.5 | 5.2 | 5 | 5 | 4.5 | 0 | 0 | 0 | 0 | 0 | NA | NA | NA |
| Proloog L | 5.3 | 5.8 | 6.8 | 8.3 | 9 | 10 | 12.3 | 14.6 | 17.6 | 19.3 | 23.1 | 23.8 | 31.7 | 32.3 | NA | NA | NA |
| Proloog L | 5.4 | 5.7 | 6.9 | 8.2 | 8.6 | 10 | 12.1 | 14.5 | 16.2 | 17.1 | 19.3 | 21.6 | 28.5 | 30 | NA | NA | NA |
| Proloog L | 5 | 5.5 | 6.3 | 7.5 | 8.3 | 10 | 12.2 | 14.4 | 16.5 | 19.9 | 21 | 22.9 | 30.4 | 31 | NA | NA | NA |
| Proloog L | 4.2 | 4.6 | 5.4 | 5.5 | 5.2 | 5.3 | 6.1 | 6.5 | 8 | 9.3 | 11 | 12.5 | 22.3 | 24.2 | NA | NA | NA |
| … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … | … |
| Katrina D | 0.62 | 0.64 | 0.65 | 0.65 | 0.74 | 0.9 | 1.07 | 1.09 | 1.18 | 1.09 | 1.13 | 1.14 | 1.3 | 1.48 | 1.75 | 2.02 | 2.97 |
| Katrina D | 0.55 | 0.59 | 0.64 | 0.62 | 0.64 | 0.7 | 0.68 | 0.66 | 0.66 | 0.7 | 0.6 | 0.6 | 0.57 | 0 | NA | NA | NA |
| Katrina D | 0.64 | 0.73 | 0.77 | 0.8 | 0.86 | 1.2 | 1.25 | 1.42 | 1.75 | 1.92 | 2.24 | 2.52 | 4.02 | 3.9 | NA | NA | NA |
| Katrina D | 0.86 | 0.87 | 0.9 | 1.1 | 1.12 | 1.4 | 1.7 | 2.1 | 2.35 | 2.4 | 2.51 | 2.57 | 4.07 | 4.3 | NA | NA | NA |
| Katrina D | 0.62 | 0.67 | 0.7 | 0.71 | 0.78 | 0.95 | 1.1 | 1.24 | 1.41 | 1.66 | 2.02 | 2.25 | 3.85 | 4.07 | NA | NA | NA |
| Tageslänge | 13.98 | 14.05 | 14.12 | 14.18 | 14.25 | 14.3 | 14.37 | 14.43 | 14.5 | 14.57 | 14.62 | 14.68 | 14.75 | 14.8 | 15.1 | 15.15 | 15.32 |
Wir haben zwei Typen von Daten für das Gurkenwachstum. Einmal messen wir den Durchmesser für jede Sorte (D im Namen der Versuchsgruppe) oder aber die Länge (L im Namen der Versuchsgruppe). Wir betrachten hier nur das Längenwachstum und deshalb filtern wir erstmal nach allen Versuchsgruppen mit einem L im Namen. Dann müssen wir die Daten noch in Long-Format bringen. Da wir dann auch noch auf zwei Arten die Daten über die Zeit darstellen wollen, brauchen wir einmal die Zeit als Faktor time_fct und einmal als numerisch time_num. Leider haben wir auch Gurken mit einer Länge von 0 cm. Diese Gurken schmeißen wir am Ende mal raus. Auch haben wir ab Woche 14 keine Messungen mehr in der Versuchsgruppe Prolong, also nehmen wir auch nur die Daten bis zur vierzehnten Woche.
gurke_time_len_tbl <- gurke_raw_tbl |>
filter(str_detect(versuchsgruppe, "L$")) |>
mutate(versuchsgruppe = factor(versuchsgruppe,
labels = c("Proloog", "Quarto", "Katrina"))) |>
pivot_longer(cols = t1:t17,
values_to = "length",
names_to = "time") |>
mutate(time_fct = as_factor(time),
time_num = as.numeric(time_fct)) |>
filter(length != 0) |>
filter(time_num <= 14)30.3 Hypothesen für die ANCOVA
Wir haben für jeden Faktor der ANCOVA ein Hypothesenpaar sowie ein Hypothesenpaar für die Kovariate. Im Folgenden sehen wir die jeweiligen Hypothesenpaare.
Einmal für animal, als Haupteffekt. Wir nennen einen Faktor den Hauptfaktor, weil wir an diesem Faktor am meisten interessiert sind. Wenn wir später einen Posthoc Test durchführen würden, dann würden wir diesen Faktor nehmen. Wir sind primär an dem Unterschied der Sprungweiten in [cm] in Gruppen Hund, Katze und Fuchs interessiert.
\[ \begin{aligned} H_0: &\; \bar{y}_{cat} = \bar{y}_{dog} = \bar{y}_{fox}\\ H_A: &\; \bar{y}_{cat} \ne \bar{y}_{dog}\\ \phantom{H_A:} &\; \bar{y}_{cat} \ne \bar{y}_{fox}\\ \phantom{H_A:} &\; \bar{y}_{dog} \ne \bar{y}_{fox}\\ \phantom{H_A:} &\; \mbox{für mindestens ein Paar} \end{aligned} \]
Du kannst mehr über Geraden sowei lineare Modelle und deren Eigenschaften im Kapitel 38 erfahren.
Für die Kovariate testen wir anders. Die Kovariate ist ja eine numerische Variable. Daher ist die Frage, wann gibt es keinen Effekt von weight auf die Sprunglänge? Wenn wir eine parallele Linie hätten. Das heißt, wenn sich der Wert von weight ändert, ändert sich der Wert von jump_length nicht. Wir schreiben, dass sich die Steigung der Geraden nicht ändert. Wir bezeichnen die Steigung einer Graden mit \(\beta\). Wenn kein Effekt vorliegt und die Nullhpyothese gilt, dann ist die Steigung der Geraden \(\beta_{weight} = 0\).
\[ \begin{aligned} H_0: &\; \beta_{weight} = 0\\ H_A: &\; \beta_{weight} \neq 0 \end{aligned} \]
Du kannst dir überlegen, ob due die Interaktion zwischen dem Faktor und der Kovariate mit ins Modell nehmen willst. Eigentlich schauen wir uns immer nur die Interaktion zwischen den Faktoren an. Generell schreiben wir eine Interaktionshypothese immer in Prosa.
\[ \begin{aligned} H_0: &\; \mbox{keine Interaktion}\\ H_A: &\; \mbox{eine Interaktion zwischen animal und site} \end{aligned} \]
Wir haben also jetzt die verschiedenen Hypothesenpaare definiert und schauen uns jetzt die ANCOVA in R einmal in der Anwendung an.
30.4 Die einfaktorielle ANCOVA in R
Du kannst mehr über Geraden sowie lineare Modelle und deren Eigenschaften im Kapitel 38 erfahren.
Wir können die ANCOVA ganz klassisch mit dem linaren Modell fitten. Wir nutzen die Funktion lm() um die Koeffizienten des linearen Modellls zu erhalten. Wir erinnern uns, wir haben haben einen Faktor \(f_1\) und eine Kovariate bezwiehungsweise ein numerisches \(c_1\). In unserem Beispiel sieht dann der Fit des Modells wie folgt aus.
fit_1 <- lm(jump_length ~ animal + weight + animal:weight, data = ancova_tbl)Nachdem wir das Modell in dem Objekt fit_1 gespeichert haben können wir dann das Modell in die Funktion anova() pipen. Die Funktion erkennt, das wir eine ANCOVA rechnen wollen, da wir in unserem Modell einen Faktor und eine Kovariate mit enthalten haben.
fit_1 |> anova ()Analysis of Variance Table
Response: jump_length
Df Sum Sq Mean Sq F value Pr(>F)
animal 2 2693.8 1346.88 204.2764 <2e-16 ***
weight 1 1918.0 1917.99 290.8961 <2e-16 ***
animal:weight 2 0.3 0.14 0.0214 0.9788
Residuals 594 3916.5 6.59
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In der ANCOVA erkennne wir nun, dass der Faktor animal signifikant ist. Der \(p\)-Wert ist mit \(<0.001\) kleiner das das Signifikanzniveau \(\alpha\) von 5%. Ebenso ist die Kovariate weight signifikant. Der \(p\)-Wert ist ebenfalls mit \(<0.001\) kleiner das das Signifikanzniveau \(\alpha\) von 5%. Wir können also schlussfolgern, dass sich mindestens eine Gruppenvergleich der Level des Faktors animal voneinander unterscheidet. Wir wissen auch, dass mit der Zunahme des Gewichts, die Sprunglänge sich ändert.
Die ANCOVA liefert keine Informationen zu der Größe oder der Richtung des Effekts der Kovariate.
Was wir nicht wissen, ist die Richtung. Wir wissen nicht, ob mit ansteigenden Gewicht sich die Sprunglänge erhöht oder vermindert. Ebenso wenig wissen wir etwas über den Betrag des Effekts. Wieviel weiter springen denn nun Flöhe mit 1 mg Gewicht mehr? Wir haben aber die Möglichkeit, den Sachverhalt uns einmal in einer Abbildung zu visualisieren. In Abbildung 30.1 sehen wir die Daten einmal als Scatterplot dargestellt.
ggplot(ancova_tbl, aes(weight, jump_length, color = animal)) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_okabeito() +
theme_minimal() +
geom_point() +
labs(color = "Tierart", shape = "Geschlecht") 
Der Abbildung 30.1 können wir jetzt die positive Steigung entnehmen sowie die Reihenfolge der Tierarten nach Sprungweiten. Die ANCOVA sollte immer visualisiert werden, da sich hier die Stärke der Methode mit der Visualiserung verbindet.
30.5 Die zweifaktorielle ANCOVA in R
Die zweifaktorielle ANCOVA erweitert die einfaktorielle ANCOVA um einen weiteren Faktor. Das ist manchmal etwas verwirrend, da wir auf einmal drei oder mehr Terme in einem Modell haben. Klassischerweise haben wir nun zwei Faktoren \(f_1\) und \(f_2\) in dem Modell. Weiterhin haben wir nur eine Kovariate \(c_1\). Damit sehe das Modell wie folgt aus.
\[ y \sim f_1 + f_2 + c_1 \]
Wir können das Modell dann in R übertragen und ergänzen noch den Interaktionsterm für die Faktoren animal und sex in dem Modell. Das Modell wird klassisch in der Funktion lm() gefittet.
fit_2 <- lm(jump_length ~ animal + sex + weight + animal:sex, data = ancova_tbl)Nach dem Fit können wir das Modell in dem Obkjekt fit_2 in die Funktion anova() pipen. Die Funktion erkennt die Struktur des Modells und gibt uns eine ANCOVA Ausgabe wieder.
fit_2 |> anova() Analysis of Variance Table
Response: jump_length
Df Sum Sq Mean Sq F value Pr(>F)
animal 2 2693.8 1346.9 359.0568 <2e-16 ***
sex 1 3608.1 3608.1 961.8568 <2e-16 ***
weight 1 0.0 0.0 0.0053 0.9422
animal:sex 2 2.2 1.1 0.2981 0.7424
Residuals 593 2224.4 3.8
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In der ANCOVA erkennne wir nun, dass der Faktor animal signifikant ist. Der \(p\)-Wert ist mit \(<0.001\) kleiner das das Signifikanzniveau \(\alpha\) von 5%. Ebenso ist der Faktor sex signifikant. Der \(p\)-Wert ist mit \(<0.001\) kleiner das das Signifikanzniveau \(\alpha\) von 5%. Die Kovariate weight ist nicht mehr signifikant. Der \(p\)-Wert ist mit \(0.94\) größer das das Signifikanzniveau \(\alpha\) von 5%. Wir können also schlussfolgern, dass sich mindestens eine Gruppenvergleich der Level des Faktors animal voneinander unterscheidet. Ebenso wie können wir schlussfolgern, dass sich mindestens eine Gruppenvergleich der Level des Faktors site voneinander unterscheidet. Da wir nur zwei Level in dem Faktor sex haben, wissenwir nun, dass sich die beiden Geschlechter der Flöhe in der Sprungweite unterscheiden. Wir wissen auch, dass mit der Zunahme des Gewichts, sich die Sprunglänge nicht ändert.
In Abbildung 30.2 sehen wir nochmal den Zusammenhang dargestellt. Wenn wir die Daten getrennt für den Faktor sex anschauen, dann sehen wir, dass das Gewicht keinen Einfluss mehr auf die Sprungweite hat.
ggplot(ancova_tbl, aes(weight, jump_length, color = animal)) +
geom_smooth(method = "lm", se = FALSE) +
scale_color_okabeito() +
theme_minimal() +
geom_point() +
labs(color = "Tierart", shape = "Geschlecht") +
facet_wrap(~ sex, scales = "free_x")
Nach einer berechneten ANCOVA können wir zwei Fälle vorliegen haben.
Wenn du in deinem Experiment keine signifikanten Ergebnisse findest, ist das nicht schlimm. Du kannst deine Daten immer noch mit der explorativen Datenanalyse auswerten wie in Kapitel 16 beschrieben.
- Wir habe eine nicht signifkante ANCOVA berechnet. Wir können die Nullhypothese \(H_0\) nicht ablehnen und die Mittelwerte über den Faktor sind vermutlich alle gleich. Wir enden hier mit unserer statistischen Analyse.
- Wir haben eine signifikante ANCOVA berechnet. Wir können die Nullhypothese \(H_0\) ablehnen und mindestens ein Gruppenvergleich über mindestens einen Faktor ist vermutlich unterschiedlich. Wir können dann in Kapitel 34 eine Posthoc Analyse rechnen.
30.6 Linearer Trendtest
“Make use of time, let not advantage slip.” — William Shakespeare
Im Folgenden schauen wir uns dann die Auswertung der Gurkendaten einmal genauer an. Für mehr Informationen zu dem Paket {emmeans} und den entsprechenden Funktionen dann bitte einmal in das Kapitel Kapitel 34 schauen. In der Abbildung 30.3 sehen wir die Daten einmal als Scatterplot. Die durchgezogene Gerade stellt den Verlauf der Mittelwerte über die Versuchsgruppen dar. Die gestrichelte Linie zeigt den Median über die Gruppen. Wir wollen jetzt der Frage nachgehen, ob es einen Unterschied zwischen den Gurkensorten versuchsgruppen über den zeitlichen Verlauf der vierzehn Wochen gibt.
ggplot(gurke_time_len_tbl, aes(time_num, length, color = versuchsgruppe)) +
theme_minimal() +
geom_point() +
stat_summary(fun = "mean", geom = "line") +
stat_summary(fun = "median", geom = "line", linetype = 2) +
scale_color_okabeito()
Unser erstes Modell was wir uns anschauen wollen ist nochmal eine klassische zweifaktorielle ANOVA mit einem Interaktionsterm. Wir wollen raus finden ob die Länge von den Gurken von den Sorten (\(f_1\)), dem zeitlichen Verlauf (\(f_2\)) und der Interaktion zwischen der Sorten und der Zeit abhängt (\(f_1:f_2\)). Wir stellen nun folgendes lineares Modell für die zweifaktorielle ANOVA auf.
\[ length \sim \overbrace{versuchsgruppe}^{f_1} + \underbrace{time\_fct}_{f_2} + \overbrace{versuchsgruppe:time\_fct}^{f_1:f_2} \]
Dieses Modell können wir dann auch in R einmal über die Funktion lm() abbilden.
time_fct_fit <- lm(length ~ versuchsgruppe + time_fct + versuchsgruppe:time_fct,
gurke_time_len_tbl)Nun wollen wir auch überprüfen, ob es eine Interaktion zwischen den Versuchsgruppen und dem zeitlichen Verlauf gibt. Das ganze schauen wir uns neben einer ANOVA auch einmal graphisch mit der Funktion emmip() an. Wenn wir keine signifikante Interaktion erwaten würden, dann müssten die drei Versuchgruppe über den zeitlichen Verlauf gleichmäßig ansteigen. Wir sehen in der Abbildung 30.4, dass dies nicht der Fall ist. Wir nehmen daher eine Interaktion zwischen den Versuchsgruppen und dem zeitlichen Verlauf an.
emmip(time_fct_fit, versuchsgruppe ~ time_fct, CIs = TRUE) +
theme_minimal() +
scale_color_okabeito()
Nochmal kurz mit der ANOVA überprüft und wir sehen eine signifikante Interaktion.
time_fct_fit |> anova() |> model_parameters()Parameter | Sum_Squares | df | Mean_Square | F | p
---------------------------------------------------------------------------
versuchsgruppe | 3296.84 | 2 | 1648.42 | 169.67 | < .001
time_fct | 4184.15 | 13 | 321.86 | 33.13 | < .001
versuchsgruppe:time_fct | 1628.12 | 26 | 62.62 | 6.45 | < .001
Residuals | 1981.94 | 204 | 9.72 | |
Anova Table (Type 1 tests)Wir haben eben einmal die Zeit als Faktor mit ins Modell genommen, da wir so für jeden Zeitpunkt einen Mittelwert schätzen können. Wenn wir eine einfaktorielle ANCOVA rechnen, dann geht die Zeit als numerische Kovariate (\(c_1\)) linear in das Modell ein. Wir haben also von jedem Zeitpunkt zum nächsten den gleichen Anstieg. Wir modellieren ja auch einen linearen Zusammenhang. Hier einmal das Modell für die ANCOVA.
\[ length \sim \overbrace{versuchsgruppe}^{f_1} + \underbrace{time\_num}_{c_1} + \overbrace{versuchsgruppe:time\_num}^{f_1:c_1} \]
Wir fitten wieder das Modell in R mit der Funktion lm().
time_num_fit <- lm(length ~ versuchsgruppe + time_num + versuchsgruppe:time_num, gurke_time_len_tbl)Auch hier einmal die Überprüfung auf eine Interaktion mit der ANCOVA. Wir sehen, dass wir eine signifikante Interaktion zwischen den Versuchsgruppen und dem zeitlichen verlauf vorliegen haben.
time_num_fit |> anova() |> model_parameters()Parameter | Sum_Squares | df | Mean_Square | F | p
---------------------------------------------------------------------------
versuchsgruppe | 3296.84 | 2 | 1648.42 | 155.05 | < .001
time_num | 3816.42 | 1 | 3816.42 | 358.98 | < .001
versuchsgruppe:time_num | 1426.30 | 2 | 713.15 | 67.08 | < .001
Residuals | 2551.50 | 240 | 10.63 | |
Anova Table (Type 1 tests)In der Abbildung 30.5 sehen wir dann nochmal das Modell im Interaktionsplot. In beiden Abbildungen sehe wir den linearen Zusammenhang. Die Interaktion drückt sich in der unterschiedlichen Steigung der Versuchsgruppen über den zeitlichen Verlauf aus. Wir können durch die Option at = entscheiden für welche Zeitpunkte wir uns die 95% Konfidenzintervalle anzeigen lassen wollen, wie in der linken Abbildung exemplarisch für die Zeitpunkte 1 Woche, 7 Wochen und 14 Wochen. Oder aber wir lassen uns alle durch die Option cov.reduce = FALSE anzeigen, wie in der rechten Abbildung gezeigt.
emmip(time_num_fit, versuchsgruppe ~ time_num, CIs = TRUE,
at = list(time_num = c(1, 7, 14))) +
theme_minimal() +
scale_x_continuous(breaks = 1:14) +
scale_color_okabeito()
emmip(time_num_fit, versuchsgruppe ~ time_num, CIs = TRUE,
cov.reduce = FALSE) +
theme_minimal() +
scale_x_continuous(breaks = 1:14) +
scale_color_okabeito()
Jetzt kommt eigentlich der spannende Teil, wir wollen jetzt über den zeitlichen Verlauf einen statistischen Test rechnen, ob wir einen Trend in den verschiedenen Versuchsgruppen haben. Wir rechnen also einen Trendtest über die Kovariate \(c_1\) für die drei Gruppen getrennt. Das können wir mit der Funktion emtrends() durchführen. Hier musst du angeben welche deine Kovariate var ist. In unserem Fall ist die Kovariate natürlich time_num.
emtrends(time_num_fit, ~ versuchsgruppe, var = "time_num", infer = TRUE) versuchsgruppe time_num.trend SE df lower.CL upper.CL t.ratio p.value
Proloog 1.857 0.0923 240 1.675 2.039 20.116 <.0001
Quarto 0.465 0.0883 240 0.291 0.639 5.267 <.0001
Katrina 0.699 0.0897 240 0.522 0.876 7.794 <.0001
Confidence level used: 0.95 Die Spalte time_num.trend zeigt uns jetzt den linearen Anstieg über die Zeit für die jeweilige Versuchsgruppe. Das heißt jede Woche wächst die Sorte Proloog um \(1.857cm\) an Länge. In der gleichen Art können wir auch die anderen Werte in der Spalte interpretieren. Jetzt stellt sich natürlich die Frage, ob dieses Längenwachstum untereinander unterschiedlich ist. Dafür könne wir dann leicht den Code abändern und setzen das Wort pairwise vor die Tilde. Dann testet die Funktion auch alle paarweisen Vergleiche. Wir nutzen jetzt noch die Funktion model_parameters() um eine schönere Ausgabe zu erhalten.
emtrends(time_num_fit, pairwise ~ versuchsgruppe, var = "time_num", infer = TRUE) |>
model_parameters()# emtrends
Parameter | Coefficient | SE | 95% CI | t(240) | p
---------------------------------------------------------------
Proloog | 1.86 | 0.09 | [1.68, 2.04] | 20.12 | < .001
Quarto | 0.46 | 0.09 | [0.29, 0.64] | 5.27 | < .001
Katrina | 0.70 | 0.09 | [0.52, 0.88] | 7.79 | < .001
# Contrasts
Parameter | Coefficient | SE | 95% CI | t(240) | p
------------------------------------------------------------------------
Proloog - Quarto | 1.39 | 0.13 | [ 1.14, 1.64] | 10.90 | < .001
Proloog - Katrina | 1.16 | 0.13 | [ 0.90, 1.41] | 9.00 | < .001
Quarto - Katrina | -0.23 | 0.13 | [-0.48, 0.01] | -1.86 | 0.153 Wir sehen wieder erst die Trends, die kennen wir schon. Dann sehen wir die Kontraste oder auch paarweisen Vergleiche. Das kannst du schnelle nachrechnen, der Unterschied von Proloog zu Quarto ist \(1.48 - 0.45 = 1.03\). Damit können wir dann auch über den gesamten zeitlichen Verlauf testen, ob ein Unterschied zwischen den Sorten vorliegt.
Wir sind jetzt schon sehr weit gekommen, aber wir könnten auch einen nicht-linearen Zusammenhang zwsichen der Zeit und dem Längenwachstum von Gurken annehmen. Das würde auch biologisch etwas mehr Sinn ergeben. Deshalb modellieren wir den Einfluss der Zeit time_num durch einen Exponenten hoch drei. Daher schreiben wir mathematisch \((time\_num)^3\).
\[ length \sim \overbrace{versuchsgruppe}^{f_1} + \underbrace{(time\_num)^3}_{c_1} + \overbrace{versuchsgruppe:(time\_num)^3}^{f_1:c_1} \]
Den Exponenten schreiben wir dann entsprechend in R mit poly(time_num, 3) in die Formel. Die Funktion poly() nutzen wir innerhalb von Formelaufrufen um einen Exponenten einzufügen.
time_num_poly_fit <- lm(length ~ versuchsgruppe * poly(time_num, 3), gurke_time_len_tbl)In der Abbildung 30.6 sehen wir dann nochmal die Anpassung des Modells und wir sehen, dass unser Modell besser zu den Daten passt. Das sieht schon sehr viel sauberer aus, als der brutale lineare Zusammenhang, den wir vorher hatten. Du musst hier etwas mit den Exponenten spielen und ausprobieren, welcher da am besten passt.
emmip(time_num_poly_fit, versuchsgruppe ~ time_num, CIs = TRUE,
cov.reduce = FALSE) +
theme_minimal() +
scale_color_okabeito()
Jetzt wollen wir die Analyse nochmal auf die Spitze treiben. Wir schauen uns jetzt an den Zeitpunkten 1 Woche, 7 Wochen und 14 Wochen den Unterschied zwischen den Sorten einmal an. Das machen wir indem wir zu der Funktion emmeans() die Optione at = list(time_num = c(1, 7, 14)) ergänzen. Am Ende lassen wir uns noch das compact letter display wiedergeben.
emmeans(time_num_poly_fit, pairwise ~ versuchsgruppe | time_num,
at = list(time_num = c(1, 7, 14))) |>
cld(Letters = letters) time_num = 1:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 2.20 1.008 234 0.209 4.18 a
Katrina 3.07 1.009 234 1.077 5.05 a
Proloog 5.10 1.011 234 3.106 7.09 a
time_num = 7:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 3.90 0.492 234 2.932 4.87 a
Katrina 5.32 0.492 234 4.348 6.29 a
Proloog 10.40 0.496 234 9.421 11.38 b
time_num = 14:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 8.82 1.008 234 6.832 10.80 a
Katrina 13.15 1.073 234 11.033 15.26 b
Proloog 30.86 1.097 234 28.698 33.02 c
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 3 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Was erkennen wir? Wir sehen, dass sich in der ersten Woche die Sorten noch nicht voneinander unterscheiden. Erst in der Woche sieben sehen wir einen Unterschied von Proloog zu dem Rest der Sorten. In der letzten Woche unterscheiden sich dann alle Sorten voneinander. Wie stark, kannst du aus der Spalte emmean entnehmen, dort steht der Mittelwert für die jeweilige Sorte zu dem jeweiligen Zeitpunkt.
Abschließend wollen wir nochmal schauen, wie sich der Trend in den verschiedenen Modellierungen der Exponenten zeigen würde. Wir müssen dafür bei der Funktion emtrends() angeben bis zu welchen maximalen Exponenten, bei uns hoch Drei, die Funktion rechnen soll. Deshalb setzen wir max.degree = 3. Dann noch pairwise vor die Tilde gesetzt, damit wir dann auch die paarweisen Vergleiche für die Sorten und Verläufe angezeigt kriegen. Achtung, jetzt kommt eine lange Ausgabe.
emtrends(time_num_poly_fit, pairwise ~ versuchsgruppe, var = "time_num", infer = TRUE,
max.degree = 3) $emtrends
degree = linear:
versuchsgruppe time_num.trend SE df lower.CL upper.CL t.ratio p.value
Proloog 1.64460 0.21158 234 1.22775 2.0615 7.773 <.0001
Quarto 0.35912 0.20333 234 -0.04148 0.7597 1.766 0.0787
Katrina 0.54881 0.20467 234 0.14559 0.9520 2.682 0.0079
degree = quadratic:
versuchsgruppe time_num.trend SE df lower.CL upper.CL t.ratio p.value
Proloog 0.15795 0.02308 234 0.11249 0.2034 6.845 <.0001
Quarto 0.03267 0.02245 234 -0.01156 0.0769 1.455 0.1469
Katrina 0.05790 0.02285 234 0.01289 0.1029 2.534 0.0119
degree = cubic:
versuchsgruppe time_num.trend SE df lower.CL upper.CL t.ratio p.value
Proloog 0.00714 0.00671 234 -0.00608 0.0204 1.064 0.2883
Quarto 0.00338 0.00644 234 -0.00931 0.0161 0.525 0.5999
Katrina 0.00505 0.00657 234 -0.00789 0.0180 0.770 0.4424
Confidence level used: 0.95
$contrasts
degree = linear:
contrast estimate SE df lower.CL upper.CL t.ratio p.value
Proloog - Quarto 1.28548 0.29345 234 0.5933 1.9776 4.381 0.0001
Proloog - Katrina 1.09579 0.29437 234 0.4014 1.7901 3.722 0.0007
Quarto - Katrina -0.18969 0.28850 234 -0.8702 0.4908 -0.658 0.7883
degree = quadratic:
contrast estimate SE df lower.CL upper.CL t.ratio p.value
Proloog - Quarto 0.12527 0.03220 234 0.0493 0.2012 3.891 0.0004
Proloog - Katrina 0.10005 0.03247 234 0.0235 0.1766 3.081 0.0065
Quarto - Katrina -0.02522 0.03203 234 -0.1008 0.0503 -0.787 0.7111
degree = cubic:
contrast estimate SE df lower.CL upper.CL t.ratio p.value
Proloog - Quarto 0.00375 0.00930 234 -0.0182 0.0257 0.404 0.9141
Proloog - Katrina 0.00208 0.00939 234 -0.0201 0.0242 0.222 0.9732
Quarto - Katrina -0.00167 0.00920 234 -0.0234 0.0200 -0.182 0.9820
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 3 estimates
P value adjustment: tukey method for comparing a family of 3 estimates Was sehen wir in diesem Fall? In der linearen Modellierung der Zeit haben wir sehr viele signifikante Ergebnisse. Leider entspricht der lineare verlauf über die Zeit nicht so den beobachteten Daten. Bei der kubischen Modellierung, also hoch Drei, haben wir dann in der Abbildung eine bessere Modellierung. Die Effekte reichen dann aber nicht aus um über den gesamten zeitlichen Verlauf, Betonung liegt auf dem gesamten zeitlichen Verlauf, einen Unterschied zeigen zu können.
Daher müsste man hier einmal überlegen, ob man nicht die frühen Wochen aus den Daten entfernt. Hier sind die Gurlen noch sehr ähnlich, so dass wir hier eigentlich auch keinen Unterschied erwarten. Sonst könntest du nochmal mit einer \(\log\)-Transformation der Länge spielen, dann verlierst du zwar die direkte biologische Interpretierbarkeit der Effektschätzer, aber dafür könnte der Verlauf über die Zeit besser aufgesplittet werden. Das Problem sind ja hier sehr kleine Werte zu Anfang und sehr große Werte zum Ende.