53  Marginal effect models

Letzte Änderung am 10. July 2025 um 13:07:11

“Life is difficult.” — Morgan Scott Peck, The Road Less Traveled

Stand des Kapitels: Konstruktion (seit 07.2025)

Dieses Kapitel wird in den nächsten Wochen geschrieben und ist damit meine aktuelle Großbaustelle. Ich plane zum Beginn des WiSe 2025/26 eine fertige Version des Kapitels erstellt zu haben. Während das Kapitel entsteht, funktioniert so manches dann nicht so wie es soll. Bitte daher hier dann abwarten.

Lange habe ich gebraucht um mich dazu durchzuringen das Kapitel zu den Marginal effect models (deu. marginale Effektmodelle, ungebräuchlich) zu schreiben. Ich werde hier bei dem englischen Begriff bleiben, den deutschen Begriff habe ich eher selten gehört und daher sind es für mich Marginal effect models. Insbesondere da der Begriff “marginal” auch sehr an gering oder minderwertig erinnert. Damit haben aber die Marginal effect models nicht im geringsten zu tun. Die Modelle sind sehr mächtig und können uns helfen wichtige Fragen an unsere Daten zu beantworten. Insbesondere die Dualität der beiden Pakete {emmeans} für experimentelle Daten und {marginaleffects} für beobachtete Daten ist spannend und möchte ich hier nochmal genauer betrachten. Neben diesen beiden Ecksteinen gibt es noch andere Pakete und ich werde auch hier einmal in die Pakete reinschauen.

Anfangen kann ich aber nicht ohne Heiss (2022) mit seinem Blogpost Marginalia: A guide to figuring out what the heck marginal effects, marginal slopes, average marginal effects, marginal effects at the mean, and all these other marginal things are zu erwähnen. Vieles entlehnt sich direkt oder indirekt an seine Ausführungen. Wie immer habe ich etwas angepasst, wenn ich der Meinung war, dass es noch besser zu verstehen ist. Teilweise entfallten die Marginal effect models ihre wahre Kraft erst bei den nicht linearen Zusammenhängen und der Interpretation von Generalized Additive Models udn somit der nicht linearen Regression. Fangen wir also an Marginal effect models zu zerforschen und arbeiten uns dann voran.. Beginnen wollen wir aber mit einem allgemeinen Hintergrund bevor wir uns dann nochmal tiefer mit den Marginal effect models beschäftigen.

53.1 Allgemeiner Hintergrund

“Statistics is all about lines, and lines have slopes, or derivatives. These slopes represent the marginal changes in an outcome. As you move an independent/explanatory variable, what happens to the dependent/outcome variable?” — Heiss (2022)

Wenn wir von Marginal effect models sprechen, dann können wir uns im Prinzip zwei Aspekte anschauen. Wir können über die Steigung einer Funktion einer Geraden sprechen oder aber über die vorhergesagten \(y\)-Werte auf der Geraden für beliebige \(x\)-Werte. Damit sind wir dann bei den beiden Aspekten Steigung und Vorhersage. Wenn wir uns in der Welt der linearen Modelle bewegen, dann ist die Steigung eigentlich kein Problem und die Vorhersage auch nicht so komplex. Spannender wird das Zusammenspiel eines nichtlinearen Modells und eben den Marginal effect models. Hier kommt dann die eigentliche Kraft der Marginal effect models zu trage. In den folgenden beiden Abbildungen habe ich dir einmal eine nichtlinere Funktion dargestellt. Wir schauen uns hier den Zusammenhang zwischen der standardisierten Enzymeaktivität und dem standardisierten pH-Wert an. Wir haben die Enzymeaktivität zu festen pH-Werten wiederholt gemessen. Auf der linken Seite betrachten wir die Steigung der Geraden an drei Punkten und auf der rechten Seite sehen wir einmal die Vorhersage für drei pH-Werte auf der Geraden. Die Gerade folgt der Funktion \(y = x^3 - 8x^2 + 10x + 10\) und ist mir somit für die Bestimmung der Steigung und der Vorhersage bekannt. Daher habe ich also die Möglichkeit die exakten Werte der Steigung und der Vorhersage zu bestimmen. Einen Luxus den wir selten mit echten Daten haben.

Abbildung 53.1— Modellierung des nichtlinearen Zusammenhangs zwischen der standardisierten Enzymeaktivität und dem standardisierten pH-Wert. Der pH-Wert ist kontinuierlich. Die Funktionsgleichung ist bekannt. (A) Berechnung der Steigung (eng. slope) für ausgewählte pH Werte anhand der Modellierung. (B) Berechnete Vorhersagewerte (eng. prediction) der standardisierten Enzymeaktivität für ausgewählte pH-Werte anhand der Modellierung. [Zum Vergrößern anklicken]

Betrachten wir also einmal die Antworten, die die Steigung und die Vorhersage liefert. Dabei haben wir hier in diesem Fall ein kontinuierliches \(X\) vorliegen und machen uns die Sachlage auch einfacher indem wir ein kontinuierliches \(Y\) mit der standardisierten Enzymeaktivität vorliegen haben.

Welche Antwort liefert die Steigung?

Wenn sich \(X\) ändert, wie ändert sich dann \(Y\) an dem Wert von \(X\)? Hierbei muss sich \(X\) nicht um eine Einheit verändern, wie wir es gerne im linearen Zusammenhang sagen, sondern wir wollen die Steigung direkt im Punkt von \((X|Y)\) haben.

Welche Antowort liefert die Vorhersage?

Welche Werte für \(Y\) sagt das Modell für \(X\) vorraus? Wir müssen hier die beobachteten Werte von \(Y\), die in unserem Beispiel für einen pH-Wert wiederholt vorkommen, von dem einen vorhergesagten \(Y\) Wert aus dem Modell für ein gegebenes \(X\) unterscheiden.

Dann können wir auch schon die Steigung und die Vorhersage einmal interpretieren. In dem linken Tab findest du einmal die Interpretation der Steigung sowie die Ausgabe der Funktion slopes() aus dem R Paket {marginaleffects}. In dem rechten Tab dann die Ergebnisse der Vorhersage und die Ausgabe der Funktion predictions(). Mehr zu den beiden Funktionen dann weiter unten in der Anwendung. Ich berechne hier die Steigung und bestimme die vorhergesagten Werte für drei ausgewählte pH-Werte.

Interpretation Steigung

Wir können einmal die Steigung (dx/dy) aus der ersten Ableitung der quadratischen Geradenfunktion für unsere ausgewählten pH-Werte bestimmen und dann entsprechend interpretieren.

Tabelle 53.1— Interpretation der Steigung der Enzymeaktivität an drei ausgewählten pH-Werten.

pH	Steigung	Interpretation
-1	29	Bei einem standardisierten pH-Wert von -1 steigt die Enzymeaktivität um 29U an.
2	-10	Bei einem standardisierten pH-Wert von 2 sinkt die Enzymeaktivität um 10U.
6	22	Bei einem standardisierten pH-Wert von 6 steigt die Enzymeaktivität um 22U an.

Wir können uns dann die Steigung auch direkt mit der Funktion slopes() bestimmen lassen und erhalten dann die folgenden Informationen. Häufig haben wir ja nicht die Geradengleichung vorliegen. Hier haben wir dann auch die p-Werte sowie einen entsprechenden Fehler. Wir haben hier eine leichte Abweichung, da ich die obige Steigung durch die Ableitung der quadratischen Funktion erstellt habe und nicht aus einem Polynomialmodell entnommen habe.

 ph Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
 -1     28.3       2.44 11.57   <0.001 100.5  23.5  33.04
  2    -10.3       1.14 -9.08   <0.001  63.0 -12.6  -8.11
  6     21.4       1.63 13.16   <0.001 129.0  18.3  24.64

Term: ph
Type: response
Comparison: dY/dX

Interpretation Vorhersage

Auch hier können wir direkt durch das Einsetzen der pH-Werte in unsere quadratische Geradenfunktion die vorgesagten Enzymeaktivitäten für die ausgewählten pH-Werte bestimmen.

Tabelle 53.2— Interpretation der Vorhersage der Enzymeaktivität an drei ausgewählten pH-Werten.

pH	Vorhersage (y)	Interpretation
-1	-9	Für einen standardisierten pH-Wert von -1 sagt die Funktion eine Enzymeaktivität von -9U vorher.
2	6	Für einen standardisierten pH-Wert von 2 sagt die Funktion eine Enzymeaktivität von 6U vorher.
6	-2	Für einen standardisierten pH-Wert von 6 sagt die Funktion eine Enzymeaktivität von -2U vorher.

Auch können wir die Steigung aus den Daten direkt mit der Funktion predictions() bestimmen. Häufig haben wir ja nicht die Geradengleichung vorliegen. Hier haben wir theoretisch noch eine riesige Auswahl an Funktionen in R, wir konzentrieren uns hier aber auf das R Paket {marginaleffects}. Wir haben auch hier eine leichte Abweichung, da ich die obigen Vorhersagen durch die quadratische Funktion bestimmt habe und nicht aus einem Polynomialmodell entnommen habe.

 ph Estimate Std. Error      z Pr(>|z|)    S 2.5 % 97.5 %
 -1   -4.545       2.48 -1.830   0.0673  3.9 -9.41  0.324
  2    8.951       2.23  4.013   <0.001 14.0  4.58 13.323
  6   -0.545       2.36 -0.231   0.8176  0.3 -5.18  4.089

Type: response

Wenn wir die Vorhersage betrachten dann können wir auch den Fall vorliegen haben, dass wir ein kategoriales \(X\) als Gruppen gemessen haben. In R wäre es dann ein Faktor und daher dann auch der Begriff des faktoriellen Designs, wenn es um das experimentelle Design geht. Wir hatten bis jetzt ein kontinuierlichen pH-Wert vorliegen. Zwar hatten wir auch Messgruppen, aber wir konnten von einem konitnuierlichen pH-Wert im Sinne der Modellierung ausgehen. Jetzt wollen wir uns einmal den Fall anschauen, dass wir auf eben wirkliche pH-Wertgruppen mit niedrigen, mittleren und hohen pH-Wertgruppen für die Enzymeaktivität vergleichen wollen. Ich schreibe hier schon gleiche von einem Vergleich, als erstes wollen wir aber die Mittelwerte pro Gruppe vorhersagen.

In der folgenden Abbildung siehst du auf der linken Seite einmal den einfaktoriellen Fall mit dem gruppierten pH-Wert als Gruppe. Auf der rechten Seite dann noch zusätzlich die beiden Gruppen der Prokaryoten sowie Eukaryoten. Wir sind jetzt an den Mittelwerten pro Gruppe interessiert und wir könnten diese Mittelwerte dann auch einfach berechnen. Dafür bräuchten wir dann erstmal kein Modell, wenn wir nur so ein simples Experiment vorliegen haben. Konkret bestimmen wir hier die marginal means aus einem Modell heraus.

Abbildung 53.2— Modellierung des Zusammenhangs zwischen der standardisierten Enzymeaktivität und den gruppierten pH-Werten nach niedrigen, mittleren und hohen pH-Werten. Der pH-Wert ist kategorial. (A) Einfaktorielle Vorhersage der Gruppenmittelwerte der Enzymeaktivität. (B) Zweifaktorielle Vorhersage der Gruppenmittelwerte der Enzymeaktivität aufgetrennt nach der Gruppe der Eukaryoten und Prokaryoten. [Zum Vergrößern anklicken]

Wir wollen jetzt für alle Gruppen die Mittelwerte vorhersagen. Nichts anderes macht dann auch die Funktion predictions(). Hier muss ich auch schon gleich einmal eine leichte Warnung aussprechen. So gut {marginaleffects} in der Bestimmung der Steigung und der Vorhersagen ist, um so viel besser ist das R Paket {emmeans} wenn es um die Auswertung von einem faktoriellen Design geht. Hier ist dann einfach {emmeans} besser in der Anwendung. Das hat auch Gründe im Algorithmus der beiden Pakete, dazu dann aber später mehr. Hier erstmal die Interpretation der beiden Vorhersagen für kategoriale \(x\)-Werte oder eben auch Faktoren in R genannt.

Vorhersage 1-faktoriell

Wir haben hier also in unserem einfaktoriellen Modell einmal die gruppierten pH-Werte in den Gruppen niedrig, mittel und hoch vorliegen. Wir wollen dann den Mittelwert für jede der Gruppen bestimmen.

Tabelle 53.3— Interpretation der Vorhersage der Enzymeaktivität an drei ausgewählten pH-Werten.

pH	Mittelwert	Interpretation
niedrig	3.2	Für die Gruppe der niedrigen pH-Werte haben wir im Mittel eine Enzymeaktivität von 3.2U vorliegen.
mittel	2.59	Für die Gruppe der mittleren pH-Werte haben wir im Durchschnitt eine Enzymeaktivität von 2.59U vorliegen.
hoch	-5.6	Für die Gruppe der hohen pH-Werte haben wir im Mittel eine Enzymeaktivität von -5.9U vorliegen.

Jetzt können wir uns auch mit der Funktion predictions() aus dem R Paket {marginaleffects} die Mittelwerte für die Gruppen vorhersagen lassen. Wir erhalten dann auch die entsprechenden Standardfehler und andere statistische Maßzahlen. Hier kriegen wir noch keinen Vergleich, hier wird nur getestet, ob der Mittelwert sich von der Null unterscheidet.

     grp Estimate Std. Error      z Pr(>|z|)    S  2.5 % 97.5 %
 niedrig    -2.50       3.99 -0.626    0.531  0.9 -10.32  5.327
 mittel     -1.35       1.06 -1.276    0.202  2.3  -3.42  0.722
 hoch       17.11       4.91  3.483   <0.001 11.0   7.48 26.745

Type: response

Vorhersage 2-faktoriell

Kommen wir nun zu der Vorhersage der Mittelwerte oder Marginal means in den gruppierten pH-Werten aufgeteilt nach den beiden Gruppen der Prokaryoten sowie Eukaryonten. Daher haben wir hier ein zweifaktorielles Design vorliegen. Wir wollen eben die Mittelwerte für alle Faktorkombinationen wissen.

Tabelle 53.4— Interpretation der Vorhersage der Enzymeaktivität an drei ausgewählten pH-Werten in den Prokaryoten und Eukaryonten

pH	Gruppe	Mittelwert	Interpretation
niedrig	Eukaryot	14.21	Für die Gruppe der niedrigen pH-Werte der Eukaryoten haben wir im Mittel eine Enzymeaktivität von 14.21U vorliegen.
niedrig	Prokaryot	-19.2	Für die Gruppe der niedrigen pH-Werte der Prokaryoten haben wir im Mittel eine Enzymeaktivität von -19.2U vorliegen.
mittel	Eukaryot	-9.22	Für die Gruppe der mittleren pH-Werte der Eukaryoten haben wir im Mittel eine Enzymeaktivität von -9.22U vorliegen.
mittel	Prokaryot	2.59	Für die Gruppe der mittleren pH-Werte der Prokaryoten haben wir im Mittel eine Enzymeaktivität von 2.59U vorliegen.
hoch	Eukaryot	26.24	Für die Gruppe der hohen pH-Werte der Eukaryoten haben wir im Mittel eine Enzymeaktivität von 26.24U vorliegen.
hoch	Prokaryot	-5.14	Für die Gruppe der hohen pH-Werte der Prokaryoten haben wir im Mittel eine Enzymeaktivität von -5.14U vorliegen.

Wenn wir dann die Funktion predictions() aus dem R Paket {marginaleffects} benutzen erhalten wir dann auch die Marginal means für alle Faktorkombinationen zurück. Hier sind die Werte die gleichen wir auch in der simplen Berechnung oben in der Abbildung, da unser Modell eben dann doch nur die beiden Faktoren enthält. Wir erhalten dann auch die Standardfehler und den p-Wert für den Test gegen einen Mittelwert von Null.

     grp      type Estimate Std. Error       z Pr(>|z|)     S   2.5 % 97.5 %
 niedrig Eukaryot     14.21     0.5498   25.84   <0.001 486.6  13.129  15.28
 niedrig Prokaryot   -19.20     4.1190   -4.66   <0.001  18.3 -27.277 -11.13
 mittel  Eukaryot     -9.22     0.0325 -283.79   <0.001   Inf  -9.286  -9.16
 mittel  Prokaryot     2.59     1.2580    2.06   0.0396   4.7   0.123   5.05
 hoch    Eukaryot     28.24    13.4263    2.10   0.0354   4.8   1.924  54.55
 hoch    Prokaryot    -5.14     0.5571   -9.22   <0.001  64.9  -6.228  -4.04

Type: response

Soviel dann einmal zu dem allgemeinen Hintergrund. Die Idee der Marginal effect models ist es also dir die Steigung einer Funktion oder die vorhergesagten Werte einer Funktion wiederzugeben. Hier kannst du dann Funktion auch mit Modell erstetzen. Für den linearen Fall auf einem normalverteilten Messwert ist die Anwendung und der Erkenntnisgewinn der Marginal effect models begrenzt. Da reichen dann auch die Ausgaben der Standardfunktionen. Aber auch dort werden wir dann noch in den entsprechenden Kapiteln sehen, dass wir hier noch was rausholen können. Die Marginal effect models entwicklen dann ihre Stärke für nicht normalverteilte Messwerte sowie eben nicht lineare Zusammenhänge.

Weitere Tutorien für die Marginal effects models

Wir oben schon erwähnt, kann dieses Kapitel nicht alle Themen der Marginal effects models abarbeiten. Daher präsentiere ich hier eine Liste von Literatur und Links, die mich für dieses Kapitel hier inspiriert haben. Nicht alles habe ich genutzt, aber vielleicht ist für dich was dabei.

• Ohne den Blogpost Marginalia: A guide to figuring out what the heck marginal effects, marginal slopes, average marginal effects, marginal effects at the mean, and all these other marginal things are von Heiss (2022) wäre dieses Kapitel nicht möglich gewesen.
• Wie alles im Leben ist nichts ohne Kritik. Is least squares means (lsmeans) statistical nonsense? ist dann auch eine gute Frage. Ich bin der Meinung nein und auch andere sind es, aber hier kannst du dann nochmal eine andere Meinung lesen.

53.2 Auf ein Wort…

“What he says?” — Asterix, Sieg über Caesar

Soweit so gut. Wenn du verstanden hast, was die Marginal effect models können, dann kannst du auch bei den Daten und deren Auswertung weitermachen. Hier geht es dann etwas tiefer und ich gehe nochmal auf einzelen Aspekte etwas ausführtlicher ein. So wollen wir nochmal verstehen, was eigentlich die Steigung nochmal war und wie wir die Steigung berechnen. Dann müssen wir nochmal über simple und multiple Modelle sprechen. Dann gehen wir nochmal auf die Vorhersage ein und ich gebe nochmal einen kurzen Überblick, was wir da eigentlich alles vorhersagen oder genauer was wie heißt. Dann besprechen wir nochmal die Unterschiede zwischen den Paketen {marginaleffects} und {emmeans}. Wie immer kann man den Teil hier auch überspringen, wenn es nur um die Anwendung geht.

53.2.1 …zur Steigung

Wenn wir mit dem Verstehen und Zerforschen der Steigung vorankommen wollen, dann können wir Heiss (2022) und Arel-Bundock et al. (2024) mit der Veröffentlichung Model to meaning — How to Interpret Statistical Models With marginaleffects for R and Python nicht ignorieren. Ich nutze jetzt eine etwas allgemeinere Erklärung der Marginal effect models und konzentriere mich erstmal auf ein normalverteiltes \(x\) sowie ein normalverteiltes \(y\) in einem simplen Modell mit einem \(x\) und einem \(y\) in der folgenden Form.

\[ y \sim \beta_0 + \beta_1 x_1 \] mit

• \(\beta_0\), dem Koeffizienten des y-Achsenabschnitt oder Intercept der Geraden.
• \(\beta_1\), dem Koeffizienten der Steigung der Geraden

Daher haben wir hier in unserem \(x\) keine Gruppen vorliegen sondern einen klassischen Scatterplot mit Punkten als Beobachtungen. Wir können die Marginal effect models auch auf beliebige kategorielle \(x\) wie eben Behandlungsgruppen sowie jedes beliebige \(y\) anwenden, aber hier fangen wir einmal einfach an.

Welche Frage wollen wir mit Marginal effect models beantworten?

Wenn sich das \(x\) um einen Wert oder eine Einheit erhöht, um wieviele Einheiten verändert sich dann der Wert von \(y\)?

In der folgenden Abbildung siehst du einmal zwei Scatterplots. In dem linken Scatterplot haben wir einen linearen Zusammenhang zwischen unseren \(x\)-Werten und den \(y\)-Werten. Wir können sagen, dass wenn sich \(x\) um einen Wert erhöht, dann erhöht sich auch \(y\) um einen konstanten Wert. Dieser konstante Wert um den sich die \(y\)-Werte mit ansteigenden \(x\) erhöhen, nennen wir auch die Steigung \(\beta_1\). In einem linearen Zusammenhang ist die Frage damit mit der Steigung der Geraden eigentlich beantwortet. Steigt \(x\) um einen Wert, dann steigt \(y\) um den Wert der Steigung \(\beta_1\) der Geraden. Diesen konstanten Zusammenhang haben wir aber nicht bei einem quadratischen Zusammenhang wie in der rechten Abbildung. Wir können hier nicht sagen, dass wenn sich \(x\) um einen Wert erhöht, sich auch \(y\) um einen konstanten Wert ändert. Hier hängt es von dem betrachteten \(x\)-Wert ab.

Abbildung 53.3— Scatterplot der kontinuierlichen x-Werte und kontinuierlichen y-Werte. In einem Modell wird die Abhängigkeit von y und x modellieren. (A) Linearer Zusammenhang. (B) Quadratischer Zusammenhang. [Zum Vergrößern anklicken]

Schauen wir mal in ein Zahlenbeispiel und lassen die Beobachtungen weg. Beginnen wir einmal mit dem linearen Zusammenhang der Funktion \(f(x) = 2x-1\). Ich habe die Gerade einmal in der folgenden Abbildung eingezeichnet. Wenn usn jetzt die Steigung an jedem beliebigen Punkt von \(x\) interessiert, dann bilden wir die erste Abbleitung \(f'(x) = 2\). Erhöht sich also der Wert von \(x\) um 1 dann steigt der Wert von \(y\) um 2 an. Wir sehen aber auch, dass für jedes beliebige Punktepaar wir eine Steigung von 2 vorliegen haben.

Abbildung 53.4— Gerade des Modells für einen linearen Zusammenhang. In einem Modell wird die Abhängigkeit von y und x modellieren. (A) Lineares Modell mit Gleichung. (B) Steigung an der Geraden für ausgewählte Punktepaare. [Zum Vergrößern anklicken]

Spanndender wird die Sachlage in einem quadratischen Zusammenhang in der folgenden Abbildung. Oder allgemeiner gesprochen, wenn wir keinen linearen Zusammenhang vorliegen haben. Wir haben hier den Zusammenhgang \(f(x) = -0.5x^2+5x\) vorliegen. Damit haben wir dann eine erste Ableitung von \(f'(x) = x+5\). Wie du siehst, ändert sich auch die Steigung in Abhänigkeit von \(x\). Wenn wir \(x\)-Werte links betrachten, dann liegt hier eher eine positive Steigung vor. Wenn wir nach rechts laufen, dann sehen wir immer stärkere negative Steigungen. Und hier kommen dann die Marginal effect models ins Spiel. Wir können allgemein gesprochen uns mit den Marginal effect models für jedes \(x\) die Steigung wiedergeben lassen.

Abbildung 53.5— Gerade des Modells für einen quadratischen Zusammenhang. In einem Modell wird die Abhängigkeit von y und x modellieren. (A) Quadratisches Modell mit Gleichung. (B) Steigung an der Geraden für ausgewählte Punktepaare. [Zum Vergrößern anklicken]

Aber moment, denkst du jetzt, in dem linearen Zusammenhang ist es ja einfach mit der Steigung für jeden beliebigen \(x\)-Wert. Wir erhalten für jeden \(x\)-Wert genau die gleiche Steigung. Aber bei den nicht-linearen Zusammenhängen hat ja jeder \(x\)-Wert seine eigene Steigung. Wenn wir viele \(x\)-Werte gemessen haben, dann haben wir ja dutzende bis hunderte Steigungen durch ein Marginal effect model ermittelt. Das stimmt und damit kommen wir auch gleich zu dem nächsten Punkt, dem Aggregieren der Daten. Oder wie im folgenden Cartoon richtig dargestellt, müssen wir uns überlegen wie wir den den Durchschnitt der Steigungen berechnen.

Abbildung 53.6— “Should I cut the red wire or the blue one!?” “WAIT! We’re going to watch ALL the action movies from the ’80s and ’90s and then calculate the average!” Quelle: wumo.com

Wir haben uns in dem obigen Beispiel nur ein koninuierliches \(x\) angeschaut. Jetzt kann es aber auch sein, dass deine \(x\)-Werte keine kontinuierlichen Messwerte wie das Gewicht oder die Zeit sind, sondern eben Gruppen. Also du hast verschiedene Düngestufen oder Behandlungsgruppen auf der \(x\)-Achse als Faktoren aufgetragen. Auch dann können wir eine lineare Regression rechnen, eine Linie durch die Punkte legen und anschließend ein Marginal effect model rechnen. Was ist also der Unterschied zwischen einem kontinuierlichen und einem kategoriellen \(x\)-Wert?

53.2.2 …zu simplen und multiplen Modelle

Unterschied zwischen kontinuierlichen und kategoriale \(x\)-Werte

Wir kennen verschiedene Namen für das Gleiche. So nennen wir dann ein kontinuierliches \(x\) dann auch gerne eine stetige Variable oder intervalskaliert. Nichts destotrotz, wir haben ein \(x\) was in kleinen, marignalen Schritten anwachsen kann. Hier kannst du eben an das Gewicht der Flöhe oder aber Zeiteinheiten sowie das Einkommen denken. Wir verändert sich das \(y\), wenn wir die \(x\)-Werte erhöhen?

Auf der anderen Seite haben wir dann kategoriale oder kategorielle \(x\)-Werte. Diese bezeichnen wir dann auch gerne diskret oder aber als Faktoren in R. Wenn wir die Werte von \(x\) ändern, dann springen wir in eine neue Gruppe und es liegt hier eigentlich kein kleiner Schritt vor. Hier haben wir dann eben Düngestufen oder aber Behandlungsgruppen vorliegen. Hier fragen wir uns, wie ändert sich der Wert von \(y\), wenn wir eine Gruppe in \(x\) weiterspringen?

In der folgenden Abbildung von Heiss (2022) siehst du nochmal schön den Unterschied dargestellt. Wir haben beider einer kategorialen Variable einen Schalter. Entweder ist der Schalter an oder eben aus. Im simpelsten Fall haben wir männliche oder eben weibliche Flöhe vorliegen. Das Geschlecht ist somit kategorial. Die Sprungweite oder das Gewicht von Flöhen ist eine kontinuierliche Variable. Wir haben einen Schieberegeler den wir ziemlich fein einstellen können.

Abbildung 53.7— Unterschied zwischen einer kategorialen Variable und einer kontinuierlichen Variable in einem statistischen Modell visualisiert als Schalter und Schieberegler. Übersetzt nach Heiss (2022)

Als wäre das nicht kompliziert genug, schauen wir uns meistens dann nicht nur eine \(x\) Variable in einem Modell an, die wir dann ändern, sondern eben mehrere. Dann kombinieren wir noch gerne kontinuierliche und kategoriale \(x\)-Werte in einem Modell miteinander und erhalten ein Mischboard. Wir können einiges an Schiebereglern und Schaltern in einem Modell betätigen und erhalten entsprechende andere \(y\)-Werte. Hier helfen dann auch Marginal effect models um mehr Erkenntnisse aus einem Modell zu erhalten.

Abbildung 53.8— Kombination verschiedener kategorialer Variablen und kontinuierlichen Variablen in einem statistischen Modell visualisiert als Mischboard. Übersetzt nach Heiss (2022)

Somit kommen wir dann hier mal zu einer Definition, wie wir dann die beiden Arten der möglichen \(x\)-Werte als kontinuierliche und kategoriale Werte sprachlich unterscheiden. Wir immer, je nach wissenschaftlichen Hintergrund können sich dann die Namen ändern und anders sein. Das ist dann eben so in der Statistik.

Marginal effect (deu. marginaler Effekt)

Ein marginaler Effekt beschreibt den statistischen Effekt für kontinuierliche erklärende Variablen; die partielle Ableitung einer Variablen in einem Regressionsmodell; der Effekt eines einzelnen Schiebereglers.

Conditional effect (deu. bedingter Effekt) oder Gruppenkontrast (eng. group contrast)

Ein bedingter Effekt beschreibt den statistischen Effekt für kategoriale erklärende Variablen; der Unterschied in den Mittelwerten, wenn eine Bedingung eingeschaltet ist und wenn sie ausgeschaltet ist; der Effekt eines einzelnen Schalters.

53.2.3 …zur Vorhersage

53.2.4 …zu {marginaleffects} und {emmeans}

Wenn wir Marginal effect models rechnen wollen, dann können wir im Prinzip auf zwei große Pakete zurückgreifen. Einmal das R Paket {marginaleffects} sowie das R Paket {emmeans}. Das R Paket {modelbased} setzt sich im Prinzip auf die beiden Pakete drauf und ist mehr oder minder ein Wrapper mit anderen Funktionsnamen. Das ist eigentlich eine gute Idee und ich zeige dann auch nochmal, wie sich das R Paket {modelbased} verhält. Kommen wir erstmal zu dem hauptsächlichen Unterschied zwischen unseren beiden Elefanten.

Wie unterscheiden sich {emmeans} und {marginaleffects}?

Das R Paket {emmeans} erstellt Durchschnittswerte der Daten und fügt diese Durchschnittswerte dann in Modelle ein. Das R Paket {marginaleffects} fügt alle Werte der Daten in ein Modell ein und erstellt dann Durchschnittswerte aus der Ausgabe des Modells. Am Ende ist es vermutlich dann auch wieder ein nur kleiner Unterschied, der was ausmachen kann. Aber da kommt es dann auf die wissenschaftliche Fragestellung an.

Dabei gibt es noch einen weiteren bedeutenden Unterschied zwischen den beiden Paketen, die sich dann direkt aus der Aggregierung der Daten ableitet. Die Frage ist ja, erst den Mittelwert bilden und dann Modellieren oder umgekehrt. Das R Paket {emmeans} hat als philosophischen Hintergrund experimentelle Daten als Basis. Das R Paket {marginaleffects} hingegen nimmt beobachtete Daten an. Hier möchte ich dann einmal die Vingette des R Pakets {emmeans} zitieren.

“To start off with, we should emphasize that the underpinnings of estimated marginal means – and much of what the {emmeans} package offers – relate more to experimental data than to observational data. In observational data, we sample from some population, and the goal of statistical analysis is to characterize that population in some way. In contrast, with experimental data, the experimenter controls the environment under which test runs are conducted, and in which responses are observed and recorded. Thus with experimentation, the population is an abstract entity consisting of potential outcomes of test runs made under conditions we enforce, rather than a physical entity that we observe without changing it.” — R Paket {emmeans}

Was will uns nun dieser Text sagen und was bedeutet der Unterschied zwischen experimentellen und beobachteten Daten?

• Wir nutzen {emmeans}, wenn wir Gruppenvergleiche aus einem experimentellen, faktoriellen Design rechnen wollen. Solche faktorielle Designs sind in den Agrarwissenschaften sehr häufig.
• Wir nutzen {marginaleffects}, wenn wir beobachtete Daten vorliegen haben. Dies ist sehr häufig bei zeitlichen Verläufen der Fall. Wenn wir also wissen wollen, wie ändert sich den Messwert über die Zeit?

53.3 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

pacman::p_load(tidyverse, gtsummary, marginaleffects, emmeans, scales,
               janitor, ggpmisc, conflicted)
conflicts_prefer(dplyr::mutate)
conflicts_prefer(dplyr::summarize)
conflicts_prefer(dplyr::filter)
conflicts_prefer(ggplot2::annotate)
cb_pal <- c("#000000", "#E69F00", "#56B4E9", 
            "#009E73", "#F0E442", "#F5C710", 
            "#0072B2", "#D55E00", "#CC79A7")
## 
nice_number <- label_number(style_negative = "minus", accuracy = 0.01)
nice_p <- label_pvalue(prefix = c("p < ", "p = ", "p > "))
find_intercept <- function(x1, y1, slope) {
  intercept <- slope * (-x1) + y1
  return(intercept)
}

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

53.4 Daten

Modellierung von Enzymen

enzyme_tbl <- read_excel("data/enzyme_kinetic.xlsx") |> 
  mutate(grp = factor(grp, levels = c("niedrig", "mittel", "hoch")))  

Tabelle 53.5— Daten für die einfaktorielle MANOVA mit der Sprungweite in [cm] und dem Gewicht der Flöhe in [mg] für drei Floharten.

ph	activity	grp	type
-2	-38.05	niedrig	Prokaryot
-2	-43.07	niedrig	Prokaryot
-2	-51.55	niedrig	Prokaryot
…	…	…	…
7.5	57.68	hoch	Eukaryot
7.5	60.68	hoch	Eukaryot
7.5	69.32	hoch	Eukaryot

Tabelle 53.6— Daten für die einfaktorielle MANOVA mit der Sprungweite in [cm] und dem Gewicht der Flöhe in [mg] für drei Floharten.

ph	n	percent
-2.0	3	7.0%
-1.0	5	11.6%
0.0	4	9.3%
1.0	4	9.3%
2.0	4	9.3%
3.0	4	9.3%
4.0	4	9.3%
5.0	3	7.0%
6.0	5	11.6%
7.0	4	9.3%
7.5	3	7.0%

\[ y = x^3 - 8x^2 + 10x + 10 \]

Funktion der Steigung als erste Ableitung nach

\[ y' = 3x^2 - 16x + 10 \]

Abbildung 53.9— Scatterplot der [Zum Vergrößern anklicken]

Abbildung 53.10— (A) Linearer Zusammenhang. (B) gg [Zum Vergrößern anklicken]

Modellierung von Flöhen

flea_model_tbl <- read_excel("data/fleas_model_data.xlsx") |> 
  mutate(feeding = as_factor(feeding),
         stage = as_factor(stage),
         bonitur = as.numeric(bonitur),
         infected = factor(infected, labels = c("healthy", "infected"))) |> 
  select(feeding, stage, jump_length, weight, hatched, count_leg,  bonitur, infected)

Tabelle 53.7— Daten für die einfaktorielle MANOVA mit der Sprungweite in [cm] und dem Gewicht der Flöhe in [mg] für drei Floharten.

feeding	stage	jump_length	weight	count_leg	hatched	bonitur	infected
sugar_water	adult	70.9	16.42	63.4	516.41	4	infected
sugar_water	adult	50.75	12.62	53.96	363.5	1	healthy
sugar_water	adult	62.22	15.57	120.44	303.01	2	healthy
…	…	…	…	…	…	…	…
ketchup	juvenile	41.08	7.18	423.07	429.18	4	infected
ketchup	juvenile	49.68	6.6	550.97	629.58	5	infected
ketchup	juvenile	43.28	4.19	873.28	192.66	5	healthy

Abbildung 53.11— foo. (A) Kombinierte Darstellung in einem Scatterplot (B) Aufgeteilte Darstellung für beide Entwicklungsstadien. [Zum Vergrößern anklicken]

Tabelle 53.8— Deskriptive Statistik des Infektionsstatus (0 = nein / 1 = ja) mit Flohschnupfen aufgeteilt nach den Faktoren und anderen Variablen.

Characteristic	healthy N = 211	infected N = 271
feeding
sugar_water	9 / 21 (43%)	7 / 27 (26%)
blood	3 / 21 (14%)	13 / 27 (48%)
ketchup	9 / 21 (43%)	7 / 27 (26%)
stage
adult	14 / 21 (67%)	10 / 27 (37%)
juvenile	7 / 21 (33%)	17 / 27 (63%)
1 n / N (%)

53.5 Visualisierung von Modellen

In dem folgenden Abschnitten wollen wir immer Modell in unsere Visualisierungen einzeichen. Nehmen wir einmal einen simplen Datensatz,d en wir uns einfach selber bauen und dann wollen wir dort eine Linie durchzeichnen. Dafür nehmen wri einmal zwanzig x-Werte und bauen uns dann die y-Werte nach \(y = 1.5 + 0.75 \cdot x\) zusammen. Dann addieren wir noch einen Fehler aus einer Standardnormalverteilung hinzu. Wenn wir keinen Fehler hinzuaddieren würden, dann lägen die Punkte wie auf einer Perlenschnur aneinandergereit.

set.seed(20250703)
modell_line_tbl <- tibble(x = rnorm(20, 2, 1),
                          y = 1.5 + 0.75 * x + rnorm(length(x), 0, 1))

Jetzt können wir einmal das Modell anpassen und schauen, ob wir die Koeffizienten des Modells wiederfinden. Dann wollen wir natürlich auch sehen, ob unser Modell durch die Punkte läuft. Also erstmal das Modell mit lm() gebaut. Dann schauen wir uns noch die Koeffizienten einmal mit an. Bei nur so wenigen Beobachtungen werden die Koeffizienten aus dem Modell nicht mit den voreingestellten übereinstimmen.

model_fit <- lm(y ~ x, modell_line_tbl)
model_fit

Call:
lm(formula = y ~ x, data = modell_line_tbl)

Coefficients:
(Intercept)            x  
     1.9574       0.5534  

In der folgenden Abbildung siehst du dann einmal den Scatterplot von unseren x-Werten und y-Werten. Wir wollen jetzt die Gerade, die wir im Modell geschätzt haben einmal durch die Punkte legen um zu schauen, ob das Modell auch die Punkte beschreibt. Dabei soll die Gerade durch die Mitte der Punkte laufen und die Punkte sollten auf beiden Seiten der Geraden gleichmäßig verteilt sein.

Abbildung 53.12— Scatterplot der x-Werte und y-Werte. Durch die Punkte soll die Gerade aus dem Modell gelegt werden. [Zum Vergrößern anklicken]

Wir haben jetzt verschiedene Möglichkeiten die Koeffizienten und damit das Modell in den obigen Plot einzuzeichnen. Ich zeige dir hier einmal die häufigsten, die ich dann auch nutze. Erstmal die Anwendung direkt in {ggplot} und dann einmal noch in dem R Paket {ggpmisc}.

…mit {ggplot}

geom_function()

In der Funktion geom_function() müssen wir die Funktion angeben, die wir dann abbilden wollen. Wenn du verstehst, was die Koeffizienten in dem Modell bedeuten, dann kannst du natürlich die mathematische Funktion wie hier entsprechend ergänzen.

ggplot(modell_line_tbl, aes(x, y)) +
  theme_minimal() +
  geom_point() +
  geom_function(fun = \(x) 1.9574 + 0.5534 * x, 
                color = "#CC79A7")

Abbildung 53.13— Scatterplot der x-Werte und y-Werte. Durch die Punkte läuft die Gerade mit den Koeffizienten aus dem Modell. [Zum Vergrößern anklicken]

geom_line(aes(y = predict(model_fit)))

Manchmal ist das Modell zu komplex, dass wir die mathematische Funktion einfach aufschreiben könnten. In dem Fall hilft die Funktion geom_line() die wir dann die vorhergesagten y-Werte mit der Funktion predict() aus dem Modell übergeben. Das funktioniert auch sehr gut.

ggplot(modell_line_tbl, aes(x, y)) +
  theme_minimal() +
  geom_point() +
  geom_line(aes(y = predict(model_fit)), 
                color = "#CC79A7")

Abbildung 53.14— Scatterplot der x-Werte und y-Werte. Durch die Punkte läuft die Gerade mit den vorhergesagten Werten aus dem Modell. [Zum Vergrößern anklicken]

geom_smooth()

Abschließend können wir auch einfach so eine Gerade durch die Punkte legen indem wir die Funktion geom_smooth() als eine Art der Glättung nutzen. Aber hier muss ich sagen, dass uns dann die Geradengleichung fehlt. So mal zum gucken ist das wunderbar. Du kannst über die Option formula auch eine Funktion übergeben. Darüber hinaus erhalten wir dann noch einen Fehlerbalken des Standardfehlers, was in manchen Fällen nützlich ist. Wenn du die Geradengleichung brauchst, dann schaue einmal in dem Paket {ggpmisc} rein.

ggplot(modell_line_tbl, aes(x, y)) +
  theme_minimal() +
  geom_point() +
  geom_smooth(method = "lm", color = "#CC79A7") +
  geom_smooth(method = "lm", formula = y ~ I(x^4), 
              color = "#0072B2")

Abbildung 53.15— Scatterplot der x-Werte und y-Werte. Durch die Punkte läuft die Gerade aus einer Glättung. [Zum Vergrößern anklicken]

…mit {ggpmisc}

Ich möchte hier nich zu sehr in die Tiefe von {ggpmisc} gehen, aber das Paket verbindet im Prinzip die Funktion geom_smooth() mit der Wiedergabe der Informationen zu den Regressionsgleichungen. Du findest bei StackOverflow einmal eine schöne Übersicht in Add regression line equation and R^2 on graph. Wenn du mehr willst, dann schaue dir einmal die Hilfeseite von {ggpmisc} mit Fitted-Model-Based Annotations näher an. Es geht echt eine Menge, von dem ich hier nur einmal den Klassiker zeige. Wir wollen einmal die Regressionsgleichung plus das Bestimmtheitsmaß einzeichnen. Das geht über drei Funktionen zusammen mit der Regressionsgeraden.

ggplot(modell_line_tbl, aes(x, y)) +
  theme_minimal() +
  geom_point() +
  stat_poly_line(color = "#CC79A7") +
  stat_poly_eq(use_label("eq")) +
  stat_poly_eq(label.y = 0.9) 

Abbildung 53.16— Scatterplot der x-Werte und y-Werte. Durch die Punkte läuft die Gerade aus einer Glättung plus die Geradengleichung und das Bestimmtheitsmaß. [Zum Vergrößern anklicken]

53.6 Datenraster

Was ist ein Datenraster (eng. data grid) eigentlich? Wir brauchen die Idee des Datenrasters um überhaupt die Vorhersagen und die Steigungen zu verstehen. Im Prinzip beinhaltet das Datenraster die Information zu welchen Beobachtungen wir eine Vorhersage machen wollen. Wenn wir nochmal kutz zu dem Enyzmebeispiel kommen, zu welchen pH-Werten möchtest du dann eine Steigung oder aber eine Vorhersage der Enzymetätigkeit haben? Beginnen wir mit einem einfachen Beispiel um zu verstehen was die einzelnen Datenraster aussagen wollen.

In Folgenden siehst du einmal einen kleinen Datensatz mit einer numerischen Variable, einer dichotomen Variable sowie einer Variable mit drei Kategorien. So ähnlich haben wir ja auch Datensätze in echt vorliegen.

set.seed(20250709)
grid_tbl <-  tibble(numerisch = rnorm(n = 8),
                    dichotom = rbinom(n = 8, size = 1, prob = 0.5),
                    kategorial = sample(c("niedrig", "mittel", "hoch"), 
                                        size = 8, replace = TRUE))

Beobachtetes Raster (eng. empirical grid)

grid_tbl |> tt()

numerisch	dichotom	kategorial
0.32260788	1	niedrig
0.64136316	1	hoch
-0.81802243	1	mittel
-1.46551819	0	mittel
-0.01230531	1	mittel
-1.61493276	0	mittel
-0.82683177	0	mittel
0.76328096	1	hoch

Abbildung 53.17— foo. [Zum Vergrößern anklicken]

Interessantes Raster (eng. interesting grid)

datagrid(dichotom = c(0, 1), newdata = grid_tbl) |> tt()

numerisch	kategorial	dichotom	rowid
-0.3762948	mittel	0	1
-0.3762948	mittel	1	2

datagrid(numerisch = range, dichotom = mean, kategorial = unique, 
         newdata = grid_tbl) |> tt()

numerisch	dichotom	kategorial	rowid
-1.614933	0.625	niedrig	1
-1.614933	0.625	hoch	2
-1.614933	0.625	mittel	3
0.763281	0.625	niedrig	4
0.763281	0.625	hoch	5
0.763281	0.625	mittel	6

Abbildung 53.18— foo. [Zum Vergrößern anklicken]

Repräsentatives Raster (eng. representative grid)

datagrid(grid_type = "mean_or_mode", newdata = grid_tbl) |> tt()

numerisch	dichotom	kategorial	rowid
-0.3762948	1	mittel	1

Abbildung 53.19— foo. [Zum Vergrößern anklicken]

Balanciertes Raster (eng. balanced grid)

datagrid(grid_type = "balanced", newdata = grid_tbl) |> tt()

numerisch	dichotom	kategorial	rowid
-0.3762948	0	hoch	1
-0.3762948	0	mittel	2
-0.3762948	0	niedrig	3
-0.3762948	1	hoch	4
-0.3762948	1	mittel	5
-0.3762948	1	niedrig	6

Kontrafaktisches Raster (eng. counterfactual grid)

cf_grid <- datagrid(
  dichotom = c(0, 1),
  grid_type = "counterfactual",
  newdata = grid_tbl[1:5,]
)
nrow(cf_grid )

[1] 10

Beobachtete Daten

grid_tbl[1:5,] |> tt()

numerisch	dichotom	kategorial
0.32260788	1	niedrig
0.64136316	1	hoch
-0.81802243	1	mittel
-1.46551819	0	mittel
-0.01230531	1	mittel

Kontrafaktisches Raster

cf_grid |> tt()

rowidcf	numerisch	kategorial	dichotom
1	0.32260788	niedrig	0
2	0.64136316	hoch	0
3	-0.81802243	mittel	0
4	-1.46551819	mittel	0
5	-0.01230531	mittel	0
1	0.32260788	niedrig	1
2	0.64136316	hoch	1
3	-0.81802243	mittel	1
4	-1.46551819	mittel	1
5	-0.01230531	mittel	1

53.7 Steigung (eng. slopes)

“Let us introduce another concept that is likely to get very popular in the near future within the world of regressions. Derivatives.” — {modelbased}

53.7.1 Marginale Effekte

Hier sprechen wir von der Ableitung (eng. derivative)

Modellierung von Flöhen

Abbildung 53.20— foo. (A) Kombinierte Darstellung in einem Scatterplot (B) Aufgeteilte Darstellung für beide Entwicklungsstadien. [Zum Vergrößern anklicken]

Linearer Zusammenhang (\(y \sim x\))

model_ln <- lm(jump_length ~ weight,
               data = flea_model_tbl)
tidy(model_ln)

# A tibble: 2 × 5
  term        estimate std.error statistic  p.value
  <chr>          <dbl>     <dbl>     <dbl>    <dbl>
1 (Intercept)    52.1      5.31       9.81 7.58e-13
2 weight          1.61     0.469      3.43 1.28e- 3

\[ \begin{aligned} \operatorname{E}[y \mid x] &= \beta_0 + \beta_1 x \\[4pt] \frac{\partial \operatorname{E}[y \mid x]}{\partial x} &= \beta_1 \end{aligned} \]

\[ \begin{aligned} \operatorname{E}[\text{Sprungweite} \mid \text{Gewicht}] &= 52.10 + 1.61 \times \text{Gewicht} \\[6pt] \frac{\partial \operatorname{E}[\text{Sprungweite} \mid \text{Gewicht}]}{\partial\ \text{Gewicht}} &= 1.61 \end{aligned} \]

Quadratischer Zusammenhang (\(y \sim x + x^2\))

R Paket {polypoly}

Polynomial Regression - An example

Fitting Polynomial Regression in R

model_sq <- lm(jump_length ~ weight + I(weight^2),
               data = flea_model_tbl)
tidy(model_sq)

# A tibble: 3 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  47.1      10.4        4.53  0.0000437
2 weight        2.62      1.87       1.40  0.169    
3 I(weight^2)  -0.0397    0.0710    -0.559 0.579    

\[ \begin{aligned} \operatorname{E}[y \mid x] &= \beta_0 + \beta_1 x + \beta_2 x^2 \\[4pt] \frac{\partial \operatorname{E}[y \mid x]}{\partial x} &= \beta_1 + 2 \beta_2 x \end{aligned} \]

\[ \begin{aligned} \operatorname{E}[\text{Sprungweite} \mid \text{Gewicht}] &= 47.10 + (2.62 \times \text{Gewicht}) + (−0.04 \times \text{Gewicht}^2) \\[6pt] \frac{\partial \operatorname{E}[\text{Sprungweite} \mid \text{Gewicht}]}{\partial\ \text{Gewicht}} &= 2.62 + (2\times −0.04 \times \text{Gewicht}) \end{aligned} \]

jump_weight_slope <- function(x) 2.62 + (2 *-0.04 * x)
jump_weight_slope(c(5, 10, 15))

[1] 2.22 1.82 1.42

Abbildung 53.21— foo. . [Zum Vergrößern anklicken]

{marginaleffects}

model_sq |> 
  slopes(newdata = datagrid(weight = c(5, 15, 25)))

 weight Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
      5    2.224      1.200 1.854   0.0637 4.0 -0.127   4.57
     15    1.431      0.569 2.514   0.0119 6.4  0.315   2.55
     25    0.638      1.800 0.354   0.7230 0.5 -2.890   4.17

Term: weight
Type: response
Comparison: dY/dX

{emmeans}

model_sq |> 
  emtrends(~ weight, var = "weight",
           at = list(weight = c(5, 15, 25))) |> 
  test()

 weight weight.trend   SE df t.ratio p.value
      5        2.223 1.20 45   1.856  0.0700
     15        1.430 0.57 45   2.510  0.0158
     25        0.637 1.80 45   0.353  0.7254

Abbildung 53.22— foo. . [Zum Vergrößern anklicken]

Modellierung von Enzymen

Nachdem wir das einmal konzeptionellm haben einmal komplexer.

gam()

 ph Estimate Std. Error     z Pr(>|z|)    S  2.5 % 97.5 %
 -1     28.7       3.95  7.27  < 0.001 41.4  20.99  36.47
  2    -11.4       4.12 -2.77  0.00552  7.5 -19.49  -3.35
  6     18.1       4.50  4.02  < 0.001 14.1   9.28  26.92

Term: ph
Type: response
Comparison: dY/dX

Abbildung 53.23— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

poly()

 ph Estimate Std. Error     z Pr(>|z|)     S 2.5 % 97.5 %
 -1     28.3       2.44 11.57   <0.001 100.5  23.5  33.04
  2    -10.3       1.14 -9.08   <0.001  63.0 -12.6  -8.11
  6     21.4       1.63 13.16   <0.001 129.0  18.3  24.64

Term: ph
Type: response
Comparison: dY/dX

Abbildung 53.24— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

loess()

 ph Estimate
 -1    22.37
  2    -4.74
  6    21.85

Term: ph
Type: response
Comparison: dY/dX

Abbildung 53.25— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

Philosophien zur Mittelwertbildung

Average marginal effects (AME)

jump_weight_slope <- function(x) 2.62 + (2 *-0.04 * x)
jump_weight_slope(c(5, 10, 15))

[1] 2.22 1.82 1.42

Abbildung 53.26— foo. Modifiziert nach Heiss (2022)

model_sq <- lm(jump_length ~ weight + I(weight^2),
               data = flea_model_tbl)
tidy(model_sq)

# A tibble: 3 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  47.1      10.4        4.53  0.0000437
2 weight        2.62      1.87       1.40  0.169    
3 I(weight^2)  -0.0397    0.0710    -0.559 0.579    

Abbildung 53.27— foo. [Zum Vergrößern anklicken]

mfx_sq <- slopes(model_sq)
head(mfx_sq)

 Estimate Std. Error    z Pr(>|z|)    S   2.5 % 97.5 %
    1.318      0.702 1.88  0.06037  4.0 -0.0575   2.69
    1.620      0.473 3.43  < 0.001 10.7  0.6935   2.55
    1.386      0.617 2.24  0.02482  5.3  0.1755   2.60
    0.686      1.717 0.40  0.68938  0.5 -2.6786   4.05
    1.316      0.705 1.87  0.06189  4.0 -0.0654   2.70
    1.447      0.554 2.61  0.00899  6.8  0.3614   2.53

Term: weight
Type: response
Comparison: dY/dX

mfx_sq |> 
  group_by(term) |> 
  summarize(avg_slope = mean(estimate))

# A tibble: 1 × 2
  term   avg_slope
  <chr>      <dbl>
1 weight      1.83

avg_slopes(model_sq)

 Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
     1.83      0.618 2.96  0.00305 8.4  0.62   3.04

Term: weight
Type: response
Comparison: dY/dX

Marginal effects at the mean (MEM)

Abbildung 53.28— foo. Modifiziert nach Heiss (2022)

model_sq <- lm(jump_length ~ weight + I(weight^2),
               data = flea_model_tbl)
tidy(model_sq)

# A tibble: 3 × 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)  47.1      10.4        4.53  0.0000437
2 weight        2.62      1.87       1.40  0.169    
3 I(weight^2)  -0.0397    0.0710    -0.559 0.579    

Abbildung 53.29— foo. [Zum Vergrößern anklicken]

avg_jump_weight <- mean(flea_model_tbl$weight)
avg_jump_weight

[1] 9.95375

jump_weight_fitted <- model_sq |> 
  augment(newdata = tibble(weight = c(avg_jump_weight, avg_jump_weight + 0.001)))
jump_weight_fitted

# A tibble: 2 × 2
  weight .fitted
   <dbl>   <dbl>
1   9.95    69.2
2   9.95    69.2

model_sq |> 
  emtrends(~ weight, var = "weight")

 weight weight.trend    SE df lower.CL upper.CL
   9.95         1.83 0.617 45    0.587     3.07

Confidence level used: 0.95 

model_sq |> 
  avg_slopes(newdata = "mean")

 Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
     1.83      0.618 2.96  0.00305 8.4  0.62   3.04

Term: weight
Type: response
Comparison: dY/dX

53.7.2 Gruppierte marginale Effekte

Abbildung 53.30— foo. Modifiziert nach Heiss (2022)

Hier ist die Anordugn wichtig. Erst kommt die gruppierende Variable, dann der Rest.

model_grp_sq <- lm(jump_length ~ stage * weight + I(weight^2),
                   data = flea_model_tbl)
tidy(model_grp_sq)

# A tibble: 5 × 5
  term                 estimate std.error statistic p.value
  <chr>                   <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)           84.4       29.7       2.84  0.00687
2 stagejuvenile        -27.9       29.4      -0.949 0.348  
3 weight                -1.68       3.93     -0.427 0.672  
4 I(weight^2)            0.0752     0.121     0.621 0.538  
5 stagejuvenile:weight   1.60       3.38      0.474 0.638  

Abbildung 53.31— foo. (A) Kombinierte Darstellung in einem Scatterplot (B) Aufgeteilte Darstellung für beide Entwicklungsstadien. [Zum Vergrößern anklicken]

Schritt-für-Schritt

mfx_grp_sq <- model_grp_sq |> 
  slopes(variables = "weight")
head(mfx_grp_sq)

 Estimate Std. Error     z Pr(>|z|)   S  2.5 % 97.5 %
    0.794      0.766 1.037    0.300 1.7 -0.707   2.29
    0.223      1.099 0.203    0.839 0.3 -1.931   2.38
    0.666      0.760 0.876    0.381 1.4 -0.824   2.16
    1.993      2.186 0.912    0.362 1.5 -2.292   6.28
    0.799      0.767 1.041    0.298 1.7 -0.704   2.30
    0.551      0.803 0.685    0.493 1.0 -1.024   2.13

Term: weight
Type: response
Comparison: dY/dX

mfx_grp_sq |> 
  group_by(stage) |> 
  summarize(stage_ame = mean(estimate))

# A tibble: 2 × 2
  stage    stage_ame
  <fct>        <dbl>
1 adult        0.438
2 juvenile     0.807

slopes()

model_grp_sq |> 
  slopes(variables = "weight",
         by = "stage")

    stage Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
 adult       0.438       0.88 0.498    0.618 0.7 -1.29   2.16
 juvenile    0.807       2.24 0.361    0.718 0.5 -3.58   5.19

Term: weight
Type: response
Comparison: dY/dX

Abbildung 53.32— foo. (A) Kombinierte Darstellung in einem Scatterplot (B) Aufgeteilte Darstellung für beide Entwicklungsstadien. [Zum Vergrößern anklicken]

53.7.3 Repräsentative Werte

Abbildung 53.33— foo. Modifiziert nach Heiss (2022)

model_grp_sq <- lm(jump_length ~ feeding * weight + I(weight^2),
                   data = flea_model_tbl)
tidy(model_grp_sq)

# A tibble: 7 × 5
  term                  estimate std.error statistic   p.value
  <chr>                    <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)            48.2      10.8        4.44  0.0000655
2 feedingblood            3.44     12.6        0.273 0.786    
3 feedingketchup        -16.9      13.2       -1.28  0.208    
4 weight                  1.93      1.84       1.05  0.301    
5 I(weight^2)            -0.0382    0.0720    -0.530 0.599    
6 feedingblood:weight     1.05      1.09       0.963 0.341    
7 feedingketchup:weight   2.11      1.18       1.79  0.0804   

Abbildung 53.34— foo. (A) Kombinierte Darstellung in einem Scatterplot (B) Aufgeteilte Darstellung für beide Entwicklungsstadien. [Zum Vergrößern anklicken]

datagrid(model = model_grp_sq,
         weight = c(5, 10, 15),
         feeding = unique)

  weight     feeding rowid
1      5 sugar_water     1
2      5       blood     2
3      5     ketchup     3
4     10 sugar_water     4
5     10       blood     5
6     10     ketchup     6
7     15 sugar_water     7
8     15       blood     8
9     15     ketchup     9

Ohne Gruppenvergleich

model_grp_sq |> 
  slopes(variables = "weight",
         newdata = datagrid(weight = c(5, 10, 15),
                            feeding = unique))

 weight     feeding Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
      5 sugar_water    1.552      1.195 1.30  0.19423 2.4 -0.791   3.89
      5 blood          2.597      1.563 1.66  0.09658 3.4 -0.466   5.66
      5 ketchup        3.660      1.416 2.59  0.00973 6.7  0.885   6.44
     10 sugar_water    1.170      0.701 1.67  0.09522 3.4 -0.204   2.54
     10 blood          2.215      1.021 2.17  0.02999 5.1  0.215   4.22
     10 ketchup        3.279      1.031 3.18  0.00147 9.4  1.259   5.30
     15 sugar_water    0.788      0.770 1.02  0.30582 1.7 -0.720   2.30
     15 blood          1.834      0.824 2.22  0.02613 5.3  0.218   3.45
     15 ketchup        2.897      1.076 2.69  0.00708 7.1  0.789   5.01

Term: weight
Type: response
Comparison: dY/dX

Mit Gruppenvergleich

comp_tbl <- as.data.frame(t(combn(levels(flea_model_tbl$feeding), 2))) |> 
    unite("z", sep = " - ")

model_grp_sq |> 
  slopes(variables = "weight",
         newdata = datagrid(weight = c(5),
                            feeding = unique))

 weight     feeding Estimate Std. Error    z Pr(>|z|)   S  2.5 % 97.5 %
      5 sugar_water     1.55       1.20 1.30  0.19423 2.4 -0.791   3.89
      5 blood           2.60       1.56 1.66  0.09658 3.4 -0.466   5.66
      5 ketchup         3.66       1.42 2.59  0.00973 6.7  0.885   6.44

Term: weight
Type: response
Comparison: dY/dX

model_grp_sq |> 
  slopes(variables = "weight",
         newdata = datagrid(weight = c(5),
                            feeding = unique),
         hypothesis = "pairwise")  |> 
  mutate(term = pluck(comp_tbl, "z"))

model_grp_sq |> 
  emtrends(~ weight + feeding, 
           var = "weight",
           at = list(weight = c(5, 10, 15),
                     feeding = c("sugar_water", "ketchup", "blood")),
           regrid = "response") 

 weight feeding     weight.trend    SE df lower.CL upper.CL
      5 sugar_water        1.551 1.190 41   -0.860     3.96
     10 sugar_water        1.169 0.701 41   -0.246     2.58
     15 sugar_water        0.787 0.771 41   -0.769     2.34
      5 ketchup            3.659 1.410 41    0.803     6.52
     10 ketchup            3.278 1.030 41    1.197     5.36
     15 ketchup            2.896 1.080 41    0.722     5.07
      5 blood              2.596 1.560 41   -0.557     5.75
     10 blood              2.215 1.020 41    0.154     4.27
     15 blood              1.833 0.825 41    0.167     3.50

Confidence level used: 0.95 

53.8 Vorhersagen (eng. predictions)

loess_fit <- loess(activity ~ ph, data = enzyme_tbl)

predictions(loess_fit)

Warning: Unable to extract a variance-covariance matrix from this model.

 Estimate
   -37.22
   -37.22
   -37.22
    -9.33
    -9.33
--- 33 rows omitted. See ?print.marginaleffects ---
    32.90
    32.90
    52.13
    52.13
    52.13
Type: response

predictions(loess_fit, by = "ph")

Warning: Unable to extract a variance-covariance matrix from this model.

   ph Estimate
 -2.0   -37.22
 -1.0    -9.33
  0.0     7.60
  1.0    13.93
  2.0     8.54
  3.0    -1.63
  4.0   -10.10
  5.0    -9.90
  6.0     4.93
  7.0    32.90
  7.5    52.13

Type: response

Kontinuierliche \(x\)-Werte

gam()

gam_fit <- gam(activity ~ s(ph), data = enzyme_tbl)

gam_pred <- predictions(gam_fit, newdata = datagrid(ph = c(-1, 2, 6)))
gam_pred

 ph Estimate Std. Error      z Pr(>|z|)   S  2.5 % 97.5 %
 -1    -5.73       3.57 -1.605   0.1085 3.2 -12.73   1.27
  2     8.94       3.73  2.398   0.0165 5.9   1.63  16.24
  6    -3.21       3.44 -0.935   0.3500 1.5  -9.95   3.53

Type: response

poly()

poly_fit <- lm(activity ~ poly(ph, 3), data = enzyme_tbl)

poly_pred <- predictions(poly_fit, newdata = datagrid(ph = c(-1, 2, 6)))
poly_pred

 ph Estimate Std. Error      z Pr(>|z|)    S 2.5 % 97.5 %
 -1   -4.545       2.48 -1.830   0.0673  3.9 -9.41  0.324
  2    8.951       2.23  4.013   <0.001 14.0  4.58 13.323
  6   -0.545       2.36 -0.231   0.8176  0.3 -5.18  4.089

Type: response

loess()

loess_fit <- loess(activity ~ ph, data = enzyme_tbl)

loess_pred <- predictions(loess_fit, newdata = datagrid(ph = c(-1, 2, 6)))

Warning: Unable to extract a variance-covariance matrix from this model.

loess_pred

 ph Estimate
 -1    -9.33
  2     8.54
  6     4.93

Type: response

Abbildung 53.35— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

Kategoriale \(x\)-Werte

Achtung, bitte beachten!

Ich persönlich finde die Implementierung des multiplen Testens in {emmeans} um Längen besser gelöst. Den Rest von {marginaleffects} dann eher nicht so. Daher würde ich dir hier davon abraten, deine Gruppenvergleiche mit predictions() zu rechnen. Es ist gut zu verstehen was die Funktion macht, aber {emmeans} hat den klaren Vorteil, dass wir das Compact letter disply berechnen können. Darüber hinaus finde ich die zweifaktorielle Analyse besser gelöst durch die beiden Zeichen Stern * und Strich | besser gelöst. In {marginaleffects} haben wir dann die Problematik mit den Hypothesen und Gruppennamen innerhalb der Ausgabe der Funktion hypotheses(). Hier nochmal die Alternativen in {emmeans}, die ich hier aber dann nicht ausführen lasse. Bitte schaue dann nochmal in das Kapitel zu den Post-hoc Tests vorbei.

enzyme_1fac_fit |> 
  emmeans(~ grp, vcov = sandwich::vcovHAC) |> 
  contrast(method = "pairwise", adjust = "bonferroni")

enzyme_2fac_fit |> 
  emmeans(~ grp | type, vcov = sandwich::vcovHAC) |> 
  contrast(method = "pairwise", adjust = "bonferroni")

“Manchmal schaue ich mir sehr lange Funktionen oder Pakete wie hier {marginaleffects} an und muss dann feststellen, dass für eien spezielle Orichideenanwendung dann doch ein anderes Paket besser ist. Ist nicht schlimm, dann musst du dir nicht die Arbeit machen. Für das faktorielle Experiment würde ich dann immer {emmeans} nehmen.” — Jochen Kruppa-Scheetz, meiner bescheidener Meinung nach.

Die Funktion hypotheses() ist mächtig.

Marginal mean nochmal erklären

Vorhersage 1-faktoriell

enzyme_1fac_fit <-  lm(activity ~ grp, data = enzyme_tbl)

predictions(enzyme_1fac_fit, by = c("grp"), vcov = "HAC")

     grp Estimate Std. Error      z Pr(>|z|)    S  2.5 % 97.5 %
 niedrig    -2.50       3.99 -0.626    0.531  0.9 -10.32  5.327
 mittel     -1.35       1.06 -1.276    0.202  2.3  -3.42  0.722
 hoch       17.11       4.91  3.483   <0.001 11.0   7.48 26.745

Type: response

predictions(enzyme_1fac_fit, by = c("grp"), vcov = "HAC",
            hypothesis = ~pairwise)  |> 
  hypotheses(multcomp = "bonferroni")

 Term           Hypothesis Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
   b1 (mittel) - (niedrig)     1.15       4.67 0.246  1.00000 -0.0 -9.68   12.0
   b2 (hoch) - (niedrig)      19.61       7.07 2.774  0.01659  5.9  3.21   36.0
   b3 (hoch) - (mittel)       18.46       5.59 3.304  0.00286  8.4  5.50   31.4

predictions(enzyme_1fac_fit, by = c("grp"), vcov = "HAC",
            hypothesis = ~pairwise) 

           Hypothesis Estimate Std. Error     z Pr(>|z|)    S 2.5 % 97.5 %
 (mittel) - (niedrig)     1.15       4.67 0.246  0.80542  0.3 -8.00   10.3
 (hoch) - (niedrig)      19.61       7.07 2.774  0.00553  7.5  5.76   33.5
 (hoch) - (mittel)       18.46       5.59 3.304  < 0.001 10.0  7.51   29.4

Type: response

Abbildung 53.36— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

Vorhersage 2-faktoriell

enzyme_2fac_fit <- lm(activity ~ grp * type, data = enzyme_tbl)

predictions(enzyme_2fac_fit, by = c("grp", "type"), vcov = "HAC",
            hypothesis = ~pairwise | type, 
            newdata = "balanced") 

      type           Hypothesis Estimate Std. Error      z Pr(>|z|)    S  2.5 %
 Eukaryot  (mittel) - (niedrig)   -23.43      0.571 -41.00  < 0.001  Inf -24.55
 Eukaryot  (hoch) - (niedrig)      14.03     13.570   1.03  0.30109  1.7 -12.56
 Eukaryot  (hoch) - (mittel)       37.46     13.408   2.79  0.00521  7.6  11.18
 Prokaryot (mittel) - (niedrig)    21.79      3.460   6.30  < 0.001 31.6  15.01
 Prokaryot (hoch) - (niedrig)      14.07      4.257   3.30  < 0.001 10.0   5.72
 Prokaryot (hoch) - (mittel)       -7.72      1.647  -4.69  < 0.001 18.5 -10.95
 97.5 %
  -22.3
   40.6
   63.7
   28.6
   22.4
   -4.5

Type: response

predictions(enzyme_2fac_fit, by = c("grp", "type"), vcov = "HAC",
            newdata = "balanced") |> 
  hypotheses(~pairwise | type, multcomp = "bonferroni")

      type  Hypothesis Estimate Std. Error      z Pr(>|z|)    S  2.5 % 97.5 %
 Eukaryot  (b3) - (b1)   -23.43      0.571 -41.00  < 0.001  Inf -24.84 -22.01
 Eukaryot  (b5) - (b1)    14.03     13.570   1.03  1.00000 -0.0 -19.59  47.66
 Eukaryot  (b5) - (b3)    37.46     13.408   2.79  0.03124  5.0   4.24  70.69
 Prokaryot (b4) - (b2)    21.79      3.460   6.30  < 0.001 29.0  13.22  30.37
 Prokaryot (b6) - (b2)    14.07      4.257   3.30  0.00571  7.5   3.52  24.62
 Prokaryot (b6) - (b4)    -7.72      1.647  -4.69  < 0.001 15.9 -11.81  -3.64

plot_predictions(enzyme_2fac_fit, by = c("grp", "type")) +
  theme_minimal()

Abbildung 53.37— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

Vergleich {marginaleffects} und {emmeans}

feeding_fit <-  lm(jump_length ~ feeding * stage, data = flea_model_tbl)

{marginaleffects}

predictions(feeding_fit, by = c("stage", "feeding"), vcov = "HAC")

    stage     feeding Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
 adult    sugar_water     68.9       4.62 14.9   <0.001 164.7  59.8   77.9
 adult    blood           82.1       5.73 14.3   <0.001 152.1  70.8   93.3
 adult    ketchup         81.1       4.77 17.0   <0.001 212.8  71.8   90.5
 juvenile sugar_water     57.3       3.20 17.9   <0.001 235.2  51.0   63.5
 juvenile blood           68.8       6.25 11.0   <0.001  91.3  56.6   81.0
 juvenile ketchup         50.3       2.37 21.2   <0.001 328.6  45.6   54.9

Type: response

predictions(feeding_fit, by = c("stage", "feeding"), vcov = "HAC",
            hypothesis = ~pairwise | stage)

    stage                Hypothesis Estimate Std. Error      z Pr(>|z|)   S
 adult    (blood) - (sugar_water)     13.185       7.36  1.791   0.0733 3.8
 adult    (ketchup) - (sugar_water)   12.246       6.64  1.843   0.0653 3.9
 adult    (ketchup) - (blood)         -0.939       7.47 -0.126   0.9001 0.2
 juvenile (blood) - (sugar_water)     11.532       7.02  1.644   0.1002 3.3
 juvenile (ketchup) - (sugar_water)   -6.981       3.99 -1.752   0.0798 3.6
 juvenile (ketchup) - (blood)        -18.514       6.68 -2.770   0.0056 7.5
   2.5 % 97.5 %
  -1.245  27.61
  -0.776  25.27
 -15.589  13.71
  -2.217  25.28
 -14.793   0.83
 -31.613  -5.41

Type: response

predictions(feeding_fit, by = c("stage", "feeding"), vcov = "HAC") |> 
  hypotheses(multcomp = "bonferroni")

 Term Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
   b1     68.9       4.62 14.9   <0.001 Inf  56.7   81.0
   b2     82.1       5.73 14.3   <0.001 Inf  67.0   97.2
   b3     81.1       4.77 17.0   <0.001 Inf  68.6   93.7
   b4     57.3       3.20 17.9   <0.001 Inf  48.8   65.7
   b5     68.8       6.25 11.0   <0.001 Inf  52.4   85.2
   b6     50.3       2.37 21.2   <0.001 Inf  44.0   56.5

Abbildung 53.38— (A) Linearer Zusammenhang. (B) ggg [Zum Vergrößern anklicken]

{emmeans}

feeding_fit |> 
  emmeans(~ feeding | stage, vcov = sandwich::vcovHAC)

stage = adult:
 feeding     emmean   SE df lower.CL upper.CL
 sugar_water   68.9 4.62 42     59.6     78.2
 blood         82.1 5.73 42     70.5     93.6
 ketchup       81.1 4.77 42     71.5     90.8

stage = juvenile:
 feeding     emmean   SE df lower.CL upper.CL
 sugar_water   57.3 3.20 42     50.8     63.7
 blood         68.8 6.25 42     56.2     81.4
 ketchup       50.3 2.37 42     45.5     55.1

Confidence level used: 0.95 

feeding_fit |> 
  emmeans(~ feeding | stage, vcov = sandwich::vcovHAC) |> 
  contrast(method = "pairwise", adjust = "bonferroni")

stage = adult:
 contrast              estimate   SE df t.ratio p.value
 sugar_water - blood    -13.185 7.36 42  -1.791  0.2416
 sugar_water - ketchup  -12.246 6.64 42  -1.843  0.2171
 blood - ketchup          0.939 7.47 42   0.126  1.0000

stage = juvenile:
 contrast              estimate   SE df t.ratio p.value
 sugar_water - blood    -11.533 7.02 42  -1.644  0.3230
 sugar_water - ketchup    6.981 3.99 42   1.752  0.2614
 blood - ketchup         18.514 6.68 42   2.770  0.0249

P value adjustment: bonferroni method for 3 tests 

53.9 Kontrafaktische Vergleiche (eng. counterfactual)

53.10 Hypothesen- und Äquivalenztests

Referenzen

Arel-Bundock, V., Greifer, N., & Heiss, A. (2024). How to Interpret Statistical Models Using marginaleffects for R and Python. Journal of Statistical Software, 111(9), 1–32. https://doi.org/10.18637/jss.v111.i09

Heiss, A. (2022, Mai 20). Marginalia: A Guide to Figuring Out What the Heck Marginal Effects, Marginal Slopes, Average Marginal Effects, Marginal Effects at the Mean, and All These Other Marginal Things Are. https://doi.org/10.59350/40xaj-4e562