45  Multiple lineare Regression

Letzte Änderung am 01. August 2025 um 22:00:45

“All models are approximations. Essentially, all models are wrong, but some are useful. However, the approximate nature of the model must always be borne in mind.” — George E. P. Box

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.

In diesem Kapitel geht es um die multiple lineare Regression. Der große Unterschied zu der simplen linearen Regression ist, dass wir in der multiplen Regression mehr als eine Einflussvariable vorliegen haben. Damit wird unser Modell um einiges komplexer. Hier gibt es jetzt gewisse Überschneidungen mit den Marginal effect models, die im Prinzip die Erweiterung der multiplen Regression sind. Wir können mit den Marginal effect models eben dann doch etwas anders lineare und nicht linare Regressionen interpretieren. Damit kannst du in dem Kapitel dann mehr über eine alternative Art der Interpretation von multiplen Modellen lernen. Ebenso gehe nicht hier so stark auf das Modellieren in R ein. Hier kannst du dann in dem entsprechenden Kapitel auch nochmal mehr lesen. Hier zeige ich dann die Basisanwendungen. Am Ende habe ich mich auch entschieden hier nur einen normalverteilten Messwert zu betrachten. Sonst wird es einfach viel zu komplex, wenn wir hier noch andere Verteilungen berücksichtigen. Dazu dann mehr in den entsprechenden Kapiteln zur statistischen Modellierung.

“In den folgenden Kapiteln gibt es immer mal wieder Redundanzen. Teilweise sind diese Redundanzen gewollt, denn eine gute Widerholung hilft ja immer. Meistens habe ich aber den Überblick verloren und fand, dass das Thema dann an der Stelle doch nochmal gut reinpasste. Ist so passiert.” — Jochen Kruppa-Scheetz, meiner bescheidener Meinung nach.

45.1 Allgemeiner Hintergrund

“This idea of ‘holding everything constant’ though can be tricky to wrap your head around.” — Andrew Heiss

Wenn wir über die multiple lineare Regression sprechen, dann sprechen wir immer über mehr als eine Einflussvariable. Dadurch können wir den Fall haben, dass wir nicht nur eine Kovariate sondern mehrere Kovariaten in einem Modell vorliegen haben. Zum Beispiel wollen wir wissen, ob die Sprungweite von dem Gewicht der Flöhe oder der Anzahl an Beinhaaren abhängt. Wir haben hier also mehrere kontinuierliche Einflussvariablen vorliegen. Ebenso können wir den Fall haben, dass uns nur kategoriale Einflussvariablen interessieren. Hat die Ernährungsform oder aber der Entwicklungsstand einen Einfluss auf die Sprungweite der Flöhe. Klassischerweise wären wir hier zwar auf dem ANOVA Pfad eines faktoriellen Experiments, aber wir können natürlich die Koeffzienten eines mehfaktoriellen Modells in der Form einer multiplen Regression interpretieren. Wenn wir dann keinen normalverteilten Messwert haben, dann bietet sich natürlich dennoch die multiple Regression als eine Analyseform an. Abschließend können wir natürlich auch kombinierte Modelle vorliegen haben. Wir haben also in unseren multiple Modell nicht nur Kovariaten oder Faktoren sondern eben eine Mischung aus beiden Formen. Das führt dann zu einer noch komplexeren Interpretation. Alle Fälle wollen wir uns hier in diesem Kapitel einmal an einem normalverteilten Messwert anschauen.

Sprachlicher Hintergrund

“In statistics courses taught by statisticians we don’t use independent variable because we use independent on to mean stochastic independence. Instead we say predictor or covariate (either). And, similarly, we don’t say”dependent variable” either. We say response.” — User berf auf r/AskStatistics

Wenn wir uns mit dem statistischen Modellieren beschäftigen wollen, dann brauchen wir auch Worte um über das Thema reden zu können. Statistik wird in vielen Bereichen der Wissenschaft verwendet und in jedem Bereich nennen wir dann auch Dinge anders, die eigentlich gleich sind. Daher werde ich mir es hier herausnehmen und auch die Dinge so benennen, wie ich sie für didaktisch sinnvoll finde. Wir wollen hier was verstehen und lernen, somit brauchen wir auch eine klare Sprache.

“Jeder nennt in der Statistik sein Y und X wie er möchte. Da ich hier nicht nur von Y und X schreiben will, führe ich eben die Worte ein, die ich nutzen will. Damit sind die Worte dann auch richtig, da der Kontext definiert ist. Andere mögen es dann anders machen. Ich mache es eben dann so. Danke.” — Jochen Kruppa-Scheetz, meiner bescheidener Meinung nach.

In dem folgenden Kasten erkläre ich nochmal den Gebrauch meiner Begriffe im statistischen Testen. Es ist wichtig, dass wir hier uns klar verstehen. Zum einen ist es angenehmer auch mal ein Wort für ein Symbol zu schreiben. Auf der anderen Seite möchte ich aber auch, dass du dann das Wort richtig einem Konzept im statistischen Modellieren zuordnen kannst. Deshalb einmal hier meine persönliche und didaktische Zusammenstellung meiner Wort im statistischen Modellieren. Du kannst dann immer zu dem Kasten zurückgehen, wenn wir mal ein Wort nicht mehr ganz klar ist. Die fetten Begriffe sind die üblichen in den folgenden Kapiteln. Die anderen Worte werden immer mal wieder in der Literatur genutzt.

Worte und Bedeutungen im statistischen Modellieren

“Die Grenzen meiner Sprache bedeuten die Grenzen meiner Welt.” — Ludwig Wittgenstein

Hier kommt einmal die Tabelle mit den wichtigsten Begriffen im statistischen Modellieren und wie ich die Worte benutzen werde. Damit wir uns verstehen und du was lernen kannst. In anderen Büchern und Quellen findest du teilweise die Worte in einem anderen Sinnzusammenhang. Das ist gut so dort. Bei mir ist es anders.

Tabelle 45.1— Zusammenfassende Tabelle des Sprachgebrauchs in meinem statistischen Modellieren. Die fett hervorgehobenen Bezeichnungen werden hier genutzt. (LHS = left hand side, RHS = right hand side)

	Symbol	Deutsch	Englisch
LHS	\(Y\) / \(y\)	Messwert / Endpunkt / Outcome	response / outcome / endpoint
RHS	\(X\) / \(x\)	Einflussvariable / Erklärende Variable / Fester Effekt	risk factor / explanatory / fixed effect
RHS	\(Z\) / \(z\)	Zufälliger Effekt	random effect
\(X\) ist kontinuierlich	\(c_1\)	Kovariate 1	covariate 1
\(X\) ist kategorial	\(f_A\)	Faktor \(A\) mit Level \(A.1\) bis \(A.j\)	factor \(A\) with levels \(A.1\) to \(A.j\)

Wie multiple Einflussvariablen interpretieren?

Wie schon im Zitat am Anfang erwähnt, heißt die Regel, dass wir jeden Koeffizienten unseres Modell nur in soweit interpretieren können, wenn sich die anderen Einflussvariablen nicht ändern. Das klingt jetzt kryptisch und macht eigentlich nur Sinn, wenn man schon alles verstanden hat. Deshalb hier einmal die Analogie mit Schaltern und Schieberegelern. Wir können eine kategoriale Variable oder Faktor durch die Level durchschalten. Daher können wir uns den Effekt von weiblichen oder männlichen Flöhen anschauen. Es gibt aber nur diese beiden Möglichkeiten. Entweder steht der Schalter auf männlichen Floh oder aber auf weiblichen Floh. Etwas anderes ist es wenn wir uns eine Kovariate vorliegen haben. Eine kontinuierliche Einflussvariable können wir beliebig einstellen. Wir können das Gewicht eines Flohs eben hin und herschieben um einen Wert zu haben.

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

In einer multiplen Regression kombinieren wir nun alles miteinander. Wir haben eben eine Art Mischpult vorliegen. Auf unserem Mischpult finden wir eben Schalter und Schieberegler. Wir haben dann auch in unserem Modell Faktoren und Kovariaten vorliegen. Dann können wir auch Schalter auf verschiedene Positionen setzen und die Regler verschieden einstellen.

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

Worher kommt nun die Interpretation, dass wir eine Einflussvaribale nur interpretieren können, wenn wir den Rest der Einflussvariablen konstant halten? Wenn wir an dem Schieberegler des Gewichts der Flöhe drehen, dann müssen die anderen Schalter und Regler still stehen. Wir interpretieren dann eben das sich änderende Gewicht, für weibliche Flöhe unter der Ernährung mit Ketchup. Wir können nicht global das gesmate Modell interpretieren. Einen Ausweg bieten die Marginal effect models, aber diese werden wir in diesem Kapitel nur am Rande behandeln. Dafür ist das Thema zu komplex. Daher behalte ich Kopf, wenn wir später einzelne Einflussvariablen betrachten, dann tuen wir das immer im Kontext fixierter anderer Einflussvariablen.

Wie verschiedene Modelle darstellen?

Häufig stellt sich dann die Frage, wie so multiple Regressionsmodelle dargestellt werden können. Daher hier gleich einmal vorweg ein Beispiel für eine Darstellung einer multiplen linearen Regression in einer wissenschaftlichen Arbeit. In der Arbeit Energy expenditure and obesity across the economic spectrum von McGrosky et al. (2025) wird die multiple lineare Regression genutzt um den Zusammenhang von Energieaufnahme und dem Energieverbauch in veschiedenen Entwicklungsstufen von Ländern zu modellieren. Dabei kommt eigentlich ein recht spannendes Ergebnis heraus, wenn ich das einmal sehr kanpp zitieren darf. Es kommt eben mehr auf die Energieaufnahme als auf den Verbrauch an. Oder anders herum, wenn du abnehmen willst achte auf dein Essen und nicht so sehr auf deine Bewegung.

“Our analyses suggest that increased energy intake has been roughly 10 times more important than declining total energy expenditure (TEE) in driving the modern obesity crisis.” — McGrosky et al. (2025)

In der folgenden Tabelle siehst du dann eimal das multiple lineare Modell mit dem prozentualen Körperfettanteil als Messwert. Es wurden dann acht verschiedene multiple Modelle gerechnet und in jedem Modell wurden andere Einflussvariablen aufgenommen. In der vorletzten Zeile findest du dann noch das Bestimmtheitsmaß \(R^2\), was dir angibt wie gut das Modell funktioniert hat. Hier wollen wir eigentlich einen Wert gegen Eins sehen, aber in Beobachtungsdaten wie hier, ist eine höhere Zahl erstmal eine bessere Zahl.

Abbildung 45.3— Lineare Modelle für den Einfluss des Anteils ultraverarbeiteter Lebensmittel in der Ernährung (Anteil an der gesamten Kalorienaufnahme, %UPF) und des Fleischkonsums (kg pro Kopf) auf den Körperfettanteil. Signifikante Effekte (*p<0,05, **p<0,01, ***p<0,001) in Fettdruck, negative Effekte in Blau. Quelle: McGrosky et al. (2025)

Was sagt uns jetzt die Tabelle aus? Zuerst einmal rechnen wir hier acht Modelle in den Spalten der Tabelle, die wir dann untereinander vergleichen. Wir nehmen verschiedene Einflussvariablen mit in das Modell, die alle in den Zeilen stehen. Teilweise ist das bei kategorialen Variablen wie dem Geschlecht die Gruppe angegeben. Wie interpretieren wir nun diese Tabelle? Hier einmal ein Auszug zu dem Einfluss des Fleischkonsum (eng. per capita meat consumption) in den letzten drei Spalten auf den Messwert Körperfetrprozent.

“[…], per capita meat consumption, another dietary change associated with development, was not significantly associated with fat percentage when added to models including percent ultraprocessed foods (UPF) [i.e., industrial formulations of five or more ingredients (29)]” — McGrosky et al. (2025)

Was lesen wir hier? Wenn wir den Fleischkonsum zusammen mit den hochverarbeiteten Lebensmittel mit ins Modell nehmen, dann wird der Fleichkonsum nicht signifikant. Dies ist eine Art eine multiple lineare Regression zu rechnen und die Einflussvariablen untereinander zu vergleichen. Mehr dazu dann in den folgenden Abschnitten und den entsprechenden Kapiteln zur statistischen Modellierung.

Defintion von hochverarbeiteten Lebensmittel (eng. ultraprocessed foods)

Wir können in der Abrbeit von Gibney (2019) nocheinmal nachlesen, was die aktuelle Definition eines hochverarbeiteten Lebensmittel (eng. ultraprocessed foods) ist. Das ist ja immer wieder ein Streitpunkt, der aktuell ja auch noch offen ist. Wir halten uns einmal an die folgende Definition.

“Industrial formulations typically with 5 or more and usually many ingredients. Besides salt, sugar, oils, and fats, ingredients of ultra-processed foods include food substances not commonly used in culinary preparations, such as hydrolyzed protein, modified starches, and hydrogenated or interesterified oils, and additives whose purpose is to imitate sensorial qualities of unprocessed or minimally processed foods and their culinary preparations or to disguise undesirable qualities of the final product, such as colorants, flavorings, nonsugar sweeteners, emulsifiers, humectants, sequestrants, and firming, bulking, de-foaming, anticaking, and glazing agents.” — Gibney (2019)

Jetzt kommen wir nochmal zu dem theoretischen Hintergrund. Wir wollen hier nochmal gleich verstehen, wie die multiplen Modelle aufgebaut sind und wie die Modelle in R dargestellt werden. Wenn du dich mit dem allgemeinen Hintergrund zufrieden gibst, dann kannst du auch gerne direkt in die Beispiele und deren Auswertung springen. Oder eben dann nochmal später zurückkommen, wenn der Bedarf besteht.

45.2 Theoretischer Hintergrund

In dem theoretischen Hintergrund der multiplen Modelle bleibe ich etwas oberflächlich und gehe nicht auf die mathematischen Hintergründe tiefer ein. Das lohnt hier auch nicht weiter. In den Anwendungen habe ich dann immer noch die theoretische Berechnung in R ergänzt, die dir dann nochmal zeigt, wie in R intern die Modelle gerechnet werden. Ganz allgemein kann ich das Buch von Kéry (2010) wärmstens empfehlen. Dort habe ich viel über die Modelle und die Berechnungen gelernt. In dem Zusammenhang kann ich auch noch Dormann (2013) empfehlen, wenn es um die Maximum Likelihood Methode geht, die dort gut beschrieben ist. Ich gebe hier dann die Zusammenfassung wieder.

Beginnen wir nochmal mit dem klassischen statistischen Modell in der folgenden Abbildung. Wir haben Daten vorliegen, die wir in ein Modell und einen Fehler zerlegen wollen. Das liest sich sehr allgemein, in der Praxis haben wir auf der linken Seite den beobachteten Messwert Y stehen und auf der linken Seite das Modell aus dem der erklärte Anteil am Y rauskommt. Die Differenz zu den beobachteten Messwerten ist dann der Fehler. Somit ist unser Modell die Kombination aus dem Messwert und unseren Einflussvariablen. Was die Einflussvariablen an dem Messwert nicht erklären können wandert in den Fehler. Somit wollen wir einen geringen Fehler in unserer Modellierung.

Abbildung 45.4— Der Zusammenhang von Daten, dem statistischen Modell und dem Fehler. Das Modell versucht durch ein statistisches Modell \(f(x)\) den Zusammenhang zwischen den Einflussvariablen dem Messwert zu erklären. Den Anteil des unerklärten Messwert geht in den Fehlerterm. [Zum Vergrößern anklicken]

Erinnern wir uns nochmal an das simple und multiple lineares Modell. Die beiden Modelle entscheiden sich durch die Anzahl an Einflussvariablen. In der folgenden Abbildung siehst du nochmal das simple lineare Modell. Wir haben hier eine Einflussvariable X. Dann kommen noch die Koeffizienten \(\beta_0\) für den Intercept und der Koeffizient \(\beta_1\) für die Steigung der eine Einflussvariable hinzu. Alles was wir nicht durch die Grae erklären können, wandert dann in die Residuen oder Fehler \(\epsilon\). Eine klassische Gradengleichung mit einer Einflussvariable.

Abbildung 45.5— Formelschreibweise der simplen linearen Regression beinhaltend die Koeffizienten \(\beta_0\) für den y-Achsenabschnitt sowie \(\beta_1\) für die Steigung der Graden für eine Einflussvariable \(x_1\). Die Residuen werden durch \(\epsilon\) abgebildet. [Zum Vergrößern anklicken]

In diesem Kapitel erweitern wir das simple lineare Modell um mehr Einflussvariablen. Dann erhalten wir ein multiples lineares Modell wie in der folgenden Abbildung. Meistens haben wir dann mehr als zwei Einflussvariablen sondern eben \(p\) Stück. Jede dieser Einflussvariablen hat dann eine eigene Stigung. Wir haben dann aber immer noch nur einen Intercept vorliegen. Auch haben wir dann am Ende nur einen Wert für die Residuen oder Fehler. Die Einflussvariablen können eine beliebige Kombination an Kovariten und Faktoren sein.

Abbildung 45.6— Formelschreibweise der multiplen linearen Regression beinhaltend die Koeffizienten \(\beta_0\) für den y-Achsenabschnitt sowie \(\beta_1\) bis \(\beta_p\) für die partielle Steigung der Graden für jede Einflussvariable \(x_1\) bis \(x_p\). Die Residuen werden durch \(\epsilon\) abgebildet. [Zum Vergrößern anklicken]

Oder konkreter können wir sagen, dass jump_length von animal und weight abhängt und wir diesen Zusammenhang modellieren wollen.

Abbildung 45.7— Statistische Modellschreibweise in R mit dem Messwert auf der linken Seite und den Einflussvariablen auf der rechten Seite der Tilde. Die Platzhalter \(Y\) und \(X\) werden durch die Spaltennamen im Datensatz ersetzt. [Zum Vergrößern anklicken]

Abbildung 45.8— Erweiterte Darstellung der statistischen Modellierung. Die Messwerte \(Y\) folgen einer Verteilung. Die Einflussvariablen liegen kontinuierlich als Kovariaten oder aber kategorial als Faktoren vor. [Zum Vergrößern anklicken]

Wenn du jetzt denkst ‘Hä? Was soll denn mehrdimensional bedeuten?’, dann besuche doch einmal die fantastische Seite Explained Visually | Ordinary Least Squares Regression um selber zu erfahren, was eine multiple lineare Regression macht. Auf der Seite findest du interaktive Abbildungen, die dir das Konzept der linearen Regression sehr anschaulich nachvollziehen lassen.

Weitere Tutorien für die multiple lineare Regression

Wir immer geht natürlich mehr als ich hier Vorstellen kann. Du findest im Folgenden Tutorien, die mich hier in dem Kapitel inspiriert haben. Ich habe mich ja in diesem Kapitel auf die Durchführbarkeit in R und die allgemeine Verständlichkeit konzentriert. Es geht aber natürlich wie immer auch mathematischer…

• Wir funktioniert nun so eine lineare Regression und was sind den jetzt eigentlich die Koeffizienten \(\beta_0\) und \(\beta_1\) eigentlich? Hier gibt es die fantastische Seite Explained Visually | Ordinary Least Squares Regression, die dir nochmal erlaubt selbe mit Punkten in einem Scatterplot zu spielen und zu sehen wie sich dann die Regressionsgleichung ändert.
• Das Buch Tidy Modeling with R gibt nochmal einen tieferen Einblick in das Modellieren in R. Wir immer, es ist ein Vorschlag aber kein Muss.

45.3 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

pacman::p_load(tidyverse, magrittr, broom,
               see, performance, car, parameters,
               conflicted)
conflicts_prefer(magrittr::set_names)
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", 
                "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

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

45.4 Daten

Modellierung von Flöhen

Abbildung 45.9— In unseren Daten zu der Modellierung der Flöhe schauen wir auf verschiedene Ernährungsformen sowie Entwicklungsstadien und fragen uns, ob sich die Sprungweiten und andere Messwerte unterscheiden.

flea_model_tbl <- read_excel("data/fleas_model_data.xlsx") |> 
  mutate(feeding = as_factor(feeding),
         stage = as_factor(stage),
         bonitur = as.numeric(bonitur),
         infected = as.numeric(infected))

Tabelle 45.2— Tabelle der foo.

.id	feeding	stage	weight	hatched	K	Mg	CRP	Hb	BSG	jump_length	count_leg	count_leg_left	count_leg_right	bonitur	infected	M
1	sugar_water	adult	16.42	516.41	33.67	17.8	9.74	108.18	50.31	77.2	63	64	62	4	1	3971.61
2	sugar_water	adult	12.62	363.5	41.93	17.54	8.22	112.81	47.63	56.25	55	56	54	1	0	2064.13
3	sugar_water	adult	15.57	303.01	39.41	17.11	7.55	96.17	47.42	73.42	112	112	111	2	0	1763.43
…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…	…
46	ketchup	juvenile	7.18	429.18	32.45	19.53	15.67	81.7	15.6	83.38	423	420	425	4	1	20623.49
47	ketchup	juvenile	6.6	629.58	28.84	21.37	15.27	80.79	19.44	104.48	548	548	547	5	1	25149.36
48	ketchup	juvenile	4.19	192.66	26.61	20.06	13.95	80.32	16.07	130.18	869	872	866	5	0	21972.23

Theoretischer Datensatz

Abbildung 45.10— In unseren theoretischen Datensatz schauen wir uns die Körperlängen von verschiedenen Schlangen mit unterschiedlichen Einflussvariablen an. Wir schauen uns die Körperfarbe sowie die Regionen der Messung an.

Abschließend kommen wir noch zu einem theoretischen Datensatz dem ich Kéry (2010) entlehnt habe. Ich möchte nochmal zeigen, wie R intern mit der Modellmatrix die Berechnungen durchführt. Das hat mich interessiert und deshalb habe ich es einmal für mein Verständnis aufgeschrieben. Wenn es dich auch interessiert, dann kannst du in den entsprechenden Tabs weiter unten einmal nachschauen. Hier also erstmal ein sehr simpler theoretischer Datensatz zu Körperlängen von Schlangen. Wir müssen hier dann noch die Faktoren einmal umwandeln, damit wir auch später die Modellmatrix sauber vorliegen haben. Deshalb hier die explizite Benennung der Level in der Funktion factor().

snake_tbl <- read_xlsx("data/regression_data.xlsx", sheet = "theory_mult") |> 
  mutate(region = factor(region, levels = c("west", "nord")),
         color = factor(color, levels = c("schwarz", "rot", "blau")))

In der folgenden Tabelle ist der Datensatz snake_tbl nochmal dargestellt. Wir haben die Schlangenlänge svl als Messwert \(y\) sowie das Gewicht der Schlangen mass als kontinuierliche Einflussvariable sowie den Durchmesser der Schlange \(diameter\), die Sammelregion region und die Farbe der Schlangen color und als kategoriale Einflussvariablen mit unterschiedlichen Anzahlen an Gruppen. Die Region region ist also ein Faktor mit zwei Leveln und die Schlangenfarbe color ein Faktor mit drei Leveln. Ich nutze den Datensatz also einmal als einen Spieldatensatz um die Modellierung theoretisch besser nachvollziehen zu können.

Tabelle 45.3— Datensatz zu den Körperlängen als Messwert von acht Schlangen mit verschiedenen Einflussvariablen. Die Einflussvariablen haben verschiedene Eigenschaften und müssen daher unterschiedliche modelliert werden.

svl	mass	diameter	region	color
40	6	9	west	schwarz
45	8	3	west	schwarz
39	5	6	west	rot
51	7	7	nord	rot
52	9	5	nord	rot
57	11	4	nord	blau
58	12	10	nord	blau
49	10	11	nord	blau

45.5 Kausales Modell

45.5.1 Mehrkovariates Modell

Abbildung 45.11— Schemantisches multiples Modell mit einem Messwert \(Y\) und zwei kontinuierlichen Einflussvariablen als Kovariate \(c_1\) und Kovariate \(c_2\) dargestellt. [Zum Vergrößern anklicken]

Praktisch in R

cov2_fit <- lm(jump_length ~ weight + count_leg, data = flea_model_tbl) 

cov2_fit |> 
  coef() |> round(2)

(Intercept)      weight   count_leg 
      52.71        1.58        0.10 

Einmal die Ausgabe

cov2_fit |> summary()

Call:
lm(formula = jump_length ~ weight + count_leg, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-24.985 -12.823  -3.065   8.259  49.525 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 52.708201   6.389140   8.250 1.50e-10 ***
weight       1.584911   0.490155   3.233  0.00229 ** 
count_leg    0.098382   0.008789  11.193 1.35e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.78 on 45 degrees of freedom
Multiple R-squared:  0.7359,    Adjusted R-squared:  0.7242 
F-statistic: 62.71 on 2 and 45 DF,  p-value: 9.728e-14

Einmal die annotierte Ausgabe der Zusammenfassung

Abbildung 45.12— Annotierte Ausgabe der Funktion summary() aus einer linearen Modellanpassung mit der Funktion lm(). Die Ausgabe der Funktion teilt sich grob in drei informative Bereiche: Informationen zu den Residuen, Informationen zu den Koeffizienten und Informationen zu der Modelgüte. [Zum Vergrößern anklicken]

Theoretisch

cov2_tbl <- tibble(c_1 = rnorm(10, 0, 1),
                   c_2 = rnorm(10, 0, 1),
                   y = 2 + 
                       1 * c_1 + 
                       2 * c_2 + 
                           rnorm(10, 0, 0.001))

lm(y ~ c_1 + c_2, cov2_tbl) |> 
  coef() |> round(2)

(Intercept)         c_1         c_2 
          2           1           2 

Theoretisch in R

45.5.2 Mehrfaktorielles Modell

Abbildung 45.13— Schemantisches simples Modell mit einem Messwert \(Y\) und einer kategorialen Einflussvariablen als Faktor \(f_A\) mit drei Leveln sowie einer kategorialen Einflussvariablen als Faktor \(f_B\) mit zwei Leveln und deren Interaktion \(f_A \times f_B\) dargestellt. [Zum Vergrößern anklicken]

Praktisch in R

fac2_fit <- lm(jump_length ~ feeding + stage + stage:feeding, data = flea_model_tbl) 

fac2_fit |> 
  coef() |> round(2)

                 (Intercept)                 feedingblood 
                       78.98                        20.73 
              feedingketchup                stagejuvenile 
                       13.17                       -17.94 
  feedingblood:stagejuvenile feedingketchup:stagejuvenile 
                       34.26                        37.03 

fac2_fit |> summary()

Call:
lm(formula = jump_length ~ feeding + stage + stage:feeding, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-55.340 -16.617  -3.921  11.003 109.901 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     78.98      10.48   7.539 2.48e-09 ***
feedingblood                    20.73      14.81   1.400   0.1690    
feedingketchup                  13.17      14.81   0.889   0.3790    
stagejuvenile                  -17.94      14.81  -1.211   0.2326    
feedingblood:stagejuvenile      34.26      20.95   1.635   0.1095    
feedingketchup:stagejuvenile    37.03      20.95   1.768   0.0844 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 29.63 on 42 degrees of freedom
Multiple R-squared:  0.3158,    Adjusted R-squared:  0.2343 
F-statistic: 3.877 on 5 and 42 DF,  p-value: 0.005594

Abbildung 45.14— Annotierte Ausgabe der Funktion summary() aus einer linearen Modellanpassung mit der Funktion lm(). Die Ausgabe der Funktion teilt sich grob in drei informative Bereiche: Informationen zu den Residuen, Informationen zu den Koeffizienten und Informationen zu der Modelgüte. Die Erklärung der Schätzer (eng. estimate) kann in der folgenden Tabelle nachvollzogen werden. [Zum Vergrößern anklicken]

text

text

Tabelle 45.4— Zusammensetzung der Mittelwerte der Sprungweiten \(\bar{Y}\) für jede Faktorkombination aus der Ernährungsart und dem Entwicklungsstand aus den Koeffizienten des zweifaktoriellen Modells.

\(f_{feeding}\)	\(f_{stage}\)	\(\bar{Y}\)	\(\boldsymbol{\beta_{0}}\)	\(\boldsymbol{\beta_{blood}}\)	\(\boldsymbol{\beta_{ketchup}}\)	\(\boldsymbol{\beta_{juvenile}}\)	\(\boldsymbol{\beta_{blood \times juvenile}}\)	\(\boldsymbol{\beta_{ketchup \times juvenile}}\)
sugar	adult	78.90	78.98
sugar	juvenile	61.03	78.98			-17.94
blood	adult	99.71	78.98	20.73
blood	juvenile	116.03	78.98	20.73		-17.94	34.26
ketchup	adult	92.15	78.98		13.17
ketchup	juvenile	111.24	78.98		13.17	-17.94		37.03

text

Abbildung 45.15— Scatterplot der Sprungweiten der Hunde-, Katzen- und Fuchsflöhe. [Zum Vergrößern anklicken]

Theoretisch

fac2_tbl <- tibble(f_a = rep(gl(3, 5), 2),
                   f_b = gl(2, 15),
                   y = 2 + 
                       2 * as.numeric(f_a) + 
                       3 * as.numeric(f_b) + 
                           rnorm(15, 0, 0.001))

lm(y ~ f_a + f_b + f_a:f_b, fac2_tbl) |> 
  coef() |> round(2)

(Intercept)        f_a2        f_a3        f_b2   f_a2:f_b2   f_a3:f_b2 
          7           2           4           3           0           0 

Theoretisch in R

45.5.3 Kombinierte Modelle

Wenn wir eine Kovariate mit einem Faktor kombinieren, dann kippen wir den Intercept der Graden um die Steigung der Kovariate. Die Level der Faktoren liegen dann auf den Mittelwerten der Kovariaten für das entsprechende Level des Fakotrs. Wenn wir dann noch eine Interaktion hinzunehmen, dann erlauben wir jedem Level des Faktors eine eigene Steigung der Kovariaten.

\(f_A + c_1\)

Hier rechnen wir eigentlich eine klassische einfaktorielle ANCOVA (eng. analysis of covariance)

Abbildung 45.16— Schemantisches simples Modell mit einem Messwert \(Y\) und einer kategorialen Einflussvariablen als Faktor \(f_A\) mit drei Leveln sowie einer kontinuierlichen Einflussvariable als Kovariate \(c_1\) dargestellt. [Zum Vergrößern anklicken]

Praktisch in R

Ohne Interaktion

cov1_fac1_fit <- lm(jump_length ~ feeding + K, data = flea_model_tbl) 

cov1_fac1_fit |> 
  coef() |> round(2)

   (Intercept)   feedingblood feedingketchup              K 
         58.48          39.31          34.64           0.34 

cov1_fac1_fit |> summary()

Call:
lm(formula = jump_length ~ feeding + K, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.952 -17.163  -5.820   8.661 115.648 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     58.4806    17.3429   3.372 0.001564 ** 
feedingblood    39.3101    10.8621   3.619 0.000759 ***
feedingketchup  34.6434    11.4092   3.036 0.004014 ** 
K                0.3420     0.4633   0.738 0.464270    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 30.22 on 44 degrees of freedom
Multiple R-squared:  0.2543,    Adjusted R-squared:  0.2035 
F-statistic: 5.002 on 3 and 44 DF,  p-value: 0.004523

flea_model_tbl |> 
  group_by(feeding) |> 
  summarise(mean(jump_length),
            mean(K))

# A tibble: 3 × 3
  feeding     `mean(jump_length)` `mean(K)`
  <fct>                     <dbl>     <dbl>
1 sugar_water                70.0      33.7
2 blood                     108.       29.5
3 ketchup                   102.       25.1

Tabelle 45.5— Zusammensetzung der Mittelwerte der Sprungweiten \(\bar{Y}\) für jede Faktorkombination aus der Ernährungsart und dem Entwicklungsstand aus den Koeffizienten des zweifaktoriellen Modells.

\(f_{feeding}\)	\(\bar{c}_{K}\)	\(\bar{Y}\)	\(\boldsymbol{\beta_{0}}\)	\(\boldsymbol{\beta_{blood}}\)	\(\boldsymbol{\beta_{ketchup}}\)	\(\boldsymbol{\bar{c}_{K}\cdot\beta_{K}}\)
sugar	33.69	70.01	58.48			\(11.45 = 33.69 \cdot 0.34\)
blood	29.47	107.87	58.48	39.31		\(10.02 = 29.47 \cdot 0.34\)
ketchup	25.06	101.69	58.48		34.64	\(8.52 = 25.06 \cdot 0.34\)

Abbildung 45.17— Scatterplot der Sprungweiten der Hunde-, Katzen- und Fuchsflöhe. [Zum Vergrößern anklicken]

Mit Interaktion

cov1_fac1_fit <- lm(jump_length ~ feeding + K + feeding:K, data = flea_model_tbl) 

cov1_fac1_fit |> summary()

Call:
lm(formula = jump_length ~ feeding + K + feeding:K, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-52.799 -15.960  -4.885  10.009 113.476 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)
(Intercept)       44.1255    48.1386   0.917    0.365
feedingblood      44.5900    51.4164   0.867    0.391
feedingketchup    81.6294    55.2411   1.478    0.147
K                  0.7681     1.4110   0.544    0.589
feedingblood:K    -0.1181     1.5168  -0.078    0.938
feedingketchup:K  -1.7283     1.7520  -0.987    0.330

Residual standard error: 30.23 on 42 degrees of freedom
Multiple R-squared:  0.2877,    Adjusted R-squared:  0.2029 
F-statistic: 3.393 on 5 and 42 DF,  p-value: 0.01155

Tabelle 45.6— Zusammensetzung der Mittelwerte der Sprungweiten \(\bar{Y}\) für jede Faktorkombination aus der Ernährungsart und dem Entwicklungsstand aus den Koeffizienten des zweifaktoriellen Modells.

\(f_{feeding}\)	\(\bar{c}_{K}\)	\(\bar{Y}\)	\(\boldsymbol{\beta_{0}}\)	\(\boldsymbol{\beta_{blood}}\)	\(\boldsymbol{\beta_{ketchup}}\)	\(\boldsymbol{\bar{c}_{K}\cdot\beta_{K}}\)	\(\boldsymbol{\bar{c}_{K}\cdot\beta_{blood\times K}}\)	\(\boldsymbol{\bar{c}_{K}\cdot\beta_{ketchup\times K}}\)
sugar	33.69	70.01	44.13			\(25.94 = 33.69 \cdot 0.77\)
blood	29.47	107.87	44.13	44.59		\(22.69 = 29.47 \cdot 0.77\)	\(-3.54 = 29.47 \cdot -0.12\)
ketchup	25.06	101.69	44.13		81.63	\(19.29 = 25.06 \cdot 0.77\)		\(-43.35 = 25.06 \cdot -1.73\)

Abbildung 45.18— Scatterplot der Sprungweiten der Hunde-, Katzen- und Fuchsflöhe. foo. [Zum Vergrößern anklicken]

Theoretisch

cov1_fac1_tbl <- tibble(c_1 = rnorm(15, 0, 1),
                        f_a = gl(3, 5),
                        y = 2 + 
                            1 * c_1 + 
                            2 * as.numeric(f_a) + 
                                rnorm(15, 0, 0.001))

lm(y ~ f_a + c_1, cov1_fac1_tbl) |> 
  coef() |> round(2)

(Intercept)        f_a2        f_a3         c_1 
          4           2           4           1 

Theoretisch in R

\(c_1 + f_A + f_B\)

Hier rechnen wir eigentlich eine klassische zweifaktorielle ANCOVA (eng. analysis of covariance)

Abbildung 45.19— Schemantisches simples Modell mit einem Messwert \(Y\) und einer kontinuierlichen Einflussvariable als Kovariate \(c_1\) sowie einer kategorialen Einflussvariablen als Faktor \(f_A\) mit drei Leveln und einer kategorialen Einflussvariablen als Faktor \(f_B\) mit zwei Leveln dargestellt. [Zum Vergrößern anklicken]

Praktisch in R

cov1_fac2_fit <- lm(jump_length ~ weight + feeding + stage, data = flea_model_tbl) 

cov1_fac2_fit |> 
  coef() |> round(2)

   (Intercept)         weight   feedingblood feedingketchup  stagejuvenile 
         33.17           2.27          41.03          34.73          24.41 

cov1_fac2_fit |> summary()

Call:
lm(formula = jump_length ~ weight + feeding + stage, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.448 -12.133  -4.554   7.315 114.062 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      33.169     20.130   1.648 0.106693    
weight            2.267      1.219   1.859 0.069850 .  
feedingblood     41.033     10.548   3.890 0.000343 ***
feedingketchup   34.735     10.538   3.296 0.001970 ** 
stagejuvenile    24.406     13.121   1.860 0.069725 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 29.44 on 43 degrees of freedom
Multiple R-squared:  0.3082,    Adjusted R-squared:  0.2439 
F-statistic:  4.79 on 4 and 43 DF,  p-value: 0.002764

Theoretisch

cov1_fac2_tbl <- tibble(c_1 = rnorm(30, 0, 1),
                        f_a = rep(gl(3, 5), 2),
                        f_b = gl(2, 15),
                        y = 2 + 
                            1 * c_1 + 
                            2 * as.numeric(f_a) + 
                            3 * as.numeric(f_b) + 
                                rnorm(15, 0, 0.001))

lm(y ~ c_1 + f_a + f_b, cov1_fac2_tbl) |> 
  coef() |> round(2)

(Intercept)         c_1        f_a2        f_a3        f_b2 
          7           1           2           4           3 

Theoretisch in R

Im letzten Schritt wollen wir uns einmal das volle Modell anschauen. Wir bauen uns ein Modell mit allen Variablen in dem Datensatz snake_tbl. Die Funktion model.matrix() gibt uns die Modelmatrix wieder.

model.matrix(svl ~ mass + region + color, data = snake_tbl) |> as_tibble() 

# A tibble: 8 × 5
  `(Intercept)`  mass regionnord colorrot colorblau
          <dbl> <dbl>      <dbl>    <dbl>     <dbl>
1             1     6          0        0         0
2             1     8          0        0         0
3             1     5          0        1         0
4             1     7          1        1         0
5             1     9          1        1         0
6             1    11          1        0         1
7             1    12          1        0         1
8             1    10          1        0         1

In der ersten Spalte ist der Intercept angegeben, danach folgt dann die Spalte mass als kontenuierliche Variable. In der Spalte regionnord steht die Dummykodierung für die Variable region. Die ersten drei Schlangen kommen nicht aus der Region nord und werden deshalb mit einen Wert von 0 versehen. Die nächsten vier Schlangen kommen aus der Region nord und erhalten daher eine 1 in der Spalte. Die nächsten beiden Spalten sind etwas komplizierter. Die Spalten colorrot und colorblau geben jeweils an, ob die Schlange das Level rot hat oder blau oder keins von beiden. Wenn die Schlange weder rot noch blau ist, dann sind beide Spalten mit einer 0 versehen. Dann ist die Schlange schwarz.

Wir können die Modellmatrix auch mathematisch schreiben und die \(y\) Spalte für das Outcome svl ergänzen. Eben so ergänzen wir die \(\beta\)-Werte als mögliche Koeefizienten aus der linearen Regression sowie die Residuen als Abweichung von der gefitteten Gerade.

\[ \begin{pmatrix} 40 \\ 45 \\ 39 \\ 50 \\ 52 \\ 57 \\ 58 \\ 59 \\ \end{pmatrix} = \begin{pmatrix} 1 & 6 & 0 & 0 & 0 \\ 1 & 8 & 0 & 0 & 0\\ 1 & 5 & 0 & 1 & 0\\ 1 & 7 & 1 & 1 & 0\\ 1 & 9 & 1 & 1 & 0\\ 1 & 11& 1 & 0 & 1\\ 1 & 12& 1 & 0 & 1\\ 1 & 10& 1 & 0 & 1\\ \end{pmatrix} \times \begin{pmatrix} \beta_0 \\ \beta_{mass} \\ \beta^{region}_{nord} \\ \beta^{color}_{rot} \\ \beta^{color}_{blau} \\ \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \epsilon_3 \\ \epsilon_4 \\ \epsilon_5 \\ \epsilon_6 \\ \epsilon_7 \\ \epsilon_8 \\ \end{pmatrix} \]

Jetzt brauche wir die Koeffizienten aus dem linearen Modell, welches wir wie folgt fitten. Wir nutzen dann die Funktion coef() um uns die Koeffizienten wiedergeben zu lassen.

fit_4 <- lm(svl ~ mass + region + color, data = snake_tbl) 
fit_4 |> coef() |> round(2)

(Intercept)        mass  regionnord    colorrot   colorblau 
      25.00        2.50        5.00        1.50       -2.83 

Die Funktion residuals() gibt uns die Residuen der Geraden aus dem Objekt fit_4 wieder.

fit_4 |> residuals() |> round(2)

    1     2     3     4     5     6     7     8 
 0.00  0.00  0.00  2.00 -2.00  2.33  0.83 -3.17 

Wir können jetzt die Koeffizienten ergänzen mit \(\beta_0 = 25\) für den Intercept. Weiter ergänzen wir die Koeffizienten für mass mit \(\beta_{mass}=2.5\), für Region und das Level nord mit \(\beta^{region}_{nord} = 5\), für die Farbe und das Level rot mit \(\beta^{color}_{rot} = 1.5\) und für die Farbe und das Level blau mit \(\beta^{color}_{blau} = -2.83\). Ebenfalls setzen wir Werte für die Residuen für jede der Beobachtungen in die Gleichung ein.

\[ \begin{pmatrix} 40 \\ 45 \\ 39 \\ 50 \\ 52 \\ 57 \\ 58 \\ 59 \\ \end{pmatrix} = \begin{pmatrix} 25 & \phantom{0}6 \cdot 2.5 & 0 \cdot 5 & 0 \cdot 1.5& 0 \cdot -2.83 \\ 25 & \phantom{0}8 \cdot 2.5 & 0 \cdot 5 & 0 \cdot 1.5& 0 \cdot -2.83\\ 25 & \phantom{0}5 \cdot 2.5 & 0 \cdot 5 & 1 \cdot 1.5& 0 \cdot -2.83\\ 25 & \phantom{0}7 \cdot 2.5 & 1 \cdot 5 & 1 \cdot 1.5& 0 \cdot -2.83\\ 25 & \phantom{0}9 \cdot 2.5 & 1 \cdot 5 & 1 \cdot 1.5& 0 \cdot -2.83\\ 25 & 11\cdot 2.5 & 1 \cdot 5 & 0 \cdot 1.5& 1 \cdot -2.83\\ 25 & 12\cdot 2.5 & 1 \cdot 5 & 0 \cdot 1.5& 1 \cdot -2.83\\ 25 & 10\cdot 2.5 & 1 \cdot 5 & 0 \cdot 1.5& 1 \cdot -2.83\\ \end{pmatrix} + \begin{pmatrix} \phantom{+}0.00 \\ \phantom{+}0.00 \\ \phantom{+}0.00 \\ +2.00 \\ -2.00 \\ +2.33 \\ +0.83 \\ -3.17 \\ \end{pmatrix} \]

Wir können jetzt diese gewaltige Sammlung an Matrixen einmal auflösen. Steht denn nun wirklich rechts das Gleiche wie links von der Gleichung? Wir bauen uns die Zahlen von der rechten Seite der Gleichung einmal in R nach und schauen auf das Ergebnis.

c(25 +  6*2.5 + 0*5 + 0*1.5 + 0*-2.83 + 0.00,
  25 +  8*2.5 + 0*5 + 0*1.5 + 0*-2.83 + 0.00,
  25 +  5*2.5 + 0*5 + 1*1.5 + 0*-2.83 + 0.00,
  25 +  7*2.5 + 1*5 + 1*1.5 + 0*-2.83 + 2.00,
  25 +  9*2.5 + 1*5 + 1*1.5 + 0*-2.83 - 2.00,
  25 + 11*2.5 + 1*5 + 0*1.5 + 1*-2.83 + 2.33,
  25 + 12*2.5 + 1*5 + 0*1.5 + 1*-2.83 + 0.83,
  25 + 10*2.5 + 1*5 + 0*1.5 + 1*-2.83 - 3.17) 

[1] 40 45 39 51 52 57 58 49

Die Zahlen, die wir rauskriegen, sind die gleichen Werte die unser Outcome \(y\) hat. Was haben wir gelernt?

• In einem Modell gibt es immer ein Faktorlevel weniger als ein Faktor Level hat. Die Information des alphanumerisch ersten Levels steckt dann mit in dem Intercept.
• In einem Modell geht eine kontinuierliche Variable als eine Spalte mit ein.
• In einem Modell gibt es immer nur eine Spalte für die Residuen und damit nur eine Residue für jede Beobachtung, egal wie viele Variablen ein Modell hat.

\(c_1 + c_2 + f_A + f_B\)

Abbildung 45.20— Schemantisches simples Modell mit einem Messwert \(Y\) und einer kontinuierlichen Einflussvariable als Kovariate \(c_1\) sowie einer kategorialen Einflussvariablen als Faktor \(f_A\) mit drei Leveln und einer kategorialen Einflussvariablen als Faktor \(f_B\) mit zwei Leveln dargestellt. [Zum Vergrößern anklicken]

Praktisch in R

cov2_fac2_fit <- lm(jump_length ~ weight + count_leg + feeding + stage, data = flea_model_tbl) 

cov2_fac2_fit |> 
  coef() |> round(2)

   (Intercept)         weight      count_leg   feedingblood feedingketchup 
         57.16           1.00           0.10          14.79           5.16 
 stagejuvenile 
         -9.37 

cov2_fac2_fit |> summary()

Call:
lm(formula = jump_length ~ weight + count_leg + feeding + stage, 
    data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-25.427 -10.507  -2.446   8.093  47.736 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)    57.15639   11.74218   4.868 1.63e-05 ***
weight          1.00219    0.70718   1.417   0.1638    
count_leg       0.09595    0.01009   9.512 4.89e-12 ***
feedingblood   14.79263    6.61231   2.237   0.0306 *  
feedingketchup  5.15633    6.76099   0.763   0.4499    
stagejuvenile  -9.37311    8.27589  -1.133   0.2638    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 16.77 on 42 degrees of freedom
Multiple R-squared:  0.7807,    Adjusted R-squared:  0.7546 
F-statistic:  29.9 on 5 and 42 DF,  p-value: 8.061e-13

Theoretisch

cov2_fac2_tbl <- tibble(c_1 = rnorm(30, 0, 1),
                        c_2 = rnorm(30, 0, 1),
                        f_a = rep(gl(3, 5), 2),
                        f_b = gl(2, 15),
                        y = 2 + 
                            1 * c_1 + 
                            -1 * c_2 + 
                            2 * as.numeric(f_a) + 
                            3 * as.numeric(f_b) + 
                                rnorm(15, 0, 0.001))

lm(y ~ c_1 + c_2 + f_a + f_b, cov2_fac2_tbl) |> 
  coef() |> round(2)

(Intercept)         c_1         c_2        f_a2        f_a3        f_b2 
          7           1          -1           2           4           3 

Theoretisch in R

45.6 Prädiktives Modell

45.7 Weiteres von Interesse

45.7.1 Die Einheit

Abbildung 45.21— foo

Manche Flöhe haben hohe Werte der Midi-Chlorianer oder einfach M-Werte. Damit sind diese Flöhe insbesondere sensitiv für die Macht und somit auch in etwa Jediflöhe.

lm(jump_length ~ M, data = flea_model_tbl) |> 
  summary()

Call:
lm(formula = jump_length ~ M, data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-56.163 -21.533  -5.392  16.543 136.037 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.041e+01  7.356e+00  10.932 2.22e-14 ***
M           1.515e-03  6.719e-04   2.254    0.029 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32.48 on 46 degrees of freedom
Multiple R-squared:  0.09949,   Adjusted R-squared:  0.07991 
F-statistic: 5.082 on 1 and 46 DF,  p-value: 0.02898

45.7.2 Korrelation

cov2_fit <- lm(jump_length ~ count_leg + count_leg_left + count_leg_right, data = flea_model_tbl) 

cov2_fit |> 
  coef() |> round(2)

    (Intercept)       count_leg  count_leg_left count_leg_right 
          75.11          -15.07            7.79            7.37 

 1.37 -1.15 -0.12 

[1] 0.1

Abbildung 45.22— Scatterplot der Sprungweite in [cm] und dem Gewicht in [mg] von sieben Hundeflöhen. Die lineare Grade und die Gradengleichung wurden ergänzt. [Zum Vergrößern anklicken]

cov2_fit |> summary()

Call:
lm(formula = jump_length ~ count_leg + count_leg_left + count_leg_right, 
    data = flea_model_tbl)

Residuals:
    Min      1Q  Median      3Q     Max 
-29.307 -14.465  -4.338  11.832  46.265 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       75.111      5.243  14.327   <2e-16 ***
count_leg        -15.069     12.130  -1.242    0.221    
count_leg_left     7.790      6.028   1.292    0.203    
count_leg_right    7.372      6.128   1.203    0.235    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.51 on 44 degrees of freedom
Multiple R-squared:  0.6892,    Adjusted R-squared:  0.668 
F-statistic: 32.52 on 3 and 44 DF,  p-value: 3.082e-11

cov2_fit <- lm(jump_length ~ count_leg, data = flea_model_tbl) 

cov2_fit |> 
  coef() |> round(2)

(Intercept)   count_leg 
      70.33        0.09 

cov2_fit |> summary()

Call:
lm(formula = jump_length ~ count_leg, data = flea_model_tbl)

Residuals:
   Min     1Q Median     3Q    Max 
-25.46 -14.07  -5.00  10.05  53.19 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 70.325976   3.663778  19.195  < 2e-16 ***
count_leg    0.091047   0.009323   9.765 8.61e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 19.52 on 46 degrees of freedom
Multiple R-squared:  0.6746,    Adjusted R-squared:  0.6675 
F-statistic: 95.36 on 1 and 46 DF,  p-value: 8.612e-13

45.7.3 Adjustierung für Confounder

Im folgenden Abschnitt wollen wir einmal auf Confounder eingehen. Was sind Confounder? Zum einen gibt es kein gutes deutsches Wort für Confounder. Du könntest Confounder in etwa mit Verzerrer oder Störfakor übersetzen. Zum Anderen können Confounder nur zuammen mit anderen Vriabln in einer multiplen Regression auftreten. Das heißt, wir brauchen mindestens zwei Variablen in einer Regression. Ein Confounder verursacht einen Effekt, den wir eigentlich nicht so erwartet hätten. Wir wollen eigentlich einen Effekt schätzen, aber der Effekt ist viel größer oder kleiner, da der Effekt eigentlich von einder anderen Variable verursacht oder aber verdeckt wird. Wir können einen Confounder in beide Richtungen haben. Wichtig ist hierbei, das wir eigentlich nicht an dem Effekt des Confounders interessiert sind. Wir wollen uns zum Beispiel den Effekt einer Düngung auf das Trockengewicht anschauen, aber der Effekt den wir beobachten wird durch die unterschiedlichen Pflanzorte verursacht.

Schauen wir uns das ganze mal für das Beispiel des Zusammenhangs von dem Flohgewicht mit der Sprungweite an. Dafür benötigen wir den Datensatz flea_dog_cat_length_weight.csv.

model_tbl <- read_csv2("data/flea_dog_cat_length_weight.csv") |>
  select(animal, sex, weight, jump_length) |> 
  mutate(animal = as_factor(animal),
         sex = as_factor(sex))

In der Tabelle 45.7 ist der Datensatz model_tbl nochmal dargestellt.

Tabelle 45.7— Datensatz für die Confounder Adjustierung. Die Variable jump_length ist das \(y\) und die Variable weight das \(x\) von Interesse.

animal	sex	weight	jump_length
cat	male	6.02	15.79
cat	male	5.99	18.33
cat	male	8.05	17.58
cat	male	6.71	14.09
cat	male	6.19	18.22
cat	male	8.18	13.49
…	…	…	…
fox	female	8.04	27.81
fox	female	9.03	24.02
fox	female	7.42	24.53
fox	female	9.26	24.35
fox	female	8.85	24.36
fox	female	7.89	22.13

Wir können uns jetzt drei verschiedene Modelle anschauen.

1. Das Modell jump_length ~ weight. Wir modellieren die Abhängigkeit von der Sprungweite von dem Gewicht der Flöhe über alle anderen Variablen hinweg.
2. Das Modell jump_length ~ weight + animal. Wir modellieren die Abhängigkeit von der Sprungweite von dem Gewicht der Flöhe und berücksichtigen die Tierart des Flohes. Ignorieren aber das Geschlecht des Flohes.
3. Das Modell jump_length ~ weight + animal + sex. Wir modellieren die Abhängigkeit von der Sprungweite von dem Gewicht der Flöhe und berücksichtigen alle Variablen, die wir gemessen haben.

Hierbei ist wichtig, dass es natürlich auch Confounder geben kann, die wir gar nicht in den Daten erhoben haben. Also, dass sich Gruppen finden lassen, die wir gar nicht als Variable mit in den Daten erfasst haben. Hier könnte eine Hauptkomponentenanalyse helfen.

(a) jump_length ~ weight

(b) jump_length ~ weight + animal

(c) jump_length ~ weight + animal + sex

Abbildung 45.23— Darstellung des counfounder Effekts anhand des Zusammenhangs der Sprungweite in [cm] und dem Gewicht von Flöhen [mg].

In Abbildung 45.23 (a) sehen wir die blaue Linie als die Ausgabe von dem Modell jump_length ~ weight. Wir sehen, dass mit dem Anstieg des Gewichtes der Flöhe auch die Sprungweite sich erhöht. Wir würden annehmen, dass wir hier einen signifikanten Unterschied vorliegen haben. Schauen wir uns die Koeffizienten des Modells aus dem lm() einmal an.

lm(jump_length ~ weight, data = model_tbl) |> 
  model_parameters()

Parameter   | Coefficient |   SE |        95% CI | t(598) |      p
------------------------------------------------------------------
(Intercept) |        9.79 | 0.77 | [8.28, 11.30] |  12.73 | < .001
weight      |        1.34 | 0.09 | [1.16,  1.53] |  14.16 | < .001

Wir sehen, dass wir einen Effekt des Gewichts auf die Sprungweite von \(1.34\) vorliegen haben. Auch ist der Effekt und damit die Steigung signifikant.

Erweitern wir nun das Modell um die Tierart und erhalten die Abbildung 45.23 (b). Zum einen sehen wir, dass der globale Effekt nicht so ganz stimmen kann. Die Tierarten haben alle einen unterschiedlichen starken Effekt von dem Gewicht auf die Sprungweite. Auch hier fitten wir einmal das lineare Modell und schauen uns die Koeffizienten an.

lm(jump_length ~ weight + animal, data = model_tbl) |> 
  model_parameters()

Parameter    | Coefficient |   SE |       95% CI | t(596) |      p
------------------------------------------------------------------
(Intercept)  |        7.82 | 0.62 | [6.60, 9.03] |  12.64 | < .001
weight       |        1.27 | 0.07 | [1.13, 1.42] |  17.08 | < .001
animal [dog] |        2.62 | 0.26 | [2.12, 3.13] |  10.23 | < .001
animal [fox] |        4.97 | 0.26 | [4.47, 5.48] |  19.38 | < .001

Der Effekt des Gewichtes ist hier in etwa gleich geblieben. Wir sehen aber auch, dass die Sprungweite anscheinend bei Hunden und Füchsen höher ist. Daher haben wir hier einen Effekt on der Tierart. Wir adjustieren für den Confoundereffekt durch die Tierart und erhalten einen besseren Effektschätzer für das Gewicht.

In der letzten Abbildung 45.23 (c) wollen wir nochmal das Geschlecht der Flöhe mit in das Modell nehmen. Wir sehen, dass in diesem Fall, das Gewicht gar keinen Einfluss mehr auf die Sprungweite hat. Den Effekt des Gewichtes war nur der Effekt des unterschiedlichen Gewichtes der weiblichen und männlichen Flöhe. Schauen wir auch hier nochmal in des Modell.

lm(jump_length ~ weight + animal + sex, data = model_tbl) |> 
  model_parameters() |> 
  mutate(Coefficient = round(Coefficient, 2))

# Fixed Effects

Parameter   | Coefficient |   SE |         95% CI | t(595) |      p
-------------------------------------------------------------------
(Intercept) |       15.42 | 0.59 | [14.27, 16.58] |  26.23 | < .001
weight      |        0.01 | 0.08 | [-0.16,  0.17] |   0.07 | 0.942 
animaldog   |        2.61 | 0.19 | [ 2.23,  2.99] |  13.50 | < .001
animalfox   |        5.19 | 0.19 | [ 4.81,  5.57] |  26.75 | < .001
sexfemale   |        4.89 | 0.23 | [ 4.44,  5.34] |  21.25 | < .001

Nun können wir sehen, dass von unserem ursprünglichen Effekt von dem Gewicht auf die Sprungweite nichts mehr übrigbleibt. Wir haben gar keinen Effekt von dem Gewicht auf die Sprungweite vorliegen. Was wir gesehen haben, war der Effekt des Gewichtes der unterschiedlichen Geschlechter. Wenn wir unser Modell für die Confounder animal und sex adjustieren, haben wir einen unverzerrten Schätzer für weight. Wichtig ist nochmal, dafür müssen wir natürlich auch alle variablen erhoben haben. Hätten wir das Geschlecht der Flöhe nicht bestimt, hätten wir hier eventuell einen falschen Schluß getroffen.

45.7.4 Variance inflation factor (VIF)

Wenn wir sehr viele Daten erheben, dann kann es sein, dass die Variablen stark miteinander korrelieren. Wir können aber stark mieinander korrelierte Variablen nicht zusammen in ein Modell nehmen. Die Effekte der beiden Variablen würden sich gegenseitig aufheben. Wir hätten eigentlich zwei getrennt signifikante Variablen. Nehmen wir aber beide Variablen mit ins Modell, sind beide Variablen nicht mehr signifikant.

Wir kürzen hier stark ab. Wenn du mehr über Variablen in einem Modell wissen willst, gibt es dann in den Kapitel zur Variablen Selektion mehr Informationen. Wenn du nur Outcomes und Blöcke vorliegen hast, dann ist das VIF für dich uninteressant.

Um Variablen zu finden, die sich sehr ähnlich verhalten, können wir den Variance inflation factor (VIF) nutzen. Wir können den VIF jedoch nicht für kategoriale Daten verwenden. Statistisch gesehen würde es keinen Sinn machen. Die Anwendung ist recht einfach. Wir fitten als erstes unser Modell und dann können wir die Funktion vif() aus dem R Paket {car} verwenden.

gummi_tbl <- read_excel("data/gummibears.xlsx")
model <- lm(height ~ semester + age + count_color + count_bears, data = gummi_tbl)

Wir sehen im Folgenden das Eregbnis der Funktion vif(). Wenn der berechnete Wert für den VIF größer als 5 ist, dann liegt mit der Variable ein Problem vor. Wir vermuten dann, dass eine andere variable sehr stark mit dieser Variable korreliert. Wir sehen an den Werten keine Aufälligkeiten.

vif(model)

   semester         age count_color count_bears 
   1.051656    1.030493    1.121956    1.173528 

Wir können uns mit der Funcktion check_model() aus dem R Paket {performance} auch die Unsicherheit mit angeben lassen. In unserem Beispiel hieft dies gerade nicht sehr viel weiter. Wir bleiben bei den geschätzen Werten und ignorieren das Intervall.

Abbildung 45.24— Graphische Darstellung des VIF mit der Funktion check_model().

In Abbildung 45.24 sehen wir nochmal die Visualisierung. Wie immer, manchmal helfen solche Abbildungen, manchmal verwirren die Abbildungen mehr. Wir konzentrieren uns hier auf die Werte des VIF’s und ignorieren die Streuung.

45.7.5 Vergleich von Modellen

Im Folgenden wollen wir einmal verschiedene Modelle miteinander Vergleichen und uns statistisch wiedergeben lassen, was das beste Modell ist. Und hier holen wir auch einmal kurz Luft, denn wir entschieden nur was das statistisch beste Modell ist. Es kann sein, dass ein Modell biologisch mehr Sinn macht und nicht auf Platz 1 der statistischen Maßzahlen steht. Das ist vollkommen in Ordnung. Du musst abweägen, was für sich das beste Modell ist. Im Zweifel komme ruhig nochmal in meine statistische Beratung oder schreibe mir eine Mail.

Wir kürzen hier stark ab bzw. gehen nicht im Detail auf jedes Gütekriterium ein. Wichtig ist, dass du Modelle statistisch vergleichen kannst. Bedenke immmer, dass die statistische Bewertung nicht immer ausreicht! Auch was die Biologie über das Modell sagt ist wichtig.

Wir bauchen uns jetzt fünf Modelle von fit_1 bis fit_5. Jedes dieser Modelle hat andere \(x\) aber häufig das gleiche Outcome y. In dem Beispiel am Ende des Kapitels nutzen wir auch verschiedene \(y\) dafür aber dann gleiche \(x\) in dem Modell. Im Weiteren sortieren wir die Modelle von einfach nach komplex. Ich versuche immmer das einfachste Modell fit_1 zu nennen bzw. eher die niedrige Nummer zu geben. Im Idealfall benennst du die Modellobjekte nach den Modellen, die in en Objekten gespeichert sind. Oft sind die Modelle aber sehr groß und die Objekte der Fits haben dann sehr lange Namen.

fit_1 <- lm(jump_length ~ animal, data = model_tbl)
fit_2 <- lm(jump_length ~ animal + sex, data = model_tbl)
fit_3 <- lm(jump_length ~ animal + sex + weight, data = model_tbl)
fit_4 <- lm(jump_length ~ animal + sex + sex:weight, data = model_tbl)
fit_5 <- lm(log(jump_length) ~ animal + sex, data = model_tbl)

Als Ergänzung zum Bestimmtheitsmaß \(R^2\) wollen wir uns noch das Akaike information criterion (abk. \(AIC\)) anschauen. Du kannst auch das \(R^2\) bzw. das \(R^2_{adj}\) für die Modellauswahl nehmen. Das \(AIC\) ist neuer und auch für komplexere Modelle geeignet. Es gilt hierbei, je kleiner das \(AIC\) ist, desto besser ist das \(AIC\). Wir wollen also Modelle haben, die ein kleines \(AIC\) haben. Wir gehen jetzt nicht auf die Berechnung der \(AIC\)’s für jedes Modell ein. Wir erhalten nur ein \(AIC\) für jedes Modell. Die einzelnen Werte des \(AIC\)’s sagen nichts aus. Ein \(AIC\) ist ein mathematisches Konstrukt. Wir können aber verwandte Modelle mit dem \(AIC\) untereinander vergleichen. Daher berechnen wir ein \(\Delta\) über die \(AIC\). Dafür nehmen wir das Modell mit dem niedrigsten \(AIC\) und berechnen die jeweiligen Differenzen zu den anderen \(i\) Modellen. In unserem Beispiel ist \(i\) dann gleich fünf, da wir fünf Modelle haben.

\[ \Delta_i = AIC_i - AIC_{min} \]

• wenn \(\Delta_i < 2\), gibt es keinen Unterschied zwischen den Modellen. Das \(i\)-te Modell ist genauso gut wie das Modell mit dem \(AIC_{min}\).
• wenn \(2 < \Delta_i < 4\), dann gibt es eine starke Unterstützung für das \(i\)-te Modell. Das \(i\)-te Modell ist immer noch ähnlich gut wie das \(AIC_{min}\).
• wenn \(4 < \Delta_i < 7\), dann gibt es deutlich weniger Unterstützung für das \(i\)-te Modell;
• Modelle mit \(\Delta_i > 10\) sind im Vergleich zu dem besten \(AIC\) Modell nicht zu verwenden.

Nehmen wir ein \(AIC_1 = AIC_{min} = 100\) und \(AIC_2\) ist \(100,7\) an. Dann ist \(\Delta_2=0,7<2\), so dass es keinen wesentlichen Unterschied zwischen den Modellen gibt. Wir können uns entscheiden, welches der beiden Modelle wir nehmen. Hier ist dann wichtig, was auch die Biologie sagt oder eben andere Kriterien, wie Kosten und Nutzen. Wenn wir ein \(AIC_1 = AIC_{min} = 100000\) und \(AIC_2\) ist \(100700\) vorliegen haben, dann ist \(\Delta_2 = 700 \gg 10\), also gibt es keinen Grund für das \(2\)-te Modell. Das \(2\)-te Modell ist substantiell schlechter als das erste Modell. Mehr dazu kannst du unter Multimodel Inference: Understanding AIC and BIC in Model Selection nachlesen.

Wir können das \(\Delta_i\) auch in eine Wahrscheinlichkeit umrechnen. Wir können \(p_i\) berechnen und damit die relative (im Vergleich zu \(AIC_{min}\)) Wahrscheinlichkeit, dass das \(i\)-te Modell den AIC minimiert.

\[ p_i = \exp\left(\cfrac{-\Delta_i}{2}\right) \]

Zum Beispiel entspricht \(\Delta_i = 1.5\) einem \(p_i\) von \(0.47\) (ziemlich hoch) und ein \(\Delta_i = 15\) entspricht einem \(p_i =0.0005\) (ziemlich niedrig). Im ersten Fall besteht eine Wahrscheinlichkeit von 47%, dass das \(i\)-te Modell tatsächlich eine bessere Beschreibung ist als das Modell, das \(AIC_{min}\) ergibt, und im zweiten Fall beträgt diese Wahrscheinlichkeit nur 0,05%. Wir können so einmal nachrechnen, ob sich eine Entscheidung für ein anderes Modell lohnen würde. Neben dem \(AIC\) gibt es auch das Bayesian information criterion (\(BIC\)). Auch beim \(BIC\) gilt, je kleiner das BIC ist, desto besser ist das BIC.

Du siehst schon, es gibt eine Reihe von Möglichkeiten sich mit der Güte oder Qualität eines Modells zu beschäftigen. Wir nutzen die Funktion model_performance() um uns die Informationen über die Güte eines Modells wiedergeben zu lassen. Im folgenden Codeblock habe ich mich nur auf das \(AIC\) und das \(BIC\) konzentriert.

model_performance(fit_1) |> 
  as_tibble() |> 
  select(AIC, BIC)

# A tibble: 1 × 2
    AIC   BIC
  <dbl> <dbl>
1 3076. 3093.

Gut soweit. Du kannst jetzt für jedes der Modelle das \(AIC\) berechnen und dann dir die Modelle entsprechend ordnen. Wir müssen das aber nicht tun. Wir können uns auch die Funktion compare_performance() zu nutze machen. Die Funktion gibt uns die \(R^2\)-Werte wieder wie auch die \(AIC\) sowie die \(s^2_{\epsilon}\) als sigma wieder. Wir haben also alles zusammen was wir brauchen. Darüber hinaus kann die Funktion auch die Modelle rangieren. Das nutzen wir natürlich gerne.

comp_res <- compare_performance(fit_1, fit_2, fit_3, fit_4, fit_5, rank = TRUE)

comp_res

# Comparison of Model Performance Indices

Name  | Model |    R2 | R2 (adj.) |  RMSE | Sigma | AIC weights | AICc weights
------------------------------------------------------------------------------
fit_2 |    lm | 0.739 |     0.738 | 1.926 | 1.933 |       0.644 |        0.650
fit_5 |    lm | 0.731 |     0.729 | 0.099 | 0.099 |    2.58e-09 |     2.61e-09
fit_3 |    lm | 0.739 |     0.737 | 1.926 | 1.935 |       0.237 |        0.235
fit_4 |    lm | 0.739 |     0.737 | 1.925 | 1.935 |       0.119 |        0.115
fit_1 |    lm | 0.316 |     0.314 | 3.118 | 3.126 |   5.42e-126 |    5.57e-126

Name  | BIC weights | Performance-Score
---------------------------------------
fit_2 |       0.959 |            82.69%
fit_5 |    3.84e-09 |            56.58%
fit_3 |       0.039 |            50.83%
fit_4 |       0.002 |            45.01%
fit_1 |   7.27e-125 |             0.00%

Anhand der Ausgabe der Funktion compare_performance() sehen wir, dass unser Modell fit_2 das beste Modell ist. Zwar ist die Streuung der Residuen nicht die Kleinste aller Modelle (Sigma = 1.933) aber wir haben ein hohes \(R^2_{adj}\) und auch ein kleines \(AIC\). Wir würden damit sagen, dass das Modell fit_2 mit den Variablen animal und sex für \(x\) das Outcome jump_length am besten statistisch beschreibt.

In Abbildung 45.25 sehen wir die Ausgabe der Funktion compare_performance() nochmal visualisiert. Wir können dann die einzelnen Modelle nochmal besser vergleichen. Auch siehst du hier, ob ein Modell in einem Bereich sehr schlecht ist oder aber nur in einem Bereich sehr gut.

Abbildung 45.25— Graphische Darstellung der Funktion compare_performance() Wir sehen hier die einzelnen Gütekriterien in einer Übersicht dargestellt.

Zum Ende stellt sich die Frage nach dem statistischen Test. Können wir auch statistisch Testen, ob das Modell fit_1 signifikant unterschiedlich ist? Ja wir können die Funktion test_vuong() nutzen um ein Model zu den anderen Modellen zu vergleichen. Wenn du mehrere Modell miteinander vergleichen möchtest, dann muss du die Funktion mehrfach ausführen.

test_vuong(fit_1, fit_2, fit_3, fit_4, fit_5)

Name  | Model | Omega2 | p (Omega2) |      LR | p (LR)
------------------------------------------------------
fit_1 |    lm |        |            |         |       
fit_2 |    lm |   0.41 |     < .001 |  -18.46 | < .001
fit_3 |    lm |   0.41 |     < .001 |  -18.45 | < .001
fit_4 |    lm |   0.41 |     < .001 |  -18.48 | < .001
fit_5 |    lm |   0.47 |     < .001 | -123.94 | < .001
Each model is compared to fit_1.

Nutzen wir die Modellselektion auch einmal an einem konkreten Beispiel. War die Transformation sinnvoll? Wir haben also ein Outcome \(y\) vorliegen und wollen wissen, ob wir nicht lieber das Modell mit einem \(y\) mit einer \(log\)-Transformation rechnen sollten. Wir nutzen dazu den Datensatz mit den Hunde- und Katzenflöhen und die Schlüpfdauer als Outcome. Also wie lange brauchen die Flöhe bis die Flöhe aus den Eiern geschlüpft sind.

model_tbl <- read_csv2("data/flea_dog_cat_length_weight.csv") |>
  mutate(log_hatch_time = round(log(hatch_time), 2))

Wir bauen uns jetzt zwei Modelle. Zum einen das Modell fit_raw mit der hatch_time und den Variablen ainmal und sex. Das zweite Modell enthält die log(hatch_time) also das Outcome \(y\) als \(log\)-Transformation. Wiederum sind die \(x\) variablen die gleichen wie im untransformierten Modell. Wir fitten beide Modelle und speichern die Objekte entsprechend ab.

fit_raw <- lm(hatch_time ~ animal + sex, data = model_tbl)
fit_log <- lm(log_hatch_time ~ animal + sex, data = model_tbl)

Die Funktion compare_performance() erlaubt uns wieder die beiden Modelle fit_raw und fit_log miteinander zu vergleichen. Wir erhalten die folgende Ausgabe der Funktion.

comp_res <- compare_performance(fit_raw, fit_log, rank = TRUE)

comp_res

# Comparison of Model Performance Indices

Name    | Model |    R2 | R2 (adj.) |    RMSE |   Sigma | AIC weights
---------------------------------------------------------------------
fit_log |    lm | 0.006 |     0.001 |   0.992 |   0.995 |        1.00
fit_raw |    lm | 0.008 |     0.003 | 789.441 | 792.086 |    0.00e+00

Name    | AICc weights | BIC weights | Performance-Score
--------------------------------------------------------
fit_log |         1.00 |        1.00 |            71.43%
fit_raw |     0.00e+00 |    0.00e+00 |            28.57%

In diesem Beispiel wäre das \(log\)-transformierte Modell das bessere statistische Modell. Wir müssen jetzt überlegen, ob wir den Preis dafür bezahlen wollen. Wenn wir nämlich alles mit der \(log\)-Transformation rechnen, dann erhalten wir auch alle Koeffizienten auf der \(log\)-Skala (eng. log scale). Hier kannst du nur von Fall zu Fall entscheiden.

Referenzen

Dormann, C. F. (2013). Parametrische Statistik. Springer.

Gibney, M. J. (2019). Ultra-processed foods: definitions and policy issues. Current developments in nutrition, 3(2), nzy077.

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

Kéry, M. (2010). Introduction to WinBUGS for ecologists: Bayesian approach to regression, ANOVA, mixed models and related analyses. Academic Press.

McGrosky, A., Luke, A., Arab, L., Bedu-Addo, K., Bonomi, A. G., Bovet, P., Brage, S., Buchowski, M. S., Butte, N., Camps, S. G., et al. (2025). Energy expenditure and obesity across the economic spectrum. Proceedings of the National Academy of Sciences, 122(29), e2420902122.