39 Multiple lineare Regression

Letzte Änderung am 02. April 2024 um 09:52:43

“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

In dem vorherigen Kapitel haben wir uns mit der simplen linearen Regression beschäftigt. In diesem Kapitel wollen wir uns nun mit den Grundlagen der multiplen linearen Regression beschäftigen. Sagen, dass wir ein multiples lineares Modell vorliegen haben, wenn wir mehrere Einflussvariable \(x_1\) bis \(x_p\) vorliegen haben. Bei dem Outcome \(y\), was einer Verteilung folgt, ändert sich erstmal nichts. Wir wollen jetzt herausfinden, welchen Einfluss oder Effekt die Variablen \(x_1, ..., x_p\) auf das Outcome \(y\) hat. Sehr simple gesprochen legen wir eine Gerade durch eine mehrdimensionale Punktewolke. Wir konzentrieren uns hier in diesem Kapitel auf nur normalverteilte Outcomes. Später kannst du dann aber auch andere Verteilungsfamilien modellieren.

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.

Wir erinnern und nochmal, das ein simples lineares Modell nur ein \(x_1\) hat:

\[ y \sim x_1 \]

Ein multiples lineares Modell hat von \(x_1\) bis \(x_p\) Einflussvariablen:

\[ y \sim x_1 + x_2 + ... + x_p \]

Wir brauchen viele der Konzepte aus den vorherigen Kapiteln zur Modellgüte sowie Konfidenzintervallen. Vieles fällt hier zusammen. Schau dir ruhig nochmal die vorherigen Kapitel an, wenn dir etwas unklar ist. Häufig wollen wir auch etwas spezifischer sein. Wir schreiben daher einen Faktor mit einem \(f\) wie folgt.

\[ y \sim f_1 \]

Wenn wir dann noch einen Block mit hinzunehmen, dann erweitern wir das Modell um einen Blockfaktor mit einem \(b\) wie folgt.

\[ y \sim f_1 + b_1 \]

Es kann auch sein, dass wir noch einen zufälligen Effekt (eng. random effect) mit in das Modell nehmen wollen. Wir bezeichnen einen zufälligen Effekt mit einem \(z\) wie folgt.

\[ y \sim f_1 + z_1 \]

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

\[ jump\_length \sim animal + sex + weight \]

Das hilft uns jetzt nur bedingt, denn wir wollen ja aus einer Modellierung in R die Koeffizienten der Regression wiederbekommen. Dafür müssen wir uns nochmal klar werden, was die Koeffizienten einer Regression sind. Das sind zum einen der y-Achsenabschnitt \(\beta_0\) und die Steigung der einzelnen Variablen mit \(\beta_1\) bis \(\beta_p\). Wir erhalten auch bei einem multiplen Regressionsmodell nur einen Vektor mit den Residuen wieder.

\[ y \sim \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_p x_p + \epsilon \]

Dennoch sind die Residuen \(\epsilon\) normalverteilt. Wir haben im Mittel einen Abstand der \(\epsilon\)’s zur geraden von 0 und eine Streuung von \(s^2_{\epsilon}\).

\[ \epsilon \sim \mathcal{N}(0, s^2_{\epsilon}) \]

Wir zeichen im Prinzip eine Gerade durch den \(p\)-dimensionellen Raum.

In diesem Kapitel wollen wir uns nun folgende Aspekte einmal genauer anschauen. Wir gehen hier nicht so sehr auf den theoretischen Hintergrund ein, wollen aber verstehen wie wir die Statistik nutzen und anwenden können.

Wir besprechen die Modellmatrix in Kapitel 39.2. Wir wollen uns Fragen, wie R eigentlich zu den Koeffizienten der Regressionsgleichung kommt.
Wir besprechen die Interpretation der Koeffizienten der linearen Regression in Kapitel 39.3.
Wir fragen uns, was sind eigentlich Confounder (deu. Verwechsler) im Kapitel 39.4.
Wir fragen uns, was ist eigentlich ein Variance Inflation Factor (VIF) für kontinuierliche \(x\) in Kapitel 39.5.
Wir besprechen, wie wir verschiedene Modelle miteinander vergleichen können in Kapitel 39.6.

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.

39.1 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.

39.2 Die Modellmatrix in R

Als erstes wollen wir verstehen, wie ein Modell in R aussieht. Dann können wir auch besser verstehen, wie die eigentlichen Koeffizienten aus dem Modell entstehen. Wozu brauche ich das? Eine gute Frage. Du brauchst das Verständnis der Modellmatrix, wenn du verstehen willst wie R in einer linearen Regression zu den Koeffizienten kommt. Wir bauen uns dafür ein sehr simplen Datensatz. Wir bauen uns einen Datensatz mit Schlangen.

snake_tbl <- tibble(
  svl = c(40, 45, 39, 51, 52, 57, 58, 49),
  mass = c(6, 8, 5, 7, 9, 11, 12, 10),
  region = as_factor(c("west","west", "west", "nord","nord","nord","nord","nord")),
  color = as_factor(c("schwarz", "schwarz", "rot", "rot", "rot", "blau", "blau", "blau"))
)

In der Tabelle 39.1 ist der Datensatz snake_tbl nochmal dargestellt. Wir haben die Schlangenlänge svl als Outcome \(y\) sowie das Gewicht der Schlangen mass, die Sammelregion region und die Farbe der Schlangen color. Dabei ist mass eine kontinuierliche Variable, region eine kategorielle Variable als Faktor mit zwei Leveln und color eine kategorielle Variable als Faktor mit drei Leveln.

Tabelle 39.1— Datensatz zu Schlangen ist entlehnt und modifiiert nach Kéry (2010, p. 77)

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

Wir wollen uns nun einmal anschauen, wie ein Modell in R sich zusammensetzt. Je nachdem welche Spalte \(x\) wir verwenden um den Zusammenhang zum \(y\) aufzuzeigen.

39.2.1 Kontinuierliches \(x\)

Im ersten Schritt wollen wir uns einmal das Modell mit einem kontinuierlichen \(x\) anschauen. Daher bauen wir uns ein lineares Modell mit der Variable mass. Wir erinnern uns, dass mass eine kontinuierliche Variable ist, da wir hier nur Zahlen in der Spalte finden. Die Funktion model.matrix() gibt uns die Modellmatrix wieder.

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

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

In der ersten Spalte ist der Intercept angegeben, danach folgt dann die Spalte mass als kontinuierliche Variable.

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 \\ 1 & 8 \\ 1 & 5 \\ 1 & 7 \\ 1 & 9 \\ 1 & 11\\ 1 & 12\\ 1 & 10\\ \end{pmatrix} \times \begin{pmatrix} \beta_0 \\ \beta_{mass} \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 mit der Funktion lm() fitten. Wir nutzen dann die Funktion coef() um uns die Koeffizienten aus dem Objekt fit_1 wiedergeben zu lassen.

fit_1 <- lm(svl ~ mass, data = snake_tbl) 
fit_1 |> coef() |> round(2)

(Intercept)        mass 
      26.71        2.61

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

fit_1 |> residuals() |> round(2)

    1     2     3     4     5     6     7     8 
-2.36 -2.57 -0.75  6.04  1.82  1.61  0.00 -3.79

Wir können jetzt die Koeffizienten in die Modellmatrix ergänzen. Wir haben den Intercept mit \(\beta_0 = 26.71\) geschätzt. Weiter ergänzen wir die Koeffizienten aus dem linearen Modell für mass mit \(\beta_{mass}=2.61\). Ebenfalls setzen wir Werte für die Residuen für jede der Beobachtungen in die Gleichung ein. Wir erhalten dann folgende ausgefüllte Gleichung mit den Matrixen.

\[ \begin{pmatrix} 40 \\ 45 \\ 39 \\ 50 \\ 52 \\ 57 \\ 58 \\ 59 \\ \end{pmatrix} = \begin{pmatrix} 26.71 & \phantom{0}6 \cdot 2.61\\ 26.71 & \phantom{0}8 \cdot 2.61\\ 26.71 & \phantom{0}5 \cdot 2.61\\ 26.71 & \phantom{0}7 \cdot 2.61\\ 26.71 & \phantom{0}9 \cdot 2.61\\ 26.71 & 11\cdot 2.61\\ 26.71 & 12\cdot 2.61\\ 26.71 & 10\cdot 2.61\\ \end{pmatrix} + \begin{pmatrix} -2.36\\ -2.57\\ -0.75 \\ +6.04\\ +1.82\\ +1.61\\ \phantom{+}0.00\\ -3.79\\ \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. Wie du siehst ergänzen wir hier noch eine Reihe von \(+\) um den Intercept mit der Steigung zu verbinden. Steht ja auch so in der Gleichung des linearen Modells drin, alles wird mit einem \(+\) miteinander verbunden.

c(26.71 +  6*2.61 - 2.36,
  26.71 +  8*2.61 - 2.57,
  26.71 +  5*2.61 - 0.75,
  26.71 +  7*2.61 + 6.04,
  26.71 +  9*2.61 + 1.82,
  26.71 + 11*2.61 + 1.61,
  26.71 + 12*2.61 + 0.00,
  26.71 + 10*2.61 - 3.79) |> round()

[1] 40 45 39 51 52 57 58 49

Oh ha! Die Zahlen, die wir rauskriegen, sind die gleichen Werte die unser Outcome \(y\) hat. Das heißt, die ganze Sache hat funktioniert.

39.2.2 Kategorielles \(x\) mit 2 Leveln

Im diesem Schritt wollen wir uns einmal das Modell mit einem kategoriellen \(x\) mit 2 Leveln anschauen. Wir bauen uns ein Modell mit der Variable region``. Die Funktionmodel.matrix()` gibt uns die Modelmatrix wieder.

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

# A tibble: 8 x 2
  `(Intercept)` regionnord
          <dbl>      <dbl>
1             1          0
2             1          0
3             1          0
4             1          1
5             1          1
6             1          1
7             1          1
8             1          1

In der ersten Spalte ist der Intercept angegeben, danach folgt dann die Spalte regionnord. In dieser Spalte 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.

\[ \begin{pmatrix} 40 \\ 45 \\ 39 \\ 50 \\ 52 \\ 57 \\ 58 \\ 59 \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 1 & 0 \\ 1 & 0\\ 1 & 1\\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ \end{pmatrix} \times \begin{pmatrix} \beta_0 \\ \beta^{region}_{nord} \\ \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_2 <- lm(svl ~ region, data = snake_tbl) 
fit_2 |> coef() |> round(2)

(Intercept)  regionnord 
      41.33       12.07

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

fit_2 |> residuals() |> round(2)

    1     2     3     4     5     6     7     8 
-1.33  3.67 -2.33 -2.40 -1.40  3.60  4.60 -4.40

Wir können jetzt die Koeffizienten ergänzen mit \(\beta_0 = 41.33\) für den Intercept. Weiter ergänzen wir die Koeffizienten für die Region und das Level nord mit \(\beta^{region}_{nord} = 12.07\). 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} 41.33 & 0 \cdot 12.07 \\ 41.33 & 0 \cdot 12.07 \\ 41.33 & 0 \cdot 12.07 \\ 41.33 & 1 \cdot 12.07 \\ 41.33 & 1 \cdot 12.07 \\ 41.33 & 1 \cdot 12.07 \\ 41.33 & 1 \cdot 12.07 \\ 41.33 & 1 \cdot 12.07 \\ \end{pmatrix} + \begin{pmatrix} -1.33\\ +3.67 \\ -2.33 \\ -2.40 \\ -1.40 \\ +3.60 \\ +4.60 \\ -4.40 \\ \end{pmatrix} \]

c(41.33 + 0*12.07 - 1.33,
  41.33 + 0*12.07 + 3.67,
  41.33 + 0*12.07 - 2.33,
  41.33 + 1*12.07 - 2.40,
  41.33 + 1*12.07 - 1.40,
  41.33 + 1*12.07 + 3.60,
  41.33 + 1*12.07 + 4.60,
  41.33 + 1*12.07 - 4.40) |> round()

[1] 40 45 39 51 52 57 58 49

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

39.2.3 Kategorielles \(x\) mit >2 Leveln

Im diesem Schritt wollen wir uns einmal das Modell mit einem kategoriellen \(x\) mit >2 Leveln anschauen. Wir bauen uns ein Modell mit der Variable color. Die Funktion model.matrix() gibt uns die Modelmatrix wieder.

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

# A tibble: 8 x 3
  `(Intercept)` colorrot colorblau
          <dbl>    <dbl>     <dbl>
1             1        0         0
2             1        0         0
3             1        1         0
4             1        1         0
5             1        1         0
6             1        0         1
7             1        0         1
8             1        0         1

In der ersten Spalte ist der Intercept angegeben, danach folgt dann die Spalten für color. 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.

\[ \begin{pmatrix} 40 \\ 45 \\ 39 \\ 50 \\ 52 \\ 57 \\ 58 \\ 59 \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 0 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 1 & 0\\ 1 & 0 & 1\\ 1 & 0 & 1\\ 1 & 0 & 1\\ \end{pmatrix} \times \begin{pmatrix} \beta_0 \\ \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_3 <- lm(svl ~ color, data = snake_tbl) 
fit_3 |> coef() |> round(2)

(Intercept)    colorrot   colorblau 
      42.50        4.83       12.17

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

fit_3 |> residuals() |> round(2)

    1     2     3     4     5     6     7     8 
-2.50  2.50 -8.33  3.67  4.67  2.33  3.33 -5.67

Wir können jetzt die Koeffizienten ergänzen mit \(\beta_0 = 25\) für den Intercept. Weiter ergänzen wir die Koeffizienten für die Farbe und das Level rot mit \(\beta^{color}_{rot} = 4.83\) und für die Farbe und das Level blau mit \(\beta^{color}_{blau} = 12.17\). 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} 42.50 & 0 \cdot 4.83& 0 \cdot 12.17 \\ 42.50 & 0 \cdot 4.83& 0 \cdot 12.17\\ 42.50 & 1 \cdot 4.83& 0 \cdot 12.17\\ 42.50 & 1 \cdot 4.83& 0 \cdot 12.17\\ 42.50 & 1 \cdot 4.83& 0 \cdot 12.17\\ 42.50 & 0 \cdot 4.83& 1 \cdot 12.17\\ 42.50 & 0 \cdot 4.83& 1 \cdot 12.17\\ 42.50 & 0 \cdot 4.83& 1 \cdot 12.17\\ \end{pmatrix} + \begin{pmatrix} -2.50 \\ +2.50 \\ -8.33 \\ +3.67 \\ +4.67 \\ +2.33 \\ +3.33 \\ -5.67 \\ \end{pmatrix} \]

c(42.50 + 0*4.83 + 0*-12.17 - 2.50,
  42.50 + 0*4.83 + 0*-12.17 + 2.50,
  42.50 + 1*4.83 + 0*-12.17 - 8.33,
  42.50 + 1*4.83 + 0*-12.17 + 3.67,
  42.50 + 1*4.83 + 0*-12.17 + 4.67,
  42.50 + 0*4.83 + 1*-12.17 + 2.33,
  42.50 + 0*4.83 + 1*-12.17 + 3.33,
  42.50 + 0*4.83 + 1*-12.17 - 5.67) |> round()

[1] 40 45 39 51 52 33 34 25

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

39.2.4 Das volle Modell

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 x 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.

\[ \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} \]

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.

39.3 Interpretation von \(x\)

Wi interpretiee wir nun das Ergebnis einer linearen Regression? Zum einen nutzen wir häufig das Modell nur um das Modell dann weiter in einem multiplen Vergleich zu nutzen. Hier wollen wir uns jetzt aber wirklich die Koeffizienten aus einer multiplen linearen Regression anschauen udn diese Zahlen einmal interpretieren. Wichtig ist, dass wir ein normalverteiltes \(y\) vorliegen haben und uns verschiedene Formen des \(x\) anschauen. Wir betrachten ein kontinuierliches \(x\), ein kategorielles \(x\) mit zwei Leveln und ein kategorielles \(x\) mit drei oder mehr Leveln.

39.3.1 Kontinuierliches \(x\)

Bauen wir uns also einmal einen Datensatz mit einem kontinuierlichen \(x\) und einem normalverteilten \(y\). Unser \(x\) soll von 1 bis 7 laufen. Wir erschaffen uns das \(y\) indem wir das \(x\) mit 1.5 multiplizieren, den \(y\)-Achsenabschnitt von 5 addieren und einen zufälligen Fehler aus einer Normalverteilung mit \(\mathcal{N}(0, 1)\) aufaddieren. Wir haben also eine klassische Regressionsgleichung mit \(y = 5 + 1.5 \cdot x\).

set.seed(20137937)
cont_tbl <- tibble(x = seq(from = 1, to = 7, by = 1),
                   y = 5 + 1.5 * x + rnorm(length(x), 0, 1))
cont_tbl

# A tibble: 7 x 2
      x     y
  <dbl> <dbl>
1     1  5.93
2     2  8.06
3     3  8.95
4     4 12.1 
5     5 11.9 
6     6 14.9 
7     7 16.3

Wenn wir keinen Fehler addieren würden, dann hätten wir auch eine Linie von Punkten wie an einer Perlschnur aufgereiht. Dann wäre das \(R^2 = 1\) und wir hätten keine Varianz in den Daten. Das ist aber in einem biologischen Setting nicht realistisch, dass unser \(y\) vollständig von \(x\) erklärt wird.

Schauen wir uns nochmal die Modellmatrix an. Hier erwartet uns aber keine Überraschung. Wir schätzen den Intercept und dann kommt der Wert für jedes \(x\) in der zweiten Spalte.

model.matrix(y ~ x, data = cont_tbl)

  (Intercept) x
1           1 1
2           1 2
3           1 3
4           1 4
5           1 5
6           1 6
7           1 7
attr(,"assign")
[1] 0 1

Im Folgenden schätzen wir jetzt das lineare Modell um die Koeffizienten der Geraden zu erhalten. Wir nutzen die Funktion model_parameter() aus dem R Paket {parameters} für die Ausgabe der Koeffizienten.

lm(y ~ x, data = cont_tbl) |> 
  model_parameters() |> 
  select(Parameter, Coefficient)

# Fixed Effects

Parameter   | Coefficient
-------------------------
(Intercept) |        4.32
x           |        1.71

Wir erhalten einen Intercept von \(4.32\) und eine Steigung von \(x\) mit \(1.71\). Wichtig nochmal, wir haben uns hier zufällige Zahlen erstellen lassen. Wenn du oben den Fehler rausnimmst, dann erhälst du auch die exakten Zahlen für den Intercept und die Steigung von \(x\) wieder. Wir sehen also, wenn wir ein kontinuierliches \(x\) haben, dann können wir das \(x\) in dem Sinne einr Steigung interpretieren. Steigt das \(x\) um 1 Einheit an, so erhöht sich das \(y\) um den Wert der Steigung von \(x\).

In Abbildung 39.1 sehen wir den Zusammenhang nochmal graphisch dargestellt. Wir sehen, dass wir die voreingestellten Parameter von \(\beta_0 = 5\) und \(\beta_1 = 1.5\) fast treffen.

Abbildung 39.1— Graphische Darstellung der Interpretation von einem kontinuierlichen \(x\).

39.3.2 Kategorielles \(x\) mit 2 Leveln

Etwas anders wird der Fall wenn wir ein kategorielles \(x\) mit 2 Leveln vorliegen haben. Wir bauen faktisch zwei Punktetürme an zwei \(x\) Positionen auf. Dennoch können wir durch diese Punkte eine Gerade zeichnen. Bauen wir uns erst die Daten mit der Funktion rnorm(). Wir haben zwei Gruppen vorliegen, die Gruppe A hat sieben Beobachtungen und einen Mittelwert von 10. Die Gruppe B hat ebenfalls sieben Beobchatungen und einen Mittelwert von 15. Der Effekt zwischen den beiden Gruppen A und B ist die Mittelwertsdifferenz \(\Delta_{A-B}\) ist somit 5. Wir erhlalten dann folgenden Datensatz wobei die Werte der Gruppe A um die 10 streuen und die Werte der Gruppe B um die 15 streuen.

set.seed(20339537)
cat_two_tbl <- tibble(A = rnorm(n = 7, mean = 10, sd = 1),
                      B = rnorm(n = 7, mean = 15, sd = 1)) |> 
  gather(key = x, value = y) |> 
  mutate(x = as_factor(x))
cat_two_tbl

# A tibble: 14 x 2
   x         y
   <fct> <dbl>
 1 A     10.0 
 2 A     10.8 
 3 A     10.7 
 4 A     11.0 
 5 A      9.25
 6 A      8.98
 7 A      9.71
 8 B     15.1 
 9 B     16.0 
10 B     14.5 
11 B     15.1 
12 B     14.2 
13 B     16.8 
14 B     15.6

Wir wollen uns wieder die Modellmatrix einmal anschauen. Wir sehen hier schon einen Unterschied. Zum einen sehen wir, dass der Intercept für alle Beobachtungen geschätzt wird und in der zweiten Spalte nur xB steht. Somit werden in der zweiten Spalte nur die Beobachtungen in der Gruppe B berücksichtigt.

model.matrix(y ~ x, data = cat_two_tbl)

   (Intercept) xB
1            1  0
2            1  0
3            1  0
4            1  0
5            1  0
6            1  0
7            1  0
8            1  1
9            1  1
10           1  1
11           1  1
12           1  1
13           1  1
14           1  1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"

Wir rechnen einmal das lineare Modell und lassen uns überraschen was als Ergebnis herauskommt. Wir nutzen die Funktion model_parameter() aus dem R Paket {parameters} für die Ausgabe der Koeffizienten.

lm(y ~ x, data = cat_two_tbl) |> 
  model_parameters() |> 
  select(Parameter, Coefficient)

# Fixed Effects

Parameter   | Coefficient
-------------------------
(Intercept) |       10.07
xB          |        5.27

Nun erhalten wir als den Intercept 10.07 und die Steigung xB mit 5.27 zurück. Wir sehen, der Intercept ist der Mittelwert der Gruppe A und die das xB ist die Änderung von dem Mittelwert A zu dem Mittelwert B. Wir erhalten die Mittelwertsdifferenz \(\Delta_{A-B}\) von 5.27 zurück. In Tabelle 39.2 siehst du den Zusammenhang von den Faktorleveln A und B, den jeweiligen Mittelwerte der Level sowie die Differenz zum Mittel von dem ersten Level A.

Tabelle 39.2— Zusammenhang von den Mittelwerten der Level des Faktores \(x\) und deren Differenz zu Level A.

Factor x	Mean of level	Difference to level A
A	10.07	0.00
B	15.33	5.27

In Abbildung 39.2 kannst du nochmal den visuellen Zusammenhang zwischen den einzelnen Beobachtungen und der sich ergebenen Geraden sehen. Die Gerade geht durch die beiden Mittelwerte der Gruppe A und B. Daher ist die Steigung der Mittlwertsunterschied zwischen der Gruppe A und der Gruppe B. Steigt das \(x\) um 1 Einheit an, also springt von Gruppe A zu Gruppe B, so erhöht sich das \(y\) um den Mittelwertsunterschied \(\Delta_{A-B}\).

Abbildung 39.2— Graphische Darstellung der Interpretation von einem kategoriellen \(x\) mit 2 Leveln.

39.3.3 Kategorielles \(x\) mit >2 Leveln

Nachdem wir das Problem der Interpretation von einem kategoriellen \(x\) mit zwei Leveln verstanden haben, werden wir uns jetzt den Fall für ein kategoriellen \(x\) mit drei oder mehr Leveln anschauen. Wir nutzen hier nur drei Level, da es für vier oder fünf Level gleich abläuft. Du kannst das Datenbeispiel gerne auch noch um eine vierte Gruppe erweitern und schauen was sich da ändert.

Unser Datensatz besteht aus drei Gruppen A, B und C mit den Mittelwerten von 10, 15 und 3. Hierbei ist wichtig, dass wir in der Gruppe C nur einen Mittelwert von 3 haben. Also wir haben keinen linearen Anstieg über die drei Gruppen.

set.seed(20339537)
cat_three_tbl <- tibble(A = rnorm(n = 7, mean = 10, sd = 1),
                        B = rnorm(n = 7, mean = 15, sd = 1),
                        C = rnorm(n = 7, mean = 3, sd = 1)) |> 
  gather(key = x, value = y) |> 
  mutate(x = as_factor(x))
cat_three_tbl

# A tibble: 21 x 2
   x         y
   <fct> <dbl>
 1 A     10.0 
 2 A     10.8 
 3 A     10.7 
 4 A     11.0 
 5 A      9.25
 6 A      8.98
 7 A      9.71
 8 B     15.1 
 9 B     16.0 
10 B     14.5 
# i 11 more rows

Schauen wir uns einmal die Modellmatrix an. Wir sehen wieder, dass der Intercept über alle Gruppen geschätzt wird und wir Koeffizienten für die Gruppen B und C erhalten. Mit der Modellmatrix wird dann auch das lineare Modell geschätzt.

model.matrix(y ~ x, data = cat_three_tbl)

   (Intercept) xB xC
1            1  0  0
2            1  0  0
3            1  0  0
4            1  0  0
5            1  0  0
6            1  0  0
7            1  0  0
8            1  1  0
9            1  1  0
10           1  1  0
11           1  1  0
12           1  1  0
13           1  1  0
14           1  1  0
15           1  0  1
16           1  0  1
17           1  0  1
18           1  0  1
19           1  0  1
20           1  0  1
21           1  0  1
attr(,"assign")
[1] 0 1 1
attr(,"contrasts")
attr(,"contrasts")$x
[1] "contr.treatment"

Wir rechnen wieder das Modell mit der Funktion lm(). Wir nutzen die Funktion model_parameter() aus dem R Paket {parameters} für die Ausgabe der Koeffizienten zu erhalten.

lm(y ~ x, data = cat_three_tbl) |> 
  model_parameters() |> 
  select(Parameter, Coefficient)

# Fixed Effects

Parameter   | Coefficient
-------------------------
(Intercept) |       10.07
xB          |        5.27
xC          |       -7.41

Wir sehen wieder, dass der Intercept der Mittelwert der Gruppe A ist. Wir hatten einen Mittelwert von 10 für die Gruppe A eingestellt und wir erhalten diesen Wert wieder. Die anderen Koeffizienten sind die Änderung zum Mittelwert von der Gruppe A. In Tabelle 39.3 ist der Zusammenhang nochmal für alle Level des Faktors \(x\) dargestellt. Wir haben für die Gruppe B einen Mittelwert von 15 und für die Gruppe C einen Mittelwert von 3 eingestellt. Wir erhalten diese Mittelwerte wieder, wenn wir die Differenzen zu dem Intercept bilden.

Tabelle 39.3— Zusammenhang von den Mittelwerten der Level des Faktores \(x\) und deren Differenz zu Level A.

Factor x	Mean of level	Difference to level A
A	10.07	0.00
B	15.33	5.27
C	2.66	-7.41

Wir sehen den Zusammenhang nochmal in der Abbildung 39.3 visualisiert. Wir sehen die Änderung on der Gruppe A zu der Gruppe B sowie die Änderung von der Gruppe A zu der Gruppe C. Die Abbildung ist etwas verwirrend da wir nicht das \(x\) um 2 Einheiten erhöhen um auf den Mittelwert von C zu kommen. Wir rechnen sozusagen ausgehend von A die Änderung zu B und C aus. Dabei nehmen wir jedesmal an, dass die Gruppe B und die Gruppe C nur eine Einheit von \(x\) von A entfernt ist.

Warning in geom_segment(aes(x = 1, y = 10.1, xend = 2, yend = 10.1 + 5.27), : All aesthetics have length 1, but the data has 21 rows.
i Did you mean to use `annotate()`?

Warning in geom_segment(aes(x = 1, y = 10.1, xend = 3, yend = 10.1 - 7.41), : All aesthetics have length 1, but the data has 21 rows.
i Did you mean to use `annotate()`?

Abbildung 39.3— Graphische Darstellung der Interpretation von einem kategoriellen \(x\) mit 3 Leveln.

39.4 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 39.4 ist der Datensatz model_tbl nochmal dargestellt.

Tabelle 39.4— 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.

Das Modell jump_length ~ weight. Wir modellieren die Abhängigkeit von der Sprungweite von dem Gewicht der Flöhe über alle anderen Variablen hinweg.
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.
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.

In Abbildung 39.4 (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 39.4 (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 39.4 (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.

39.5 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.039437    1.030421    1.095417    1.140199

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 39.5— Graphische Darstellung des VIF mit der Funktion `check_model()`.

In Abbildung 39.5 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.

39.6 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 x 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 | BIC weights | Performance-Score
----------------------------------------------------------------------------------------------------------------
fit_2 |    lm | 0.739 |     0.738 | 1.926 | 1.933 |       0.644 |        0.650 |       0.959 |            82.69%
fit_5 |    lm | 0.731 |     0.729 | 0.099 | 0.099 |    2.58e-09 |     2.61e-09 |    3.84e-09 |            56.58%
fit_3 |    lm | 0.739 |     0.737 | 1.926 | 1.935 |       0.237 |        0.235 |       0.039 |            50.83%
fit_4 |    lm | 0.739 |     0.737 | 1.925 | 1.935 |       0.119 |        0.115 |       0.002 |            45.01%
fit_1 |    lm | 0.316 |     0.314 | 3.118 | 3.126 |   5.42e-126 |    5.57e-126 |   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 39.6 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 39.6— 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 | AICc weights | BIC weights | Performance-Score
----------------------------------------------------------------------------------------------------------------------
fit_log |    lm | 0.006 |     0.001 |   0.992 |   0.995 |        1.00 |         1.00 |        1.00 |            71.43%
fit_raw |    lm | 0.008 |     0.003 | 789.441 | 792.086 |    0.00e+00 |     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

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