pacman::p_load(tidyverse, magrittr, conflicted, broom, quantreg,
see, performance, emmeans, multcomp, janitor,
parameters, effectsize, MASS, modelsummary,
robustbase, multcompView, conflicted)
conflicts_prefer(dplyr::select)
conflicts_prefer(dplyr::filter)
#cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
# "#F0E442", "#0072B2", "#D55E00", "#CC79A7")56 Robuste und Quantilesregression
Letzte Änderung am 02. April 2024 um 09:52:43
“Methods for robust statistics, a state of the art in the early 2000s, notably for robust regression and robust multivariate analysis.” — robust: Port of the S+ “Robust Library”
In diesem Kapitel geht es um zwei Arten der Regression, die immer wieder genannt werden, aber dennoch eine Art Nischendasein fristen. Zum einen möchte ich hier die robuste Regression (eng. robust regression) und zum anderen die Quantilsregression (eng. quantile regression) vorstellen. Die robuste Regression ist faktisch statistisch tot. Das heißt, die Implementierungen werden kaum weiterentwickelt und auch findet methodische Forschung nur in der theoretischen Nische statt. Zwar wird die robuste Regression in ihrer ursprünglichen Form als Regression angewendet, aber das reicht meistens nicht. Selten wollen wir nur durch ein paar Punkte eine Gerade ziehen und uns über das gute Modell erfreuen. Wir haben mit dem Modell meist mehr vor. Wir wollen eine ANOVA rechnen und dann auch einen wie auch immer gearteten Mittelwertsvergleich. Wenn dies zwar theoretisch möglich ist, praktisch aber nicht implementiert, dann wollen und können wir die Methoden nur eingeschränkt verwenden. Bei der Quantilsregression sieht es etwas anders aus, hier können wir dann schon den ein oder anderen Mittelwertsvergleich rechnen. Was bei der Quantilsregression eher problematisch ist, ist das das Modell nicht immer konvergiert oder aber nicht algorithmisch eine Lösung für einen spezifischen Datensatz findet. In diesem Datensatz dann den einen Grund zu finden, ist dann meist so aufwendig, dass wir es auch gleich mit der Quantilsregression lassen können.
Der Charme der Quantilesregression ist ja am Ende, dass wir auch den Median als Quantile auswählen können. So haben wir dann die Möglichkeit eine Regression auf den Medianen zu rechnen. Wir vergleichen damit dann auch die Mediane und sind so nicht mehr auf die Normalverteilung unseres Outcomes wie in der gewöhnlichen Regression angewiesen. Also eigentlich eine tolle Sache, wenn wir nur an den Mittelwertsvergleichen interessiert sind. Eine klassische ANOVA geht leider nicht auf einer Medianregression. Eine klassische ANOVA wäre aber mit nicht parametrischen Methoden sowieso nicht möglich gewesen. Wir verlieren also nicht so viel, gewinnen aber etwas, wenn das Modell konvergiert und ein Ergebnis liefert.
In diesem Kapitel brechen wir etwas die bisherige Struktur der Regressionskapitel auf. Wir schauen uns hier zuerst die beiden Modelle an und entscheiden dann, ob wir die Modelle für die ANOVA oder den Gruppenvergleich überhaupt nutzen können. Am Ende vergleichen wir dann einmal alle Modell mit dem fantastische Paket {modelsummary} mit der gleichnamigen Funktion. Hier hilft dann wie immer die tolle Hilfsseite von modelsummary zu besuchen.
56.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.
56.2 Daten
Für unser erstes Beispiel nutzen wir die Daten aus einem Wachstumsversuch mit Basilikum mit vier Bodenbehandlungen und vier Blöcken, die Blöcke sind eigentlich Gewächshaustische. Wir wollen also einen klassischen Gruppenvergleich mit Berücksichtigung der Blockstruktur rechnen.
basi_tbl <- read_excel("data/keimversuch_basilikum_block.xlsx") |>
clean_names() |>
mutate(versuchsgruppe = as_factor(versuchsgruppe)) |>
select(versuchsgruppe, block_1:block_4)In Tabelle 56.1 sehen wir einmal die Daten im Wide-Format. Wir haben also das Frischgewicht der Basilikumpflanzen gemessen und wollen wissen, ob die verschiedenen Bodenarten einen Einfluss auf das Wachstum haben.
| versuchsgruppe | block_1 | block_2 | block_3 | block_4 |
|---|---|---|---|---|
| Erde | 16 | 21 | 23 | 23 |
| Erde | 17 | 19 | 18 | 24 |
| Erde | 16 | 22 | 23 | 24 |
| Erde | 9 | 17 | 18 | 21 |
| Erde | 17 | 21 | 22 | 24 |
| Erde+Fließ | 18 | 22 | 21 | 21 |
| … | … | … | … | … |
| Erde+Perlite | 20 | 25 | 25 | 25 |
| Perlite+Fließ | 22 | 25 | 25 | 24 |
| Perlite+Fließ | 25 | 25 | 26 | 26 |
| Perlite+Fließ | 15 | 19 | 19 | 19 |
| Perlite+Fließ | 17 | 22 | 22 | 22 |
| Perlite+Fließ | 22 | 22 | 22 | 22 |
Da wir die Daten im Wide-Format vorliegen haben, müssen wir die Daten nochmal in Long-Format umwandeln. Wie immer nutzen wir dafür die Funktion pivot_longer().
basi_block_tbl <- basi_tbl |>
pivot_longer(cols = block_1:block_4,
values_to = "weight",
names_to = "block") |>
mutate(block = as_factor(block))In der Abbildung 56.1 siehst du einmal die Daten als Dotplots mit Mittelwert und Standardabweichung. Wir machen hier mal einen etwas komplizierteren Plot, aber immer nur Barplot ist ja auch langweilig.
ggplot(basi_block_tbl, aes(versuchsgruppe, weight, color = block)) +
theme_minimal() +
scale_color_okabeito() +
geom_point(position = position_dodge(0.5), shape = 4, size = 2.5) +
stat_summary(fun.data="mean_sdl", fun.args = list(mult = 1),
geom="pointrange", position = position_dodge(0.5)) +
stat_summary(fun = "mean", fun.min = "min", fun.max = "max", geom = "line",
position = position_dodge(0.5)) 
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. Wir ignorieren hier einmal die abhängige Datenstruktur über die Zeitpunkte.
Mit einer abhängigen Datenstruktur müssten wir eigentlich ein lineares gemischtes Modell rechnen. Aber wir nutzen hier die Daten einmal anders.
Im Weiteren haben wir 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. Am Ende schmeißen wir noch Spalten raus, die wir nicht weiter brauchen.
gurke_raw_tbl <- read_excel("data/wachstum_gurke.xlsx") |>
clean_names() |>
filter(str_detect(versuchsgruppe, "L$")) |>
select(-pfl, -erntegewicht) |>
mutate(versuchsgruppe = factor(versuchsgruppe,
labels = c("Katrina", "Proloog", "Quarto"))) 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 | 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 | 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 | 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 | 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 | 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 | 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 | 3.7 | 3.9 | 3.9 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4.7 | 5.2 | 5.4 | 5.3 | 5.4 | 5 |
| Katrina | 3 | 3.2 | 3.3 | 3.5 | 3.6 | 4.4 | 5.2 | 5.8 | 6 | 6.2 | 6.1 | 6.2 | 6.8 | 7.9 | 9.4 | 10.4 | 13.2 |
| Katrina | 3.3 | 3.3 | 3.4 | 3.6 | 3.5 | 3.5 | 3.4 | 3.7 | 3.5 | 3.6 | 3.6 | 3.6 | 2.8 | 0 | NA | NA | NA |
| Katrina | 3.2 | 3.4 | 3.7 | 4 | 4.5 | 5.9 | 6.9 | 8.4 | 9.4 | 11 | 12.2 | 13.5 | 17.9 | 18 | NA | NA | NA |
| Katrina | 3.3 | 3.4 | 3.9 | 4.6 | 5.2 | 6.5 | 7.9 | 9.5 | 10.5 | 11.7 | 13 | 13.4 | 17.9 | 18 | NA | NA | NA |
| Katrina | 3.1 | 3.6 | 3.6 | 3.8 | 4.2 | 4.9 | 5.7 | 6.8 | 7.7 | 9.1 | 10.7 | 12.3 | 17.7 | 18.6 | NA | NA | NA |
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 Gruken lassen wir jetzt mal drin, da wir ja eine robuste Regression noch rechnen wollen. 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 |>
pivot_longer(cols = t1:t17,
values_to = "length",
names_to = "time") |>
mutate(time_fct = as_factor(time),
time_num = as.numeric(time_fct)) |>
filter(time_num <= 14)In der Abbildung 56.2 sehen wir dann nochmal den Scatterplot für das Gurkenwachstum. Die gestrichtelten Linien stellen den Median und die durchgezogene Line den Mittelwert der Gruppen dar.
ggplot(gurke_time_len_tbl, aes(time_num, length, color = versuchsgruppe)) +
theme_minimal() +
geom_point2(position = position_dodge(0.5)) +
stat_summary(fun = "mean", geom = "line") +
stat_summary(fun = "median", geom = "line", linetype = 2) +
scale_x_continuous(breaks = 1:14) +
scale_color_okabeito()
56.3 Gewöhnliche lineare Regression
Damit wir einen Vergleich zu der robusten Regression und der Quantilsregression haben wollen wir hier zu Anfang nochmal schnell die gewöhnliche Regression (eng. ordinary linear regression), die du schon durch die Funktion lm() kennst. Wir fitten also einmal das Modell mit Interaktionsterm für den Basilikumdatensatz.
basi_lm_fit <- lm(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)56.3.1 ANOVA
Dann schauen wir uns nochmal den Interaktionsplot in Abbildung 56.3 mit der Funktion emmip() aus dem R Paket {emmeans} an. Wir sehen, das wir eventuell eine leichte Interaktion vorliegen haben könnten, da sich einige der Geraden überschneiden.
emmip(basi_lm_fit, versuchsgruppe ~ block, CIs = TRUE,
cov.reduce = FALSE) +
theme_minimal() +
scale_color_okabeito()
Schauen wir einmal in eine zweifaktorielle ANOVA, ob sich unser leichter Verdacht validieren lässt.
basi_lm_fit |>
anova() |>
model_parameters()Parameter | Sum_Squares | df | Mean_Square | F | p
----------------------------------------------------------------------
versuchsgruppe | 158.84 | 3 | 52.95 | 6.53 | < .001
block | 355.34 | 3 | 118.45 | 14.60 | < .001
versuchsgruppe:block | 83.51 | 9 | 9.28 | 1.14 | 0.346
Residuals | 519.20 | 64 | 8.11 | |
Anova Table (Type 1 tests)Wir sehen, dass wir mindestens einen paarweisen Unterschied zwischen den Versuchsgruppen und den Blöcken erwarten. Die Interaktion ist nicht signifikant. Betrachten wir noch kurz das \(\eta^2\) um zu sehen, wie viel Varianz jeweils die Versuchsgruppe und Blöcke erklären.
basi_lm_fit |>
eta_squared()# Effect Size for ANOVA (Type I)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------------
versuchsgruppe | 0.23 | [0.08, 1.00]
block | 0.41 | [0.24, 1.00]
versuchsgruppe:block | 0.14 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Leider erklären hier die Blöcke sehr viel der Varianz, das ist nicht so schön, aber in diesem Kapitel lassen wir es dabei.
56.3.2 Gruppenvergleich
Wir haben keine Interaktion vorliegen, alos mitteln wir einmal über alle Blöcke. Wenn du den Code | block zu dem Aufruf der emmeans() Funktion hinzufügst, dann hast du die Analyse für die Blöcke getrennt durchgeführt.
basi_lm_fit |>
emmeans(specs = ~ versuchsgruppe) |>
cld(Letters = letters, adjust = "none") versuchsgruppe emmean SE df lower.CL upper.CL .group
Erde+Fließ 19.4 0.637 64 18.1 20.7 a
Erde 19.8 0.637 64 18.5 21.0 a
Perlite+Fließ 22.1 0.637 64 20.8 23.3 b
Erde+Perlite 22.6 0.637 64 21.4 23.9 b
Results are averaged over the levels of: block
Confidence level used: 0.95
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. Wenn wir uns die unadjustierten \(p\)-Werte anschauen, dann sehen wir einen leichten Effekt zwischen den einzelnen Versuchsgruppen. Warum sage ich leichten Effekt? Die Mittelwerte der Gruppen in der spalte emmean unterscheiden sich kaum. Aber soviel zu der gewöhnlichen Regression, das war ja hier nur der Vergleich und die Erinnerung.
56.4 Robuste Regression
Nur Normalverteilung? Geht auch mehr?
Nein es geht auch mit allen anderen Modellen. So ist auch das glm() für die robuste Regression implementiert. Wir schauen uns aber nur die Grundlagen für die klassische robuste Regression unter der Annahme eines nromalverteilten Outcomes an.
Bei einer robusten Regression werden jeder Beobachtung auf der Grundlage ihres Residuums unterschiedliche Gewichte (eng. robustness weights) von 0 bis 1 zugewiesen. Wir kennen ja ein Residuum als die Differenz zwischen den beobachteten und den vorhergesagten Werten der Daten. Die vorhergesagten Daten sind ja die Punkte auf der Geraden. Je kleiner also das Residuum ist, desto größer ist die Gewichtung und desto näher liegt eine Beobachtung an der geschätzten Geraden. Wir bestrafen also Punkte, die weit weg von unserer potenziellen Gerade liegen.
Ein englisches Tutorium gibt es dann nochmal ausführlicher unter R demo | Robust Regression (don’t depend on influential data!)
Ebenso liefert auch das Tutorium Robust regression | R data analysis example einen Überblick.
Schauen wir uns das Problem an einem kleinen Spieldatensatz einmal an, den es so in vielen Geschmacksrichtungen gibt. Wir haben acht Beobachtungen mit jeweils einem \(x\) und einem \(y\) Wert in unserem Datensatz reg_tbl. Wie folgt einmal dargestellt.
reg_tbl <- tibble(x = c(0, 1, 2, 3, 4, 5, 6, 7),
y = c(1, 1.4, 2.5, 2.7, 4.3, 5.2, 0, 6.7))Wir rechnen jetzt einmal eine gewöhnliche lineare Regression mit der Funktion lm() und eine robuste Regression mit dem Paket {MASS} und der Funtkion rlm(). Es gibt noch das R Paket {robustbase} und der Funktion rq() aber wir können die Ergebnisse dann nicht in {emmeans} weiter nutzen, was uns dann etwas den Weg mit versperrt. Wir können aber für die multiplen Mittelwertsvergleich mit {robustbase} das R Paket {multcomp} nutzen. Aber fangen wir hier erstmal mit dem Beispiel an und vergleichen die beiden Implementierungen der robusten Regression zu der gewöhnlichen Regression
reg_lm_fit <- lm(y ~ x, reg_tbl)
reg_rlm_fit <- rlm(y ~ x, reg_tbl)
reg_lmrob_fit <- lmrob(y ~ x, reg_tbl)In der Abbildung 56.4 sehen wir einmal den Unterschied zwischen den beiden robusten Regressionen und der gewöhnlichen Regression.
ggplot(reg_tbl, aes(x, y, label = x)) +
theme_minimal() +
geom_label() +
geom_line(aes(y = predict(reg_lm_fit), colour = "lm")) +
geom_line(aes(y = predict(reg_rlm_fit), colour = "rlm")) +
geom_line(aes(y = predict(reg_lmrob_fit), colour = "lmrob")) +
scale_colour_okabeito() +
labs(colour = "Method")
In unserem Beispiel ist die sechste Beobachtung der Ausreißer. Die lm-Regression wird durch den Wert 0 der sechsten Beobachtung nach unten gezogen. Die sechste Beobachtung erhält ein viel zu starkes Gewicht in der gewöhnlichen Regression im Bezug auf den Verlauf der Geraden. Die robusten Regressionen geben der sechsten Beobachtung hingegen ein viel geringeres Gewicht, so dass die Beobachtung keinen Einfluss auf den Verlauf der Geraden hat. In einer robusten Regression erhält eine Beobachtung mit einem großen Residuum sehr wenig Gewicht, in diesem konstruierten Beispiel fast 0. Daher können wir sagen, das unser Modell robust gegen Ausreißer ist.
Okay, soweit lernt es jeder und die Idee der robusten Regression ist auch ziemlich einleuchtend. Aber in der Wissenschaft zeichnen wir keine Gerade durch Punkte und freuen uns drüber, dass die eine gerade besser passt als die andere Gerade. Wir nutzen dann ja das Modell weiter in eine ANOVA oder aber eben in einem paarweisen Gruppenvergleich. Und hier fangen dann die Probleme der Implementierung an. Leider ist es so, dass nicht alle Funktionalitäten implementiert sind. Das heißt, die Funktion anova() oder emmeans() kann mit der Ausgabe der robusten Regressionen mit MASS::rlm() und robustbase::lmrob() nichts anfangen. Das ist problematisch und nicht anwenderfreundlich. Und wir sind hier nunmal Anwender…
Warum so unbeliebt? Das ist eine gute Frage. Zu der robusten Regression gibt es eine lange Geschichte mit vielen Wirrungen und Irrungen. Am Ende kannst du in Wikipedia noch etwas zur History and unpopularity of robust regression nachlesen. Nach Stromberg (2004) hat es auch etwas mit der schlechten und mangelhaften Implementierung in gängiger Statistiksoftware zu tun. Das merke ich hier auch gerade, als ich versuche die Funktionen zusammen zubauen und miteinander zu vernetzen. Sicherlich gibt es auch die biologische Diskussion. Gibt es in biologischen Daten eigentlich überhaupt Ausreißer? Oder müssten wir nicht Ausreißer besonders betrachten, weil die Ausreißer dann doch mehr über die Daten verraten? Macht es also Sinn die Ausreißer einfach platt zu bügeln und nicht mehr weiter zu beachten oder ist nicht die Beobachtung Nr.6 oben von Interesse? Denn vielleicht ist was mit dem Experiment schief gelaufen und diese Beobachtung ist gerade ein Anzeichen dafür. Dann müssten auch die anderen Messwerte kritisch hinterfragt werden.
Wie immer es auch sein. Wir fitten jetzt mal die beiden robusten Regression in der Standardausstattung mit der Funktion MASS::rlm() aus dem Paket {MASS} wie folgt.
basi_rlm_fit <- rlm(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)Dann natürlich noch die Funktion robustbase::lmrob() aus dem Paket {robustbase}. Hinter dem Paket {robusbase} steht ein ganzen Buch von Maronna u. a. (2019), aber leider keine gute Internetseite mit Beispielen und Anwendungen in R, wie es eigentlich heutzutage üblich ist. Aber dennoch, hier einmal das Modell.
basi_lmrob_fit <- lmrob(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)56.4.1 ANOVA
Dann wollen wir einmal eine zweifaktorielle ANOVA auf den Modell der robusten Regression mit der Funktion MASS::rlm() rechnen. Leider kann ich schon gleich sagen, dass wir hier ein Problem kriegen. Es ist nicht ganz klar unter Theoretikern welche Anzahl an Beobachtungen für die Residuen in der robusten Regression genutzt werden soll. Welche Beobachtungen werden den im Modell berücksichtigt und welche nicht? Daher ist dann der Wert für die Freiheitsgrade \(df\) leer und wir erhalten keine \(F\)-Statistik und dann auch keine \(p\)-Werte. Für mich super unbefriedigend, so kann man dann auch keine Funktionen richtig nutzen. Die Funktion robustbase::lmrob() bietet die ANOVA nicht an.
basi_rlm_aov <- basi_rlm_fit |>
anova() |>
tidy()
basi_rlm_aov# A tibble: 4 x 6
term df sumsq meansq statistic p.value
<chr> <int> <dbl> <dbl> <dbl> <dbl>
1 versuchsgruppe 3 114. 37.9 NA NA
2 block 3 318. 106. NA NA
3 versuchsgruppe:block 9 63.7 7.07 NA NA
4 Residuals NA 527. NA NA NAWas als Freiheitsgrad für die Residuen nehmen? In unserem Datensatz haben wir 80 Basilikumpflanzen und somit auch Beobachtungen vorliegen. Wir setzen jetzt einfach brutal die Freiheitsgrade df für die Residuen auf 80 und können dann die Mean Squares für die Residuen berechnen. Also einmal flott die meansq für die Residuen durch die Anzahl an Beobachtungen geteilt.
ms_resid <- 527 / 80Mit den Mean Square der Residuen können wir dann auch die \(F\)-Werte berechnen. Den rechnerischen Weg kennen wir ja schon aus der zweifaktoriellen ANOVA.
f_calc <- basi_rlm_aov$meansq / ms_residWie kriegen wir jetzt aus den ganzen berechneten \(F\)-Werten die entsprechenden \(p\)-Werte? Die \(p\)-Werte sind die Fläche rechts von den berechneten Werten. Wir können die Funktion pf() nutzen um diese Flächen zu berechnen. Wir machen das einmal beispielhaft für einen \(F_{calc} = 3.41\) für 4 Gruppen mit \(df_1 = 4 - 1 = 3\) und 12 Beobachtungen \(df_2 = 12\). Die Berechnung von \(df_2\) ist statistisch falsch aber für unseren MacGyver Hack hinreichend.
pf(3.41, 3, 12, lower.tail = FALSE) |> round(3)[1] 0.053Und siehe da, wir erhalten einen \(p\)-Wert von \(0.053\). Nun können wir das nicht nur für einzelnen \(F\)-Werte machen sondern auch für unseren ganzen Vektor f_calc. Als \(df_2\) setzen wir die Anzahl an Basilikumpflanzen.
p_vec <- pf(f_calc, c(3, 3, 9), c(80, 80, 80), lower.tail = FALSE) |>
scales::pvalue()Die berechneten \(F\)-Werte und \(p\)-Werte können wir jetzt in die ANOVA Tabelle ergänzen.
basi_rlm_aov |>
mutate(statistic = f_calc,
p.value = p_vec)# A tibble: 4 x 6
term df sumsq meansq statistic p.value
<chr> <int> <dbl> <dbl> <dbl> <chr>
1 versuchsgruppe 3 114. 37.9 5.75 0.001
2 block 3 318. 106. 16.1 <0.001
3 versuchsgruppe:block 9 63.7 7.07 1.07 0.391
4 Residuals NA 527. NA NA <NA> Das war jetzt ein ganz schöner Angang und den \(p\)-Werten ist nur approximativ zu trauen, aber wenn wir weit vom Signifikanzniveau \(\alpha\) gleich 5% sind, dann macht eine numerische Ungenauigkeit auch nicht viel aus. Achtung deshalb bei \(p\)-Werten nahe des Signifikanzniveau \(\alpha\), hier ist dann die Aussage mit Vorsicht zu treffen oder eher nicht möglich.
Glücklicherweise können wir aus den Sum of squares dann unsere \(\eta^2\)-Werte für die erklärte Varianz unserer ANOVA berechnen. Das ist dann doch schon mal was. Wir sehen, dass die Blöcke am meisten der Varianz erklären und die Versuchsgruppen sehr wenig. Das ist für das Modell dann nicht so schön, deckt sich aber mit den Ergebnissen gewöhnlichen Regression. Wir erhalten eine Konfidenzintervalle, woran du schon erkennen kannst, dass in der Ausgabe der ANOVA der robusten Regression was fehlt.
basi_rlm_fit |>
eta_squared()# Effect Size for ANOVA (Type I)
Parameter | Eta2 (partial) | 95% CI
----------------------------------------------------
versuchsgruppe | 0.23 | [0.08, 1.00]
block | 0.41 | [0.24, 1.00]
versuchsgruppe:block | 0.14 | [0.00, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Nochmal als Erinnerung, mit der Funktion robustbase::lmrob() ist die ANOVA nicht möglich. Die Funktionalität ist nicht implementiert.
56.4.2 Gruppenvergleich
Dafür geht der Gruppenvergleich für die robuste Regression mit MASS::rlm() wie von alleine. Wir haben keine Interaktion vorliegen, daher müssen wir auch nicht den Block berücksichtigen. Das wir in Block 1 ein Problem haben, steht dann aber auf einem anderen Blatt. Da aber alle Basilikumpflanzen immer im Block 1 kleiner sind, passt es dann wieder im Sinne keiner Interaktion. Mit Interaktion hätten wir sonst versuchsgruppe | block statt nur versuchsgruppe hinter die Tilde geschrieben. Dann wollen wir noch das compact letter disply und alles ist wie es sein sollte.
basi_rlm_fit |>
emmeans(specs = ~ versuchsgruppe) |>
cld(Letters = letters, adjust = "none") versuchsgruppe emmean SE df asymp.LCL asymp.UCL .group
Erde+Fließ 19.7 0.668 NA 18.4 21.0 a
Erde 19.9 0.668 NA 18.6 21.2 a
Perlite+Fließ 22.2 0.668 NA 20.9 23.5 b
Erde+Perlite 22.4 0.668 NA 21.1 23.7 b
Results are averaged over the levels of: block
Confidence level used: 0.95
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. Wir sehen also, dass sich Erde+Fließ und Erde nicht unterscheiden. Die beiden aber dann von Perlite+Fließ sowie Erde+Perlite. Am Ende unterscheiden sich Perlite+Fließ und Erde+Perlite nicht.
Für die Funktion robustbase::lmrob() müssen wir dann auf das R Paket {multcomp} umswitchen. Hier kriegen wir dann ein Problem, wenn wir eine Interaktion vorliegen ahben, da ist {multcomp} nicht so schön anzuwenden, wie dann {emmeans}. Aber hier schauen wir uns dann mal die einfache Implementierung ohne Interaktion an. Bei einer Interaktion müssten wir dann händisch über die Funktion filter() die Analysen für die einzelnen Blöcke aufspalten.
Wir lassen jetzt aber erstmal auf dem Modell aus der Funktion robustbase::lmrob() den paarweisen Vergleich nach Tukey laufen.
mult_lmrob <- glht(basi_lmrob_fit, linfct = mcp(versuchsgruppe = "Tukey"))Das erhaltende Objekt können wir dann mit tidy() aufräumen und uns dann die gewünschten Spalten wiedergeben lassen.
mult_lmrob_tbl <- mult_lmrob |>
tidy(conf.int = TRUE) |>
select(contrast, estimate, conf.low, conf.high, adj.p.value)
mult_lmrob_tbl# A tibble: 6 x 5
contrast estimate conf.low conf.high adj.p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 Erde+Fließ - Erde 1.68 -1.78 5.13 0.592
2 Erde+Perlite - Erde 1.63 -2.38 5.64 0.717
3 Perlite+Fließ - Erde 4.84 -1.35 11.0 0.181
4 Erde+Perlite - Erde+Fließ -0.0439 -2.86 2.78 1.00
5 Perlite+Fließ - Erde+Fließ 3.17 -2.52 8.85 0.474
6 Perlite+Fließ - Erde+Perlite 3.21 -2.79 9.21 0.509Wir sehen, dass wir mit der Adjustierung keinen signifikanten Unterschied zwischen den Bodensorten finden. Im Folgenden dann noch die Darstellung im compact letter display, wo wir etwas rumfrickeln müssen damit die Funktion multcompLetters() auch die Kontraste aus der Funktion glht() akzeptiert. Aber das haben wir dann alles zusammen.
mult_lmrob_tbl |>
mutate(contrast = str_replace_all(contrast, "\\s", "")) |>
pull(adj.p.value, contrast) |>
multcompLetters() |>
pluck("Letters") Erde+Fließ Erde+Perlite Perlite+Fließ Erde
"a" "a" "a" "a" Auch hier sehen wir keinen Unterschied im compact letter display für die adjustierten \(p\)-Werte. Für die unadjustierten rohen \(p\)-Werte nutzen wir dann den Umweg über die summary() Funktion. Die unadjustierten \(p\)-Werte können wir dann auch in das compact letter display umwandeln lassen.
summary(mult_lmrob, test = adjusted("none")) |>
tidy() |>
mutate(contrast = str_replace_all(contrast, "\\s", "")) |>
pull(p.value, contrast) |>
multcompLetters() |>
pluck("Letters") Erde+Fließ Erde+Perlite Perlite+Fließ Erde
"ab" "ab" "a" "b" Und wenn man eine andere Funktion nutzt, dann kommen auch leicht andere Ergebnisse raus. Bei den rohen \(p\)-Werten haben wir jetzt etwas andere Unterschiede. Aber das liegt auch an den sehr ähnlichen Mittelwerten. So groß sind die Unterschiede nicht, so dass wir hier natürlich durch unterschiedliche Algorithmen leicht andere Ergebnisse kriegen. Wieder ein Grund die \(p\)-Werte zu adjustieren um dann auch konsistente Ergebnisse zu haben.
56.5 Quantilsregression
Jetzt wollen wir uns als Spezialfall der robusten Regression einmal die Quantilesregression anschauen. Wir der Name schon sagt, optimiert die Quantilregression nicht über den Mittelwert und die Standardabweichung als quadratische Abweichungen sondern über den Median und die absoluten Abweichungen. Wir haben es also hier mit einer Alternativen zu der gewöhnlichen Regression zu tun. Das schöne an der Quantilesregression ist, dass der Median nicht so anfällig gegenüber Ausreißern ist. Hier haben wir einen Vorteil. Wir können theoretisch auch andere Quantile als Referenzwert nehmen und die absoluten Abstände berechnen. Wenn wir aber später dann mit emmeans die Gruppenvergleiche rechnen wollen, geht das nur über den Median.
Ein englisches Tutorium gibt es dann nochmal ausführlicher unter Quantile Regression as an useful Alternative for Ordinary Linear Regression
Um die Quantilesregression zu rechnen nutzen wir das R Paket {quantreg} und die Funktion rq(). Wichtig ist hier, dass wir das tau auf \(0.5\) setzen und somit dann eine Median-Regression rechnen. Wir können theoretisch auch andere Quantile als Referenz wählen, aber dann passt es nicht mehr mit emmeans(). Wir nutzen hier auch nochmal die Gelegenheit zwei Modelle mit anzupassen. Wir rechnen einmal ein Modell mit einer linearen Komponente der Zeit time_num und einmal ein Modell mit der quadratischen Komponente der Zeit mit poly(time_num, 2) was nicht anderes als time_num\(^2\) ist. Neben dieser Variante ist mit nlrq() auch eine Möglichkeit implementiert um noch komplexere nicht lineare Zusammenhänge zu modellieren. Hier dann einfach mal die Hilfeseite von ?nlrq anschauen.
time_rq_lin_fit <- rq(length ~ versuchsgruppe + time_num + versuchsgruppe:time_num, tau = 0.5,
gurke_time_len_tbl)
time_rq_quad_fit <- rq(length ~ versuchsgruppe * poly(time_num, 2), tau = 0.5,
gurke_time_len_tbl)In der Abbildung 56.5 plotten wir einmal die Modelle aus der linearen und der quadratischen Anpassung des zeitlichen Verlaufs in die Abbildung. Wir sehen, dass das Modell mit der quadratischen Anpassung besser zu den Daten passt. Daher nutzen wir jetzt im weiteren das Modell time_rq_quad_fit.
ggplot(gurke_time_len_tbl, aes(time_num, length, color = versuchsgruppe)) +
theme_minimal() +
geom_point2(position = position_dodge(0.5)) +
scale_color_okabeito() +
geom_line(aes(y = predict(time_rq_quad_fit), linetype = "Quadratic")) +
geom_line(aes(y = predict(time_rq_lin_fit), linetype = "Linear")) +
scale_x_continuous(breaks = 1:14) +
labs(linetype = "", color = "")
56.5.1 ANOVA
Machen wir es kurz. Die ANOVA ist für die Quantilsregression nicht implementiert und damit auch nicht anwendbar. Damit können wir aber leben, wenn wir die ANOVA nur als Vortest ansehen. Wir müssen dann eine mögliche Interaktion visuell überprüfen. Um die Interaktion visuell zu überprüfen nutzen wir die Funktion emimp() aus dem R Paket {emmeans}. Wir sehen in der Abbildung 56.6, dass sich die Gerade für die Versuchsgruppen über die Zeit nicht schneiden, was gegen eine starke Interaktion spricht. Die Steigung ist aber für alle drei Versuchsgruppen über die Zeit nicht gleich, wir haben zumindest eine mittlere Interaktion.
emmip(time_rq_lin_fit, versuchsgruppe ~ time_num, CIs = TRUE,
cov.reduce = FALSE) +
theme_minimal() +
scale_color_okabeito()
emmip(time_rq_quad_fit, versuchsgruppe ~ time_num, CIs = TRUE,
cov.reduce = FALSE) +
theme_minimal() +
scale_color_okabeito()
Wir würden also in unseren Gruppenvergleich der Mediane auf jeden Fall die Interaktion mit rein nehmen und die emmeans() Funktion entsprechend mit dem | anpassen.
56.5.2 Gruppenvergleich
Wenn wir den Gruppenvergleich in emmeans() rechnen wollen, dann geht es nur dem Spezialfall der Median-Regression. Wir müssen also zwangsweise in der Quantilesregression rq() das tau = 0.5 setzen um dann eine Median-Regression zu rechnen. Sonst können wir nicht emmeans nutzen.
time_rq_quad_fit |>
emmeans(~ versuchsgruppe | time_num,
adjust = "none", at = list(time_num = c(1, 7, 14))) |>
cld(Letters = letters, adjust = "none") time_num = 1:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 2.30 0.0513 243 2.20 2.40 a
Katrina 3.30 0.2879 243 2.73 3.87 b
Proloog 5.00 0.2243 243 4.56 5.44 c
time_num = 7:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 3.83 0.2467 243 3.35 4.32 a
Katrina 5.20 0.3184 243 4.57 5.83 b
Proloog 11.97 0.2938 243 11.39 12.55 c
time_num = 14:
versuchsgruppe emmean SE df lower.CL upper.CL .group
Quarto 7.30 1.7815 243 3.79 10.81 a
Katrina 11.36 3.1038 243 5.25 17.47 a
Proloog 30.06 1.6609 243 26.79 33.34 b
Confidence level used: 0.95
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. Das sieht ja schon mal ganz gut aus. Interessant ist, dass wir an dann in der vierzehnten Woche dann keinen signifikanten Unterschied mehr vorliegen zwischen Quarto und Katrina vorliegen haben. Das liegt dann aber alleinig an dem hohen Standardfehler SE der Sorte Katrina mit \(3.1\) gegenüber den anderen Sorten. Vermutlich sehe das Ergebnis leicht anders aus, wenn wir die unsinnigen Gurken mit einer Länge von 0 cm wirklich aus den Dten entfernen würden und nicht zu Demonstrationszwecken wie hier drin lassen würden.
Aber Achtung, die Spalte emmean beschriebt hier die Mediane. Wir haben also hier die Mediane des Längenwachstums für die Gurken vorliegen. Wenn wir die paarweisen Vergleich rechnen wollen würden dann können wir noch pairwise ~ versuchsgruppe statt ~ versuchsgruppe schreiben. Auch liefert die Funktion pwpm() die Medianunterschiede wie wir im Folgenden einmal sehen. Ich habe hier mal das at entfernt und da die Zeit numrisch ist, haben wir dann auch nur noch eine Tabelle. Die Diagonal: [Estimates] (emmean) sind die Mediane des Längenwachstums.
time_rq_quad_fit |>
emmeans(~ versuchsgruppe | time_num, adjust = "none") |>
pwpm() Katrina Proloog Quarto
Katrina [ 5.50] <.0001 0.0095
Proloog -7.41 [12.91] <.0001
Quarto 1.48 8.89 [ 4.02]
Row and column labels: versuchsgruppe
Upper triangle: P values adjust = "tukey"
Diagonal: [Estimates] (emmean)
Lower triangle: Comparisons (estimate) earlier vs. later56.6 Modellvergleich
Im Folgenden wollen wir einmal verschiedene Modelle miteinander vergleichen und schauen, welches Modell hier das beste ist. Das machen wir dann einmal für die Basilikumdaten sowie die Wachstumsdaten für die Gurken.
56.6.1 Basilikumdaten
Für den Modellvergleich der gewöhnlichen, robusten und medianen Regression nutzen wir nochmal den Datensatz für das Basilikumwachstum. In einem ersten Schritt fitten wir wieder alle Modelle und achten darauf, dass wir bei der Quantilesregression angeben welches Quantile wir wollen. Wir wählen mit tau = 0.5 dann den Median und rechnen so eine Medianregression.
basi_lm_fit <- lm(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)
basi_rlm_fit <- rlm(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)
basi_lmrob_fit <- lmrob(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl)
basi_rq_fit <- rq(weight ~ versuchsgruppe + block + versuchsgruppe:block, basi_block_tbl,
tau = 0.5)Im Folgenden nutzen wir dann das fantastische Paket {modelsummary} mit der gleichnamigen Funktion um uns einmal die Modelle im Vergleich anzuschauen. Hier hilft dann wie immer die tolle Hilfsseite von {modelsummary} zu besuchen. Ich möchte nur nicht den Intercept und die Schätzer für die Blöcke haben, deshalb fliegen die hier einmal raus.
modelsummary(lst("Ordinary" = basi_lm_fit,
"MASS::rlm" = basi_rlm_fit,
"robustbase::lmrob" = basi_lmrob_fit,
"Quantile" = basi_rq_fit),
estimate = "{estimate}",
statistic = c("conf.int",
"s.e. = {std.error}",
"t = {statistic}",
"p = {p.value}"),
coef_omit = "Intercept|block")| Ordinary | MASS::rlm | robustbase::lmrob | Quantile | |
|---|---|---|---|---|
| versuchsgruppeErde+Fließ | 2.200 | 1.515 | 1.675 | 1.000 |
| [-1.399, 5.799] | [-2.257, 5.286] | [-1.053, 4.403] | [-0.374, 2.374] | |
| s.e. = 1.801 | s.e. = 1.888 | s.e. = 1.366 | s.e. = 0.701 | |
| t = 1.221 | t = 0.802 | t = 1.227 | t = 1.426 | |
| p = 0.226 | p = 0.425 | p = 0.224 | p = 0.159 | |
| versuchsgruppeErde+Perlite | 2.200 | 1.515 | 1.631 | 0.000 |
| [-1.399, 5.799] | [-2.257, 5.286] | [-1.529, 4.791] | [-3.073, 3.073] | |
| s.e. = 1.801 | s.e. = 1.888 | s.e. = 1.582 | s.e. = 1.568 | |
| t = 1.221 | t = 0.802 | t = 1.031 | t = 0.000 | |
| p = 0.226 | p = 0.425 | p = 0.306 | p = 1.000 | |
| versuchsgruppePerlite+Fließ | 5.200 | 4.683 | 4.841 | 6.000 |
| [1.601, 8.799] | [0.911, 8.455] | [-0.037, 9.718] | [1.045, 10.955] | |
| s.e. = 1.801 | s.e. = 1.888 | s.e. = 2.442 | s.e. = 2.528 | |
| t = 2.887 | t = 2.480 | t = 1.982 | t = 2.373 | |
| p = 0.005 | p = 0.016 | p = 0.052 | p = 0.021 | |
| Num.Obs. | 80 | 80 | 80 | 80 |
| R2 | 0.535 | 0.501 | 0.480 | |
| R2 Adj. | 0.426 | 0.384 | ||
| AIC | 410.7 | 411.8 | 404.5 | |
| BIC | 451.1 | 452.3 | 442.7 | |
| Log.Lik. | -188.326 | -188.914 | ||
| F | 4.912 | 3.903 | ||
| RMSE | 2.55 | 2.57 | 2.55 | 2.69 |
Wie wir sehen sehen haben nicht alle Modelle die gleichen Informationen zurück gegeben. Insbesondere das fehlen des Bestimmtheitsmaßes \(R^2\) bei der robusten Regression MASS::rlm ist schmerzlich, da wir hier dann nicht die Möglichkeit haben eine Aussage über die erklärte Varianz der robusten Regression zu treffen. Auch fehlt bei der Medianregression das adjustierte \(R^2\), was die Nutzung bei Modellen mit mehr als einer Einflussvariable \(x\) im Modell erschwert bis nutzlos macht. Da ist dann die Funktion aus dem Paket {robustbase} besser, hier haben wir dann ein \(R^2\) vorliegen. Daher bleibt uns am Ende nur das AIC oder BIC, wobei wir dort den kleinsten Wert als besser erachten. Die AIC und BIC Werte sind somit am besten für die Quantilesregression. Für die Funktion robustbase::lmrob haben wir dann kein AIC. Dafür ist dann aber der Fehler RMSE bei der gewöhnlichen Regressionzusammen mit der Funktion robustbase::lmrob am niedrigsten. Und so stehe ich wieder davor und weiß nicht was das beste Modell ist. Hier müsste ich dann nochmal überlegen, ob ich lieber über Mittelwerte oder Mediane berichten möchte und das ist ohne die Forschungsfrage nicht hier zu lösen.
56.6.2 Gurkendaten
Für den Modellvergleich der gewöhnlichen, robusten und medianen Regression nutzen wir nochmal den Datensatz für das Längenwachstum der Gurken. In einem ersten Schritt fitten wir wieder alle Modelle und achten darauf, dass wir bei der Quantilesregression angeben welches Quantile wir wollen. Wir wählen mit tau = 0.5 dann den Median und rechnen so eine Median-Regression. Darüber hinaus schauen wir uns nochmal die quadratische Anpassung der Zeit in dem letzten Modell an. Mal schauen, ob das Modell dann auch bei den Kennzahlen das beste Modell ist. Visuell sah das quadratische Modell des zeitlichen Verlaufs schon mal sehr gut aus.
time_lm_fit <- lm(length ~ versuchsgruppe + time_num + versuchsgruppe:time_num, gurke_time_len_tbl)
time_rlm_fit <- rlm(length ~ versuchsgruppe + time_num + versuchsgruppe:time_num, gurke_time_len_tbl)
time_rq_lin_fit <- rq(length ~ versuchsgruppe + time_num + versuchsgruppe:time_num, tau = 0.5,
gurke_time_len_tbl)
time_rq_quad_fit <- rq(length ~ versuchsgruppe * poly(time_num, 2), tau = 0.5,
gurke_time_len_tbl)Auch hier nutzen wir dann das fantastische Paket {modelsummary} mit der gleichnamigen Funktion um uns einmal die Modelle im Vergleich anzuschauen.
modelsummary(lst("Ordinary" = time_lm_fit,
"Robust" = time_rlm_fit,
"Quantile linear" = time_rq_lin_fit,
"Quantile quadratic" = time_rq_quad_fit),
estimate = "{estimate}",
statistic = c("conf.int",
"s.e. = {std.error}",
"t = {statistic}",
"p = {p.value}"),
coef_omit = "Intercept|time_num")| Ordinary | Robust | Quantile linear | Quantile quadratic | |
|---|---|---|---|---|
| versuchsgruppeProloog | -0.218 | -1.028 | -0.823 | 8.483 |
| [-3.056, 2.619] | [-2.693, 0.638] | [-2.852, 1.206] | [6.612, 10.353] | |
| s.e. = 1.440 | s.e. = 0.846 | s.e. = 1.035 | s.e. = 0.954 | |
| t = -0.152 | t = -1.215 | t = -0.795 | t = 8.890 | |
| p = 0.880 | p = 0.225 | p = 0.427 | p = <0.001 | |
| versuchsgruppeQuarto | -0.412 | -0.641 | -0.986 | -1.883 |
| [-3.249, 2.426] | [-2.307, 1.025] | [-3.045, 1.074] | [-3.817, 0.050] | |
| s.e. = 1.440 | s.e. = 0.846 | s.e. = 1.051 | s.e. = 0.986 | |
| t = -0.286 | t = -0.758 | t = -0.938 | t = -1.909 | |
| p = 0.775 | p = 0.449 | p = 0.349 | p = 0.057 | |
| Num.Obs. | 252 | 252 | 252 | 252 |
| R2 | 0.581 | 0.548 | 0.537 | |
| R2 Adj. | 0.572 | |||
| AIC | 1472.0 | 1488.9 | 1345.9 | 1316.6 |
| BIC | 1496.7 | 1513.6 | 1367.1 | 1348.4 |
| Log.Lik. | -729.020 | -737.435 | ||
| F | 68.209 | 269.196 | ||
| RMSE | 4.37 | 4.51 | 4.53 | 4.59 |
Hier sieht es ähnlich aus wie bei den Modellvergleichen von den Basilikumdaten. Nur hier wissen wir, dass wir Ausreißer in der Form von Gurken mit einer Länge von 0 cm in den Daten haben. Daher sehen wir auch, dass die robuste Regression und die Median-Regression ein niedrigeres AIC und BIC haben. Für mich interessant ist, dass der Fehler RMSE wieder am kleinsten bei der gewöhnlichen Regression ist, das mag aber auch an der Berechnung liegen. Da wir Ausreißer haben, sind natürlich dann auch die robuste Regression und die Median-Regression vorzuziehen. Da die Median-Regression den kleineren AIC Wert hat, nehmen wir dann die Median-Regression als die beste Regression an. Der RMSE ist am kleinsten für die quadratische Anpassung des zeitlichen Verlaufs. Deshalb wäre dann das vierte Modell mit dem niedrigsten AIC und dem niedrigsten RMSE das Modell der Wahl.
Ich wiederhole mich hier wirklich, aber vermutlich wäre es schlauer zuerst die Daten zu bereinigen und die Gurken mit einem Wachstum von 0 cm zu entfernen. Auch die anderen Wachstumskurven von anderen Gurken sind etwas wirr. Da müsste man auch nochmal ran und schauen, ob nicht die Gurken lieber aus der Analyse raus müssen. Am Ende ist dann natürlich die Frage, ob man die Daten dann nicht doch lieber über ein gemischtes Modell auswertet, da natürlich die Zeitpunkte voneinander abhängig sind, wenn wir immer die glecihen Gurken messen und nicht ernten. Aber wie immer im Leben, alles geht nicht.
Referenzen
Maronna RA, Martin RD, Yohai VJ, Salibián-Barrera M. 2019. Robust statistics: theory and methods (with R). John Wiley & Sons.
Stromberg A. 2004. Why write statistical software? The case of robust statistical methods. Journal of Statistical Software 10: 1–8.