41 Simple lineare Regression

Letzte Änderung am 03. September 2024 um 11:27:12

“Farewell to the Fairground; These rides aren’t working anymore.” — White Lies, Farewell to the Fairground

In diesem Kapitel wollen wir uns mit den Grundlagen der simplen linearen Regression beschäftigen. Damit meine ich erstmal die Idee eine Gerade durch eine Punktewolke zu zeichnen. Wir haben auf der $x$-Achse kontinuierliche Werte wie auch auf der $y$-Achse vorliegen. Das ist erstmal die simpelste Anwendung. Wir lernen hier einmal die Grundbegriffe und erweitern diese dann auf komplexere Modelle. Das heißt, wir haben ein Outcome $y$, was eine normalverteilte und kontinuierliche Variable ist, sowie eine Einflussvariable $x_1$, die eine kontinuierliche Variable ist. Wir wollen jetzt herausfinden, welchen Einfluss oder Effekt die Einflussvariable $x_1$ auf das normalverteilte Outcome $y$ hat. Wir schreiben hier $x_1$, da wir damit ausdrücken, dass wir später auch noch mehr Einflussvariablen $x$ in der multiplen linearen Regression nutzen werden.

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

Sehr simpel gesprochen legen wir also eine Gerade durch eine Punktewolke. Du kannst dich aber davon gedanklich lösen, dass die lineare Regression nur eine Methode ist um eine Gerade durch eine Punktewolke zu legen. Die lineare Regression und damit auch das statistische Modellieren kann viel mehr. Wir fassen also zusammen. Ein simples lineares Modell hat ein kontinuierliches $x$ als Einflussvariable. Später können wir dann auch andere $x$ nehmen, wie zum Beispiel einen Faktor $f$ als kategoriale Einflussvariable.

\[ y \sim x_1 \]

Damit beginnen wir hier einmal mit den Konzept der linearen Regression in dem wir wirklich sehr einfach uns eine Punktewolke vorstellen und durch diese Punktewolke die optimale Gerade zeichen wollen. Dann wollen wir wissen wie gut das geklappt hat. Liegen die Punkte weit von der Geraden entfernt oder eher näher dran? In der Abbildung 41.1 siehst du einmal eine Punktewolke die sich aus den Werten der $x$-Achse und den Werten auf der $y$-Achse ergibt. Wir wollen jetzt mathematisch verstehen, wie wir die rote Gerade entsteht. Das ist die Idee der simplen linearen Regression.

Abbildung 41.1— Scatterplot von fünfundzwanzig Beobachtungen die sich aus den paarweisen Werten für $x_1$ und $y$ ergeben. Das Ziel einer linearen Regression ist es die rote Gerade zu bestimmen.

Aber in diesem Kapitel konzentrieren wir uns erstmal darauf die Begriffe der simplen linearen Regression und deren Durchführung in R zu verstehen. Dann wollen wir noch die Ausgabe einer linearen Regression in R sowie den Unterschied zwischen einem kausalen und einen prädiktiven Modell ergründen.

Weitere Tutorien für die simple 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.
Du kannst auf der Seite Manual linear regression analysis using R nochmal weiter über die lineare Regression lesen. Der Blogpost ist sehr umfangreich und erklärt nochmal schrittweise, wie die lineare Regression in R per hand funktioniert.
Simple Regression auf Wikipedia mit weit mehr Informationen zu den Formeln und den Zusammenhängen. Ein toller zusammenfassender Artikel.

41.1 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

R Code [zeigen / verbergen]

pacman::p_load(tidyverse, magrittr, broom,
               readxl, conflicted)
conflicts_prefer(magrittr::set_names)

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

41.2 Daten

Wir wollen uns erstmal mit einem einfachen Datenbeispiel beschäftigen. Wir können die lineare Regression auf sehr großen Datensätzen anwenden, wie auch auf sehr kleinen Datensätzen. Prinzipiell ist das Vorgehen gleich. Wir nutzen jetzt aber erstmal einen kleinen Datensatz mit $n=7$ Beobachtungen. In der Tabelle 41.1 ist der Datensatz simplel_tbl dargestellt. Wir wollen den Zusammenhang zwischen der Sprungweite in [cm] und dem Gewicht in [mg] für sieben Beobachtungen modellieren.

R Code [zeigen / verbergen]

simple_tbl <- tibble(jump_length = c(1.2, 1.8, 1.3, 1.7, 2.6, 1.8, 2.7),
                     weight = c(0.8, 1, 1.2, 1.9, 2, 2.7, 2.8))

Tabelle 41.1— Datensatz mit einer normalverteilten Variable jump_length und der normalverteilten Variable weight.

jump_length	weight
1.2	0.8
1.8	1.0
1.3	1.2
1.7	1.9
2.6	2.0
1.8	2.7
2.7	2.8

In Abbildung 41.2 sehen wir die Visualisierung der Daten simple_tbl in einem Scatterplot mit einer geschätzen Gerade.

Abbildung 41.2— Scatterplot der Beobachtungen der Sprungweite in [cm] und dem Gewicht in [mg]. Die Gerade verläuft mittig durch die Punkte.

Wir schauen uns in diesem Kapitel nur eine simple lineare Regression mit einem $x_1$ an. In unserem Fall ist das $x_1$ gleich dem weight. Später schauen wir dann multiple lineare Regressionen mit mehreren $x_1,..., x_p$ an.

Bevor wir mit dem Modellieren anfangen können, müssen wir verstehen, wie ein simples Modell theoretisch aufgebaut ist. Danach können wir uns das lineare Modell in R anschauen.

41.3 Simple lineare Regression theoretisch

Wir haben nun die ersten sieben Beobachtungen in dem Objekt simple_tbl vorliegen. Wie sieht nun theoretisch eine lineare Regression aus? Wir wollen eine Grade durch Punkte legen, wie wie wir es in Abbildung 41.2 sehen. Die blaue Gerade wir durch eine Geradengleichung beschreiben. Du kenst vermutlich noch die Form $y = mx + b$. In der Statistik beschreiben wir eine solche Gerade aber wie folgt.

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

mit

$\beta_0$ als den y-Achsenabschnitt.
$\beta_1$ als der Steigung der Geraden.
$\epsilon$ als Residuen oder die Abweichungen von den $y$-Werten auf Geraden zu den einzelnen $y$-Werten der Beobachtungen.

In Tabelle 41.2 siehst du nochmal in einer Tabelle den Vergleich von der Schreibweise der linearen Regression in der Schule und in der Statistik. Darüber hinaus sind die deutschen Begriffe den englischen Begriffen gegenüber gestellt. Warum schreiben wir die Gleichung in der Form? Damit wir später noch weitere $\beta_px_p$-Paare ergänzen könen und so multiple Modelle bauen können.

Tabelle 41.2— Vergleich und Übersicht der schulischen vs. statistischen Begriffe in den linearen Regression sowie die deutschen und englischen Begriffe.

$\boldsymbol{y = mx +b}$	$\boldsymbol{y \sim \beta_0 + \beta_1 x_1 + \epsilon}$	Deutsch	Englisch
$m$	$\beta_1$	Steigung	Slope
$x$	$x_1$	Einflussvariable	Risk factor
$b$	$\beta_0$	y-Achsenabschnitt	Intercept
	$\epsilon$	Residuen	Residual

In Abbildung 41.3 sehen wir die Visualisierung der Gleichung in einer Abbildung. Die Gerade läuft durch die Punktewolke und wird durch die statistischen Maßzahlen bzw. Parameter $\beta_0$, $\beta_1$ sowie den $\epsilon$ beschrieben. Wir sehen, dass das $\beta_0$ den Intercept darstellt und das $\beta_1$ die Steigung der Geraden. Wenn wir $x$ um 1 Einheit erhöhen $x+1$, dann steigt der $y$ Wert um den Wert von $\beta_1$. Die einzelnen Abweichungen der beobachteten $y$-Wert zu den $y$-Werten auf der Gerade ($\hat{y}$) werden als Residuen oder auch $\epsilon$ bezeichnet.

Abbildung 41.3— Visualisierung der linearen Regression. Wir legen eine Gerade durch eine Punktewolke. Die Gerade wird durch die statistischen Maßzahlen bzw. Parameter $\beta_0$, $\beta_1$ sowie den $\epsilon$ beschrieben.

Schauen wir uns einmal den Zusammenhang von $y$, den beobachteten Werten, und $\hat{y}$, den geschätzen Werten auf der Gerade in unserem Beispiel an. In Tabelle 41.3 sehen wir die Berechnung der einzelnen Residuen für die Gerade aus der Abbildung 41.2. Wir nehmen jedes beobachtete $y$ und ziehen den Wert von $y$ auf der Gerade, bezeichnet als $\hat{y}$, ab. Diesen Schritt machen wir für jedes Wertepaar $(y_i; \hat{y}_i)$. In R werden die $\hat{y}$ auch fitted values genannt. Die $\epsilon$ Werte werden dann residuals bezeichnet.

Tabelle 41.3— Zusammenhang zwischen den $y$, den beobachteten Werten, und $\hat{y}$, den geschätzen Werten auf der Gerade. Wir nennen den Abstand $y_i - \hat{y}_i$ auch Residuum oder Epsilon $\epsilon$.

x	y	$\boldsymbol{\hat{y}}$	Residuen ($\boldsymbol{\epsilon}$)	Wert
0.8	1.2	1.38	$\epsilon_1 = y_1 - \hat{y}_1$	$\epsilon_1 = 1.2 - 1.38 = -0.18$
1.0	1.8	1.48	$\epsilon_2 = y_2 - \hat{y}_2$	$\epsilon_2 = 1.8 - 1.48 = +0.32$
1.2	1.3	1.58	$\epsilon_3 = y_3 - \hat{y}_3$	$\epsilon_3 = 1.3 - 1.58 = -0.28$
1.9	1.7	1.94	$\epsilon_4 = y_4 - \hat{y}_4$	$\epsilon_4 = 1.7 - 1.94 = -0.24$
2.0	2.6	1.99	$\epsilon_5 = y_5 - \hat{y}_5$	$\epsilon_5 = 2.6 - 1.99 = +0.61$
2.7	1.8	2.34	$\epsilon_6 = y_6 - \hat{y}_6$	$\epsilon_6 = 1.8 - 2.34 = -0.54$
2.8	2.7	2.40	$\epsilon_7 = y_7 - \hat{y}_7$	$\epsilon_7 = 2.7 - 2.40 = +0.30$

Die Abweichungen $\epsilon$ oder auch Residuen genannt haben einen Mittelwert von $\bar{\epsilon} = 0$ und eine Varianz von $s^2_{\epsilon} = 0.17$. Wir schreiben, dass die Residuen normalverteilt sind mit $\epsilon \sim \mathcal{N}(0, s^2_{\epsilon})$. Wir zeichnen die Gerade also so durch die Punktewolke, dass die Abstände zu den Punkten, die Residuen, im Mittel null sind. Die Optimierung erreichen wir in dem wir die Varianz der Residuuen minimieren. Folglich modellieren wir die Varianz.

Simple lineare Regression händisch

Gut, das war jetzt die theoretische Abhandlung ohne eine mathematische Formel. Es geht natürlich auch mit den nackten Zahlen. In der Tabelle 41.4 siehst du einmal sieben Beobachtungen mit dem Körpergewicht als $y$ sowie der Körpergröße als $x$. Wir wollen jetzt einmal die Regressionsgleichung bestimmen. Wir sehen also unsere Werte für $\beta_0$ und $\beta_1$ aus?

Tabelle 41.4— Sieben Messungen der Körpergröße $x$ und dem zugehörigen Körpergewicht $y$.

height	weight
167	70
188	83
176	81
186	90
192	94
205	100
198	106

Damit wir einmal wissen, was wir als Lösung erhalten würden, hier einmal die lineare Regression mit der Funktion $lm()$ und die entsprechenden Werte für den (Intercept) und der Steigung.

R Code [zeigen / verbergen]

lm(weight ~ height, data = by_hand_tbl) |> 
  coef()

(Intercept)      height 
 -75.390377    0.877845

Wir suchen dann damit die folgende Regressionsgleichung mit der Körpergröße als $x$ und dem zugehörigen Körpergewicht als $y$.

\[ weight = \beta_0 + \beta_1 \cdot height \]

Da es dann immer etwas schwer ist, sich den Zusammenhang zwischen Körpergewicht und Körpergröße vorzustellen, habe ich nochmal in der Abbildung 41.4 den Scatterplot der Daten erstellt. Die rote Gerade stellt die Regressiongleichung dar. Wir erhalten ein y-Achsenabschnitt mit $\beta_0$ von $-75.39$ sowie eine Steigung mit $\beta_1$ von $0.88$ aus unseren Daten.

R Code [zeigen / verbergen]

ggplot(by_hand_tbl, aes(height, weight)) +
  theme_minimal() +
  geom_point() +
  geom_function(fun = \(x) -75.39 + 0.88 * x, color = "red")

Abbildung 41.4— Scatterplot der sieben Messungen der Körpergröße $x$ und dem zugehörigen Körpergewicht $y$ sowie der Regressionsgerade mit $y = -75.39 + 0.88 \cdot x$. Die gerade verlauf wie erwartet mittig durch die Punktewolke.

Jetzt stellt sich die Frage, wie wir händisch die Werte für den y-Achsenabschnitt mit $\beta_0$ sowie der Steigung mit $\beta_1$ berechnen. Dafür gibt es jeweils eine Formel. Hier müssen wir dann sehr viele Summen berechnen, was ich dann gleich einmal in einer Tabelle zusammenfasse.

Formel für y-Achsenabschnitt mit $\beta_0$: \[ \beta_0 = \cfrac{(\Sigma Y)(\Sigma X^2) - (\Sigma X)(\Sigma XY)}{n(\Sigma X^2) - (\Sigma X)^2} \]
Formel für Steigung mit $\beta_1$: \[ \beta_1 = \cfrac{n(\Sigma XY) - (\Sigma X)(\Sigma Y)}{n(\Sigma X^2) - (\Sigma X)^2} \]

In der Tabelle 41.5 siehst du nochmal die originalen Datenpunkte und dann die entsprechenden Werte für das Produkt von weight und height mit $XY$ und dann die jeweiligen Quadrate der beiden mit $X^2$ und $Y^2$. Wir brauchen dann aber nicht diese Werte sondern die Summen der Werte. Das Summieren lagere ich dann nochmal in eine weitere Tabelle aus.

Tabelle 41.5— Berechnungen des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$.

height	weight	$XY$	$X^2$	$Y^2$
167	70	11690	27889	4900
188	83	15604	35344	6889
176	81	14256	30976	6561
186	90	16740	34596	8100
192	94	18048	36864	8836
205	100	20500	42025	10000
198	106	20988	39204	11236

In der abschließenden Tabelle findest du dann einmal die Summen der beobachteten Werte $X$ und $Y$ sowie des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$. Damit haben wir dann alles zusammen um die Formel oben zu füllen.

Tabelle 41.6— Summe der Datenpunkte für $X$ und $Y$ sowie des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$

height $(\Sigma X)$	weight $(\Sigma Y)$	$\Sigma XY$	$\Sigma X^2$	$\Sigma Y^2$
1312	624	117826	246898	56522

Ich habe dann die ganzen Summen einmal händisch berechnet und dann in den Formeln von oben eingesetzt. Wir erhalten dann für den y-Achsenabschnitt $\beta_0$ folgenden Wert.

\[ \beta_0 = \cfrac{624 \cdot 246898 - 1312 \cdot 117826}{7\cdot 246898 - 1312^2} = -75.39038 \]

Die ganze Berechnung habe ich dann auch einmal für die Steigung $\beta_1$ ebenfalls einmal durchgeführt.

\[ b_1 = \cfrac{7\cdot 117826 - 1312\cdot624}{7\cdot246898 - 1312^2} = 0.877845 \]

Wir sehen, es kommen die gleichen Werte für den y-Achsenabschnitt $\beta_0$ und die Steigung $\beta_1$ raus. Das hat ja schonmal sehr gut geklappt. Eine andere Art die gleiche Werte effizienter zu berechnen ist die Matrixberechnung der Koeffizienten der linearen Regression. Wir könnten dann auch komplexere Modelle mit mehr als nur einem $x$ und einem $\beta_1$ berechnen. Die grundlegende Formel siehst du einmal im Folgenden dargestellt.

\[ \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix} = \mathbf{(X^T X)^{−1}(X^T Y)} \]

Wir brauchen jetzt einiges an Matrixrechnung um die jeweiligen Formelteile zu berechnen. Ich habe dir in den folgenden Tabs einmal Schritt für Schritt die einzelnen Teile berechnet. Wir immer machen wir das eigentlich nicht so richtig per Hand, sondern nutzen einen Computer. Prinzipiell wäre eine händische Lösung natürlich möglich.

R Code [zeigen / verbergen]

X <- as.matrix(c(167, 188, 176, 186, 192, 205, 198))
X <- cbind(rep(1, 7), X) 
X

     [,1] [,2]
[1,]    1  167
[2,]    1  188
[3,]    1  176
[4,]    1  186
[5,]    1  192
[6,]    1  205
[7,]    1  198

R Code [zeigen / verbergen]

Y <- as.matrix(c(70, 83, 81, 90, 94, 100, 106))
Y

     [,1]
[1,]   70
[2,]   83
[3,]   81
[4,]   90
[5,]   94
[6,]  100
[7,]  106

R Code [zeigen / verbergen]

Xt <- t(X) 
Xt

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    1    1    1    1    1    1    1
[2,]  167  188  176  186  192  205  198

R Code [zeigen / verbergen]

XtX <- Xt %*% X
XtX

     [,1]   [,2]
[1,]    7   1312
[2,] 1312 246898

R Code [zeigen / verbergen]

XtY <- Xt %*% Y
XtY

       [,1]
[1,]    624
[2,] 117826

R Code [zeigen / verbergen]

XtXinv <- solve(XtX)
XtXinv

           [,1]         [,2]
[1,] 35.5658312 -0.188994526
[2,] -0.1889945  0.001008355

Am Ende müssen wir dann alle Teile in der Form $\mathbf{(X^T X)^{−1}(X^T Y)}$ einmal zusammenbringen. Das siehst dann in R wie folgt aus. Wir erhalten dann eine Matrix wieder wobei die erste Zeile der y-Achsenabschnitt $\beta_0$ und die zweite Zeile die Steigung $\beta_1$ ist. Wir erhalten fast die gleichen Werte wie auch schon oben.

R Code [zeigen / verbergen]

XtXinv %*% Xt %*% Y

           [,1]
[1,] -75.390377
[2,]   0.877845

Wenn du dich tiefer in die Thematik einlesen willst, dann sind hier weitere Quellen zu der Thematik unter den folgenden Links und Tutorien.

41.4 Simples lineare Regression in R

Im Allgemeinen können wir ein Modell in R wie folgt schreiben. Wir brauchen ein y auf der linken Seite und in der simplen linearen Regressione ein $x$ auf der rechten Seite der Gleichung. Wir brauchen also zwei Variablen $y$ und $x$, die natürlich nicht im Datensatz in R so heißen müssen. Im Folgenen dann einmal die Modellschreibweise für $y$ hängt ab von $x$. Das $y$ repräsentiert eine Spalte im Datensatz und das $x$ repräsentiert ebenso eine Spalte im Datensatz.

\[ \Large y \sim x \]

Konkret würden wir in unserem Beispiel das Modell wie folgt benennen. Das $y$ wird zu jump_length und das $x$ wird zu weight. Wir haben dann das Modell in der simpelsten Form definiert. Im Folgenden siehst du dann eimal die Modellschreibweise in R.

\[ \Large\overbrace{\mbox{jump\_length}}^{\Large y} \sim \overbrace{\mbox{weight}}^{\Large x} \]

Nachdem wir das Modell definiert haben, setzen wir dieses Modell jump_length ~ weight in die Funktion lm() ein um das lineare Modell zu rechnen. Wie immer müssen wir auch festlegen aus welcher Datei die Spalten genommen werden sollen. Das machen wir mit der Option data = simple_tbl. Wir speichern dann die Ausgabe der Funktion lm() in dem Objekt fit_1 damit wir die Ausgabe noch in andere Funktionen pipen können.

R Code [zeigen / verbergen]

fit_1 <- lm(jump_length ~ weight, data = simple_tbl)

Wir können jetzt mir dem Modell mehr oder minder drei Dinge tun. Abhängig von der Fragestellung liefert uns natürlich jedes der drei Möglichkeiten eine andere Antwort. Wir können auch mehrere dieser Fragen gleichzeitig beantworten, was aber meistens zu einer konfusen Analyse führt. Im Folgenden einmal eine Flowchart als grobe Übersicht deiner Möglichkeiten nach einem beispielhaften lm()-Modell als simple Regression.

flowchart TD
    A("Lineares Modell
       zum Beispiel mit lm()"):::factor --> B1
    A("Lineares Modell
       zum Beispiel mit lm()"):::factor --> B2  
    A("Lineares Modell zum Beispiel mit lm()"):::factor --> B3
    subgraph B1["Inferenzmodell"]
    B[/"anova()"\]:::factor --> E[\"emmeans()"/]:::factor
    end
    subgraph B2["Kausales Modell"]
    C("summary()"):::factor --> F("augment()"):::eval
    C("summary()"):::factor --> G("glance()"):::eval
    end   
    subgraph B3["Prädiktives Modell"]
    D("predict()"):::factor --> H("conf_mat()"):::eval
    end 
    classDef factor fill:#56B4E9,stroke:#333,stroke-width:0.75px
    classDef eval fill:#E69F00,stroke:#333,stroke-width:0.75px

Abbildung 41.5— Flowchart der Möglichkeiten nach einer Modellierung mit der Funktion lm() als eine simple lineare Regression. Es gibt natürlich noch andere Modellierungen und damit Funktion. Der generelle Ablauf bleibt jedoch gleich. Die orangen Kacheln stellen optionale Funktionen zur Güte der Regression dar.

Hier nochmal die Möglichkeiten aus der Abbildung 41.5 mit mehr Details zu der obigen Flowchart. Für mehr Informationen musst du dann die entsprechenden Kapitel besuchen, hier passt dann alles wirklich nicht hin für ein so komplexen Thema wie die Modellierung von Daten.

Wir rechnen mit dem Fit des Modells einen Gruppenvergleich oder eine Inferenzstatistik (siehe Kapitel 30 zu der ANOVA oder Kapitel 40 zu multiplen Vergleichen) – wir wollen also einen Vergleich zwischen Gruppen rechnen. Wir haben damit dann kein kontinuierliches $x$ vorliegen sondern einen Faktor mit verschiedenen Leveln als Gruppen.
Wir rechnen ein kausales Modell, uns interessieren die Effekte (siehe Kapitel 41.4.2) – wir sind an den Effekten von verschiedenen Einflussvariablen interessiert. Dabei geht es dann weniger um eine Behandlungsgruppe sondern um viele verschiedene Einflussvariablen, die verstanden werden wollen in ihrem Einfluss auf das Outcome $y$. Wir können uns dann die Ergebnisse durch die Funktionen augment() und glance() aus dem R Paket {broom} näher anschauen.
Wir rechnen ein prädiktives Modell, uns interessiert der Wert neuer Werte (siehe Kapitel 41.4.3) – wir wollen dann also eine Vorhersage rechnen. Wir bauen uns ein komplexes Modell und wollen mit diesem Modell zukünftige Werte vom Outcome $y$ vorhersagen. Wenn wir also neue Einflussvariablen $x$ messen, können wir dann anhand des Modells unbekannte Outcomes $y$ bestimmen. Wir können uns dann mehr Informationen zu der Güte der Vorhersage mit der Funktion conf_mat() aus dem R Paket {tidymodels} anschauen.

41.4.1 Inferenzmodell

In diesem Kapitel betrachten wir nicht die Analyse von Gruppenvergleichen oder der Inferenzstatistik. Das machen wir dann in eigenen Kapiteln. Bitte besuche das Kapitel 30 zu mehr Informationen zu der einfaktoriellen und zweifaktoriellen ANOVA oder Kapitel 40 zu multiplen Vergleichen. Es würde dieses Kapitel sprengen, wenn wir uns dann hier die Sachlage nochmal hier aufrollen würde.

41.4.2 Kausales Modell

Im Folgenden rechnen wir ein kausales Modell, da wir an dem Effekt des $x$ interessiert sind. Wenn also das $x_1$ um eine Einheit ansteigt, um wie viel verändert sich dann das $y$? Der Schätzer $\beta_1$ gibt uns also den Einfluss oder den kausalen Zusammenhang zwischen $y$ und $x_1$ wieder. Im ersten Schritt schauen wir uns die Ausgabe der Funktion lm() in der Funktion summary() an. Daher pipen wir das Objekt fit_1 in die Funktion summary().

R Code [zeigen / verbergen]

fit_1 |> summary

Wir erhalten folgende Ausgabe dargestellt in Abbildung 41.6.

Abbildung 41.6— Die `summary()` Ausgabe des Modells `fit_1`. Hier nicht erschrecken, wir kriegen hier sehr viele Informationen aufeinmal, die wir teilweise nicht alle benötigen. Es hängt dann eben sehr stark von der Fragestelung ab.

Was sehen wir in der Ausgabe der summary() Funktion? Als erstes werden uns die Residuen wiedergegeben. Wenn wir nur wenige Beobachtungen haben, dann werden uns die Residuen direkt wiedergegeben, sonst die Verteilung der Residuen. Mit der Funktion augment() aus dem R Paket {broom} können wir uns die Residuen wiedergeben lassen. Die Residuen schauen wir uns aber nochmal im Kapitel 43 genauer an.

R Code [zeigen / verbergen]

fit_1 |> augment()

# A tibble: 7 × 8
  jump_length weight .fitted .resid  .hat .sigma .cooksd .std.resid
        <dbl>  <dbl>   <dbl>  <dbl> <dbl>  <dbl>   <dbl>      <dbl>
1         1.2    0.8    1.38 -0.176 0.388  0.496  0.0778     -0.496
2         1.8    1      1.48  0.322 0.297  0.471  0.151       0.844
3         1.3    1.2    1.58 -0.280 0.228  0.483  0.0725     -0.701
4         1.7    1.9    1.94 -0.237 0.147  0.492  0.0275     -0.564
5         2.6    2      1.99  0.612 0.156  0.384  0.199       1.47 
6         1.8    2.7    2.34 -0.545 0.367  0.376  0.656      -1.51 
7         2.7    2.8    2.40  0.304 0.417  0.467  0.276       0.877

Im zweiten Block erhalten wir die Koeffizienten (eng. coefficients) der linearen Regression. Das heißt, wir kriegen dort $\beta_0$ als y-Achsenabschnitt sowie die Steigung $\beta_1$ für das Gewicht. Dabei ist wichtig zu wissen, dass immer als erstes der y-Achsenabschnitt (Intercept) auftaucht. Dann die Steigungen der einzelnen $x$ in dem Modell. Wir haben nur ein kontinuierliches $x$, daher ist die Interpretation der Ausgabe einfach. Wir können die Gradengleichung wie folgt formulieren.

\[ jump\_length \sim 0.97 + 0.51 \cdot weight \]

Was heißt die Gleichung nun? Wenn wir das $x$ um eine Einheit erhöhen dann verändert sich das $y$ um den Wert von $\beta_1$. Wir haben hier eine Steigung von $0.51$ vorliegen. Ohne Einheit keine Interpretation! Wir wissen, dass das Gewicht in [mg] gemessen wurde und die Sprungweite in [cm]. Damit können wir aussagen, dass wenn ein Floh 1 mg mehr wiegt der Floh 0.51 cm weiter springen würde.

Schauen wir nochmal in die saubere Ausgabe der tidy() Funktion. Wir sehen nämlich noch einen $p$-Wert für den Intercept und die Steigung von weight.

R Code [zeigen / verbergen]

fit_1 |> tidy()

# A tibble: 2 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)    0.969     0.445      2.18  0.0813
2 weight         0.510     0.232      2.20  0.0790

Wenn wir einen $p$-Wert sehen, dann brauchen wir eine Nullhypothese, die wir dann eventuell mit der Entscheidung am Signifikanzniveau $\alpha$ von 5% ablehnen können. Die Nullhypothese ist die Gleichheitshypothese. Wenn es also keinen Effekt von dem Gewicht auf die Sprungweite gebe, wie groß wäre dann $\beta_1$? Wir hätten dann keine Steigung und die Grade würde parallel zur x-Achse laufen. Das $\beta_1$ wäre dann gleich null.

\[ \begin{align*} H_0: \beta_i &= 0\\ H_A: \beta_i &\neq 0 \\ \end{align*} \]

Wir haben für jedes $\beta_i$ ein eigenes Hypothesenpaar. Meistens interessiert uns der Intercept nicht. Ob der Intercept nun durch die Null geht oder nicht ist eher von geringem Interessen.

Spannder ist aber wie sich der $p$-Wert berechnet. Der $p$-Wert basiert auf einer t-Statistik, also auf dem t-Test. Wir rechnen für jeden Koeffizienten $\beta_i$ einen t-Test. Das machen wir in dem wir den Koeffizienten estimate durch den Fehler des Koeffizienten std.error teilen.

\[ \begin{align*} T_{(Intercept)} &= \cfrac{\mbox{estimate}}{\mbox{std.error}} = \cfrac{0.969}{0.445} = 2.18\\ T_{weight} &= \cfrac{\mbox{estimate}}{\mbox{std.error}} = \cfrac{0.510}{0.232} = 2.20\\ \end{align*} \]

Wir sehen in diesem Fall, dass weder der Intercept noch die Steigung von weight signifikant ist, da die $p$-Werte mit $0.081$ und $0.079$ leicht über dem Signifikanzniveau von $\alpha$ gleich 5% liegen. Wir haben aber einen starkes Indiz gegen die Nullhypothese, da die Wahrscheinlichkeit die Daten zu beobachten sehr gering ist unter der Annahme das die Nullhypothese gilt.

Zun Abschluß noch die Funktion glance() ebenfalls aus dem R Paket {broom}, die uns erlaubt noch die Qualitätsmaße der linearen Regression zu erhalten. Wir müssen nämlich noch schauen, ob die Regression auch funktioniert hat. Die Überprüfung geht mit einem $x$ sehr einfach. Wir können uns die Grade ja anschauen. Das geht dann mit einem Model mit mehreren $x$ nicht mehr und wir brauchen andere statistische Maßzahlen.

R Code [zeigen / verbergen]

fit_1 |> glance()

# A tibble: 1 × 12
  r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
      <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
1     0.492         0.391 0.455      4.84  0.0790     1  -3.24  12.5  12.3
# ℹ 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

41.4.3 Prädiktives Modell

Neben dem kausalen Modell gibt es auch die Möglichkeit ein prädiktives Modell zu rechnen. Im Prinzip ist die Sprache hier etwas ungenau. Wir verwenden das gefittete Modell nur anders. Anstatt das Modell fit_1 in die Funktion summary() zu pipen, pipen wir die das Modell in die Funktion predict(). Die Funktion predict() kann dann für neue Daten über die Option newdata = das $y$ vorhersagen.

In unserem Fall müssen wir uns deshalb ein tibble mit einer Spalte bauen. Wir haben ja oben im Modell auch nur ein $x_1$ mit aufgenommen. Später können wir natürlich auch für multiple Modelle die Vorhersage machen. Wichtig ist, dass die Namen gleich sind. Das heißt in dem neuen Datensatz müssen die Spalten exakt so heißen wir in dem alten Datensatz in dem das Modell gefittet wurde.

R Code [zeigen / verbergen]

simple_new_tbl <- tibble(weight = c(1.7, 1.4, 2.1, 3.0)) 

predict(fit_1, newdata = simple_new_tbl) |> round(2)

   1    2    3    4 
1.84 1.68 2.04 2.50

Wie wir sehen ist die Anwendung recht einfach. Wir haben die vier jump_length Werte vorhergesagt bekommen, die sich mit dem Fit des Modells mit den neuen weight Werten ergeben. In Abbildung 41.7 sehen wir die Visualisierung der vier vorhergesagten Werte. Die Werte müssen auf der Geraden liegen.

Abbildung 41.7— Scatterplot der *alten* Beobachtungen der Sprungweite in [cm] und dem Gewicht in [mg]. Sowie der *neuen* vorhergesagten Beobachtungen auf der Geraden.

Wir werden später in der Klassifikation, der Vorhersage von $0/1$-Werten, sowie in der multiplen Regression noch andere Prädktionen und deren Maßzahlen kennen lernen. Im Rahmen der simplen Regression soll dies aber erstmal hier genügen.

```{r echo = FALSE} pacman::p_load(tidyverse, readxl, knitr, kableExtra, see) ``` # Simple lineare Regression {#sec-modeling-simple-stat} *Letzte Änderung am `r format(fs::file_info("stat-linear-reg-basic.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`* > *"Farewell to the Fairground; These rides aren't working anymore." --- White Lies, Farewell to the Fairground* In diesem Kapitel wollen wir uns mit den Grundlagen der simplen linearen Regression beschäftigen. Damit meine ich erstmal die Idee eine Gerade durch eine Punktewolke zu zeichnen. Wir haben auf der $x$-Achse kontinuierliche Werte wie auch auf der $y$-Achse vorliegen. Das ist erstmal die simpelste Anwendung. Wir lernen hier einmal die Grundbegriffe und erweitern diese dann auf komplexere Modelle. Das heißt, wir haben ein Outcome $y$, was eine normalverteilte und kontinuierliche Variable ist, sowie eine Einflussvariable $x_1$, die eine kontinuierliche Variable ist. Wir wollen jetzt herausfinden, welchen Einfluss oder Effekt die Einflussvariable $x_1$ auf das normalverteilte Outcome $y$ hat. Wir schreiben hier $x_1$, da wir damit ausdrücken, dass wir später auch noch mehr Einflussvariablen $x$ in der multiplen linearen Regression nutzen werden. ::: {layout="[15,85]" layout-valign="center"} ![](images/angel_01.png){fig-align="center" width="100%"} > Wenn du jetzt denkst 'Hä?', dann besuche doch einmal die fantastische Seite [Explained Visually \| Ordinary Least Squares Regression](https://setosa.io/ev/ordinary-least-squares-regression/) um selber zu erfahren, was eine lineare Regression macht. Auf der Seite findest du interaktive Abbildungen, die dir das Konzept der linearen Regression sehr anschaulich nachvollziehen lassen. ::: Sehr simpel gesprochen legen wir also eine Gerade durch eine Punktewolke. Du kannst dich aber davon gedanklich lösen, dass die lineare Regression *nur* eine Methode ist um eine Gerade durch eine Punktewolke zu legen. Die lineare Regression und damit auch das statistische Modellieren kann viel mehr. Wir fassen also zusammen. Ein simples lineares Modell hat *ein* kontinuierliches $x$ als Einflussvariable. Später können wir dann auch andere $x$ nehmen, wie zum Beispiel *einen* Faktor $f$ als kategoriale Einflussvariable. $$ y \sim x_1 $$ Damit beginnen wir hier einmal mit den Konzept der linearen Regression in dem wir wirklich sehr einfach uns eine Punktewolke vorstellen und durch diese Punktewolke die optimale Gerade zeichen wollen. Dann wollen wir wissen wie gut das geklappt hat. Liegen die Punkte weit von der Geraden entfernt oder eher näher dran? In der @fig-scatter-lin-00 siehst du einmal eine Punktewolke die sich aus den Werten der $x$-Achse und den Werten auf der $y$-Achse ergibt. Wir wollen jetzt mathematisch verstehen, wie wir die rote Gerade entsteht. Das ist die Idee der simplen linearen Regression. ```{r} #| echo: false #| message: false #| label: fig-scatter-lin-00 #| fig-align: center #| fig-height: 5 #| fig-width: 6 #| fig-cap: "Scatterplot von fünfundzwanzig Beobachtungen die sich aus den paarweisen Werten für $x_1$ und $y$ ergeben. Das Ziel einer linearen Regression ist es die rote Gerade zu bestimmen." tibble(x_1 = rnorm(25, 5, 1), y = 2.5 * x_1 + rnorm(25, 0, 3)) |> ggplot(aes(x_1, y)) + geom_point() + stat_smooth(method = "lm", se = FALSE, fullrange = TRUE, color = "red") + theme_minimal() ``` Aber in diesem Kapitel konzentrieren wir uns erstmal darauf die Begriffe der simplen linearen Regression und deren Durchführung in R zu verstehen. Dann wollen wir noch die Ausgabe einer linearen Regression in R sowie den Unterschied zwischen einem kausalen und einen prädiktiven Modell ergründen. ::: callout-tip ## Weitere Tutorien für die simple 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](https://setosa.io/ev/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. - Du kannst auf der Seite [Manual linear regression analysis using R](https://davetang.org/muse/2012/02/10/manual-linear-regression-analysis-using-r/) nochmal weiter über die lineare Regression lesen. Der Blogpost ist sehr umfangreich und erklärt nochmal schrittweise, wie die lineare Regression in R per hand funktioniert. - [Simple Regression auf Wikipedia](https://en.wikipedia.org/wiki/Simple_linear_regression#Fitting_the_regression_line) mit weit mehr Informationen zu den Formeln und den Zusammenhängen. Ein toller zusammenfassender Artikel. ::: ## Genutzte R Pakete Wir wollen folgende R Pakete in diesem Kapitel nutzen. ```{r echo = TRUE} #| message: false pacman::p_load(tidyverse, magrittr, broom, readxl, conflicted) conflicts_prefer(magrittr::set_names) ``` An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst. ## Daten Wir wollen uns erstmal mit einem einfachen Datenbeispiel beschäftigen. Wir können die lineare Regression auf sehr großen Datensätzen anwenden, wie auch auf sehr kleinen Datensätzen. Prinzipiell ist das Vorgehen gleich. Wir nutzen jetzt aber erstmal einen kleinen Datensatz mit $n=7$ Beobachtungen. In der @tbl-model-0 ist der Datensatz `simplel_tbl` dargestellt. Wir wollen den Zusammenhang zwischen der Sprungweite in \[cm\] und dem Gewicht in \[mg\] für sieben Beobachtungen modellieren. ```{r} #| message: false simple_tbl <- tibble(jump_length = c(1.2, 1.8, 1.3, 1.7, 2.6, 1.8, 2.7), weight = c(0.8, 1, 1.2, 1.9, 2, 2.7, 2.8)) ``` ```{r} #| message: false #| echo: false #| tbl-cap: Datensatz mit einer normalverteilten Variable `jump_length` und der normalverteilten Variable `weight`. #| label: tbl-model-0 simple_tbl <- tibble(jump_length = c(1.2, 1.8, 1.3, 1.7, 2.6, 1.8, 2.7), weight = c(0.8, 1, 1.2, 1.9, 2, 2.7, 2.8)) simple_tbl |> kable(align = "c", "pipe") ``` In @fig-scatter-lin-01 sehen wir die Visualisierung der Daten `simple_tbl` in einem Scatterplot mit einer geschätzen Gerade. ```{r} #| echo: false #| message: false #| label: fig-scatter-lin-01 #| fig-align: center #| fig-height: 5 #| fig-width: 6 #| fig-cap: "Scatterplot der Beobachtungen der Sprungweite in \\[cm\\] und dem Gewicht in \\[mg\\]. Die Gerade verläuft mittig durch die Punkte." ggplot(simple_tbl, aes(weight, jump_length)) + geom_point() + stat_smooth(method = "lm", se = FALSE, fullrange = TRUE) + theme_minimal() + xlim(0, 3.5) + ylim(0, 3.5) ``` Wir schauen uns in diesem Kapitel nur eine *simple* lineare Regression mit einem $x_1$ an. In unserem Fall ist das $x_1$ gleich dem `weight`. Später schauen wir dann *multiple* lineare Regressionen mit mehreren $x_1,..., x_p$ an. Bevor wir mit dem Modellieren anfangen können, müssen wir verstehen, wie ein *simples* Modell theoretisch aufgebaut ist. Danach können wir uns das lineare Modell in R anschauen. ## Simple lineare Regression theoretisch Wir haben nun die ersten sieben Beobachtungen in dem Objekt `simple_tbl` vorliegen. Wie sieht nun theoretisch eine lineare Regression aus? Wir wollen eine Grade durch Punkte legen, wie wie wir es in @fig-scatter-lin-01 sehen. Die blaue Gerade wir durch eine Geradengleichung beschreiben. Du kenst vermutlich noch die Form $y = mx + b$. In der Statistik beschreiben wir eine solche Gerade aber wie folgt. $$ y \sim \beta_0 + \beta_1 x_1 + \epsilon $$ mit - $\beta_0$ als den y-Achsenabschnitt. - $\beta_1$ als der Steigung der Geraden. - $\epsilon$ als Residuen oder die Abweichungen von den $y$-Werten auf Geraden zu den einzelnen $y$-Werten der Beobachtungen. In @tbl-reg-deu-eng siehst du nochmal in einer Tabelle den Vergleich von der Schreibweise der linearen Regression in der Schule und in der Statistik. Darüber hinaus sind die deutschen Begriffe den englischen Begriffen gegenüber gestellt. Warum schreiben wir die Gleichung in der Form? Damit wir später noch weitere $\beta_px_p$-Paare ergänzen könen und so *multiple* Modelle bauen können. | $\boldsymbol{y = mx +b}$ | $\boldsymbol{y \sim \beta_0 + \beta_1 x_1 + \epsilon}$ | Deutsch | Englisch | |:------------------------:|:------------------------------------------------------:|:-----------------:|:-----------:| | $m$ | $\beta_1$ | Steigung | Slope | | $x$ | $x_1$ | Einflussvariable | Risk factor | | $b$ | $\beta_0$ | y-Achsenabschnitt | Intercept | | | $\epsilon$ | Residuen | Residual | : Vergleich und Übersicht der schulischen vs. statistischen Begriffe in den linearen Regression sowie die deutschen und englischen Begriffe. {#tbl-reg-deu-eng} In @fig-lin-reg-01 sehen wir die Visualisierung der Gleichung in einer Abbildung. Die Gerade läuft durch die Punktewolke und wird durch die statistischen Maßzahlen bzw. Parameter $\beta_0$, $\beta_1$ sowie den $\epsilon$ beschrieben. Wir sehen, dass das $\beta_0$ den *Intercept* darstellt und das $\beta_1$ die Steigung der Geraden. Wenn wir $x$ um 1 Einheit erhöhen $x+1$, dann steigt der $y$ Wert um den Wert von $\beta_1$. Die einzelnen Abweichungen der beobachteten $y$-Wert zu den $y$-Werten auf der Gerade ($\hat{y}$) werden als Residuen oder auch $\epsilon$ bezeichnet. ```{r} #| echo: false #| message: false #| warning: false #| label: fig-lin-reg-01 #| fig-align: center #| fig-height: 5 #| fig-width: 6 #| fig-cap: "Visualisierung der linearen Regression. Wir legen eine Gerade durch eine Punktewolke. Die Gerade wird durch die statistischen Maßzahlen bzw. Parameter $\\beta_0$, $\\beta_1$ sowie den $\\epsilon$ beschrieben." simple_fit <- lm(jump_length ~ weight, data = simple_tbl) |> augment() ggplot(simple_tbl, aes(weight, jump_length)) + theme_minimal() + xlim(0, 3.5) + ylim(0, 3.5) + geom_vline(xintercept = 0) + geom_hline(yintercept = 0) + geom_segment(x = simple_fit$weight, y = simple_fit$jump_length, xend = simple_fit$weight, yend = simple_fit$.fitted, color = "#009E73", linewidth = 0.5, linetype = 1) + geom_segment(x = 1.035, y = 1.5, xend = 2, yend = 1.5, color = "#D55E00", linewidth = 0.5, linetype = 1) + geom_segment(x = 2, y = 1.5, xend = 2, yend = 2, color = "#D55E00", linewidth = 0.5, linetype = 1) + annotate("text", x = 0.3, y = 0.5, label = expression(beta[0]), size = 7, color = "#E69F00") + geom_curve(x = 0.2, y = 0.5, xend = 0, yend = 0.95, arrow = arrow(length = unit(0.03, "npc"), type = "closed"), curvature = -0.3, alpha = 0.3, color = "#E69F00") + annotate("text", x = 2.3, y = 1.7, label = expression(beta[1]), size = 7, color = "#E69F00") + geom_curve(x = 2.2, y = 1.7, xend = 2.02, yend = 1.75, arrow = arrow(length = unit(0.03, "npc"), type = "closed"), curvature = 0.3, alpha = 0.3, color = "#E69F00") + annotate("text", x = 1, y = 1.35, label = expression(x), size = 7, color = "#D55E00") + annotate("text", x = 2, y = 1.35, label = expression(x+1), size = 7, color = "#D55E00") + annotate("text", x = simple_fit$weight + 0.06, y = simple_fit$.fitted + simple_fit$.resid/2, label = c(expression(epsilon[1]), expression(epsilon[2]), expression(epsilon[3]), expression(epsilon[4]), expression(epsilon[5]), expression(epsilon[6]), expression(epsilon[7])), color = "#009E73") + stat_smooth(method = "lm", se = FALSE, fullrange = TRUE, color = "#56B4E9") + geom_point(size = 3) + labs(x = "Gewicht der Hundeflöhe in [mg]", y = "Sprungeweite der Hundeflöhe in [cm]") ``` Schauen wir uns einmal den Zusammenhang von $y$, den beobachteten Werten, und $\hat{y}$, den geschätzen Werten auf der Gerade in unserem Beispiel an. In @tbl-lin-reg-epsilon sehen wir die Berechnung der einzelnen Residuen für die Gerade aus der @fig-scatter-lin-01. Wir nehmen jedes beobachtete $y$ und ziehen den Wert von $y$ auf der Gerade, bezeichnet als $\hat{y}$, ab. Diesen Schritt machen wir für jedes Wertepaar $(y_i; \hat{y}_i)$. In R werden die $\hat{y}$ auch *fitted values* genannt. Die $\epsilon$ Werte werden dann *residuals* bezeichnet. | x | y | $\boldsymbol{\hat{y}}$ | Residuen ($\boldsymbol{\epsilon}$) | Wert | |:---:|:---:|:----------------------:|:----------------------------------:|:---------------------------------:| | 0.8 | 1.2 | 1.38 | $\epsilon_1 = y_1 - \hat{y}_1$ | $\epsilon_1 = 1.2 - 1.38 = -0.18$ | | 1.0 | 1.8 | 1.48 | $\epsilon_2 = y_2 - \hat{y}_2$ | $\epsilon_2 = 1.8 - 1.48 = +0.32$ | | 1.2 | 1.3 | 1.58 | $\epsilon_3 = y_3 - \hat{y}_3$ | $\epsilon_3 = 1.3 - 1.58 = -0.28$ | | 1.9 | 1.7 | 1.94 | $\epsilon_4 = y_4 - \hat{y}_4$ | $\epsilon_4 = 1.7 - 1.94 = -0.24$ | | 2.0 | 2.6 | 1.99 | $\epsilon_5 = y_5 - \hat{y}_5$ | $\epsilon_5 = 2.6 - 1.99 = +0.61$ | | 2.7 | 1.8 | 2.34 | $\epsilon_6 = y_6 - \hat{y}_6$ | $\epsilon_6 = 1.8 - 2.34 = -0.54$ | | 2.8 | 2.7 | 2.40 | $\epsilon_7 = y_7 - \hat{y}_7$ | $\epsilon_7 = 2.7 - 2.40 = +0.30$ | : Zusammenhang zwischen den $y$, den beobachteten Werten, und $\hat{y}$, den geschätzen Werten auf der Gerade. Wir nennen den Abstand $y_i - \hat{y}_i$ auch Residuum oder *Epsilon* $\epsilon$. {#tbl-lin-reg-epsilon} ```{r} #| echo: false epsilon <- c(-0.18, 0.32, -0.28, -0.24, 0.61, -0.54, 0.30) mean_e <- round(mean(epsilon), 2) var_e <- round(var(epsilon), 2) ``` Die Abweichungen $\epsilon$ oder auch Residuen genannt haben einen Mittelwert von $\bar{\epsilon} = `r mean_e`$ und eine Varianz von $s^2_{\epsilon} = `r var_e`$. Wir schreiben, dass die Residuen normalverteilt sind mit $\epsilon \sim \mathcal{N}(0, s^2_{\epsilon})$. Wir zeichnen die Gerade also so durch die Punktewolke, dass die Abstände zu den Punkten, die Residuen, im Mittel null sind. Die Optimierung erreichen wir in dem wir die Varianz der Residuuen minimieren. Folglich modellieren wir die Varianz. ::: callout-note ## Simple lineare Regression händisch Gut, das war jetzt die theoretische Abhandlung ohne eine mathematische Formel. Es geht natürlich auch mit den nackten Zahlen. In der @tbl-example-by-hand-01 siehst du einmal sieben Beobachtungen mit dem Körpergewicht als $y$ sowie der Körpergröße als $x$. Wir wollen jetzt einmal die Regressionsgleichung bestimmen. Wir sehen also unsere Werte für $\beta_0$ und $\beta_1$ aus? ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-example-by-hand-01 #| tbl-cap: "Sieben Messungen der Körpergröße $x$ und dem zugehörigen Körpergewicht $y$." by_hand_tbl <- tibble(height = c(167, 188, 176, 186, 192, 205, 198), weight = c(70, 83, 81, 90, 94, 100, 106)) by_hand_tbl |> kable(align = "c", "pipe") ``` Damit wir einmal wissen, was wir als Lösung erhalten würden, hier einmal die lineare Regression mit der Funktion $lm()$ und die entsprechenden Werte für den `(Intercept)` und der Steigung. ```{r} lm(weight ~ height, data = by_hand_tbl) |> coef() ``` Wir suchen dann damit die folgende Regressionsgleichung mit der Körpergröße als $x$ und dem zugehörigen Körpergewicht als $y$. $$ weight = \beta_0 + \beta_1 \cdot height $$ Da es dann immer etwas schwer ist, sich den Zusammenhang zwischen Körpergewicht und Körpergröße vorzustellen, habe ich nochmal in der @fig-example-by-hand-03 den Scatterplot der Daten erstellt. Die rote Gerade stellt die Regressiongleichung dar. Wir erhalten ein y-Achsenabschnitt mit $\beta_0$ von $-75.39$ sowie eine Steigung mit $\beta_1$ von $0.88$ aus unseren Daten. ```{r} #| echo: true #| message: false #| label: fig-example-by-hand-03 #| fig-align: center #| fig-height: 5 #| fig-width: 6 #| fig-cap: "Scatterplot der sieben Messungen der Körpergröße $x$ und dem zugehörigen Körpergewicht $y$ sowie der Regressionsgerade mit $y = -75.39 + 0.88 \\cdot x$. Die gerade verlauf wie erwartet mittig durch die Punktewolke." ggplot(by_hand_tbl, aes(height, weight)) + theme_minimal() + geom_point() + geom_function(fun = \(x) -75.39 + 0.88 * x, color = "red") ``` Jetzt stellt sich die Frage, wie wir *händisch* die Werte für den y-Achsenabschnitt mit $\beta_0$ sowie der Steigung mit $\beta_1$ berechnen. Dafür gibt es jeweils eine Formel. Hier müssen wir dann sehr viele Summen berechnen, was ich dann gleich einmal in einer Tabelle zusammenfasse. Formel für y-Achsenabschnitt mit $\beta_0$ : $$ \beta_0 = \cfrac{(\Sigma Y)(\Sigma X^2) - (\Sigma X)(\Sigma XY)}{n(\Sigma X^2) - (\Sigma X)^2} $$ Formel für Steigung mit $\beta_1$ : $$ \beta_1 = \cfrac{n(\Sigma XY) - (\Sigma X)(\Sigma Y)}{n(\Sigma X^2) - (\Sigma X)^2} $$ In der @tbl-example-by-hand-02 siehst du nochmal die originalen Datenpunkte und dann die entsprechenden Werte für das Produkt von `weight` und `height` mit $XY$ und dann die jeweiligen Quadrate der beiden mit $X^2$ und $Y^2$. Wir brauchen dann aber nicht diese Werte sondern die Summen der Werte. Das Summieren lagere ich dann nochmal in eine weitere Tabelle aus. ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-example-by-hand-02 #| tbl-cap: "Berechnungen des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$." by_hand_tbl |> mutate(a = height * weight, b = height^2, c = weight^2) |> set_names(c("height", "weight", "$XY$", "$X^2$", "$Y^2$")) |> kable(align = "c", "pipe") ``` In der abschließenden Tabelle findest du dann einmal die Summen der beobachteten Werte $X$ und $Y$ sowie des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$. Damit haben wir dann alles zusammen um die Formel oben zu füllen. ```{r} #| echo: false #| message: false #| warning: false #| label: tbl-example-by-hand-03 #| tbl-cap: "Summe der Datenpunkte für $X$ und $Y$ sowie des Produkts von $X$ und $Y$ sowie deren Quadrate mit $X^2$ und $Y^2$" tibble(height = sum(by_hand_tbl$height), weight = sum(by_hand_tbl$weight), a = sum(by_hand_tbl$height * by_hand_tbl$weight), b = sum(by_hand_tbl$height^2), c = sum(by_hand_tbl$weight^2),) |> set_names(c("height $(\\Sigma X)$", "weight $(\\Sigma Y)$", "$\\Sigma XY$", "$\\Sigma X^2$", "$\\Sigma Y^2$")) |> kable(align = "c", "pipe") ``` Ich habe dann die ganzen Summen einmal händisch berechnet und dann in den Formeln von oben eingesetzt. Wir erhalten dann für den y-Achsenabschnitt $\beta_0$ folgenden Wert. $$ \beta_0 = \cfrac{624 \cdot 246898 - 1312 \cdot 117826}{7\cdot 246898 - 1312^2} = -75.39038 $$ Die ganze Berechnung habe ich dann auch einmal für die Steigung $\beta_1$ ebenfalls einmal durchgeführt. $$ b_1 = \cfrac{7\cdot 117826 - 1312\cdot624}{7\cdot246898 - 1312^2} = 0.877845 $$ Wir sehen, es kommen die gleichen Werte für den y-Achsenabschnitt $\beta_0$ und die Steigung $\beta_1$ raus. Das hat ja schonmal sehr gut geklappt. Eine andere Art die gleiche Werte effizienter zu berechnen ist die Matrixberechnung der Koeffizienten der linearen Regression. Wir könnten dann auch komplexere Modelle mit mehr als nur einem $x$ und einem $\beta_1$ berechnen. Die grundlegende Formel siehst du einmal im Folgenden dargestellt. $$ \begin{pmatrix} \beta_0 \\ \beta_1 \end{pmatrix} = \mathbf{(X^T X)^{−1}(X^T Y)} $$ Wir brauchen jetzt einiges an Matrixrechnung um die jeweiligen Formelteile zu berechnen. Ich habe dir in den folgenden Tabs einmal Schritt für Schritt die einzelnen Teile berechnet. Wir immer machen wir das eigentlich nicht so richtig per Hand, sondern nutzen einen Computer. Prinzipiell wäre eine händische Lösung natürlich möglich. ::: panel-tabset ## $X$ ```{r} X <- as.matrix(c(167, 188, 176, 186, 192, 205, 198)) X <- cbind(rep(1, 7), X) X ``` ## $Y$ ```{r} Y <- as.matrix(c(70, 83, 81, 90, 94, 100, 106)) Y ``` ## $X^T$ ```{r} Xt <- t(X) Xt ``` ## $X^T X$ ```{r} XtX <- Xt %*% X XtX ``` ## $X^T Y$ ```{r} XtY <- Xt %*% Y XtY ``` ## $(X^T X)^{−1}$ ```{r} XtXinv <- solve(XtX) XtXinv ``` ::: Am Ende müssen wir dann alle Teile in der Form $\mathbf{(X^T X)^{−1}(X^T Y)}$ einmal zusammenbringen. Das siehst dann in R wie folgt aus. Wir erhalten dann eine Matrix wieder wobei die erste Zeile der y-Achsenabschnitt $\beta_0$ und die zweite Zeile die Steigung $\beta_1$ ist. Wir erhalten fast die gleichen Werte wie auch schon oben. ```{r} XtXinv %*% Xt %*% Y ``` Wenn du dich tiefer in die Thematik einlesen willst, dann sind hier weitere Quellen zu der Thematik unter den folgenden Links und Tutorien. - [Hands-On Machine Learning with R \| Linear Regression](https://bradleyboehmke.github.io/HOML/linear-regression.html) - [Manual linear regression analysis using R](https://davetang.org/muse/2012/02/10/manual-linear-regression-analysis-using-r/) - [Linear Regression by Hand](https://towardsdatascience.com/linear-regression-by-hand-python-and-r-79994d47f68) - [How to Perform Linear Regression by Hand](https://www.statology.org/linear-regression-by-hand/) - [Matrix Approach to Simple Linear Regression in R](https://lasanthiwatagoda.github.io/STT4830ClassRepo/RCodes/Ch5Script.html) ::: ## Simples lineare Regression in R Im Allgemeinen können wir ein Modell in R wie folgt schreiben. Wir brauchen *ein* y auf der linken Seite und in der simplen linearen Regressione ein $x$ auf der rechten Seite der Gleichung. Wir brauchen also zwei Variablen $y$ und $x$, die natürlich nicht im Datensatz in R so heißen müssen. Im Folgenen dann einmal die Modellschreibweise für $y$ hängt ab von $x$. Das $y$ repräsentiert eine Spalte im Datensatz und das $x$ repräsentiert ebenso eine Spalte im Datensatz. $$ \Large y \sim x $$ Konkret würden wir in unserem Beispiel das Modell wie folgt benennen. Das $y$ wird zu `jump_length` und das $x$ wird zu `weight`. Wir haben dann das Modell in der simpelsten Form definiert. Im Folgenden siehst du dann eimal die Modellschreibweise in R. $$ \Large\overbrace{\mbox{jump\_length}}^{\Large y} \sim \overbrace{\mbox{weight}}^{\Large x} $$ Nachdem wir das Modell definiert haben, setzen wir dieses Modell `jump_length ~ weight` in die Funktion `lm()` ein um das lineare Modell zu rechnen. Wie immer müssen wir auch festlegen aus welcher Datei die Spalten genommen werden sollen. Das machen wir mit der Option `data = simple_tbl`. Wir speichern dann die Ausgabe der Funktion `lm()` in dem Objekt `fit_1` damit wir die Ausgabe noch in andere Funktionen pipen können. ```{r} fit_1 <- lm(jump_length ~ weight, data = simple_tbl) ``` Wir können jetzt mir dem Modell mehr oder minder drei Dinge tun. Abhängig von der Fragestellung liefert uns natürlich jedes der drei Möglichkeiten eine andere Antwort. Wir können auch mehrere dieser Fragen gleichzeitig beantworten, was aber meistens zu einer konfusen Analyse führt. Im Folgenden einmal eine Flowchart als grobe Übersicht deiner Möglichkeiten nach einem beispielhaften `lm()`-Modell als simple Regression. ```{mermaid} %%| label: fig-mermaid-simple-reg-01 %%| fig-width: 6 %%| fig-cap: "Flowchart der Möglichkeiten nach einer Modellierung mit der Funktion `lm()` als eine simple lineare Regression. Es gibt natürlich noch andere Modellierungen und damit Funktion. Der generelle Ablauf bleibt jedoch gleich. Die orangen Kacheln stellen optionale Funktionen zur Güte der Regression dar." flowchart TD A("Lineares Modell zum Beispiel mit lm()"):::factor --> B1 A("Lineares Modell zum Beispiel mit lm()"):::factor --> B2 A("Lineares Modell zum Beispiel mit lm()"):::factor --> B3 subgraph B1["Inferenzmodell"] B[/"anova()"\]:::factor --> E[\"emmeans()"/]:::factor end subgraph B2["Kausales Modell"] C("summary()"):::factor --> F("augment()"):::eval C("summary()"):::factor --> G("glance()"):::eval end subgraph B3["Prädiktives Modell"] D("predict()"):::factor --> H("conf_mat()"):::eval end classDef factor fill:#56B4E9,stroke:#333,stroke-width:0.75px classDef eval fill:#E69F00,stroke:#333,stroke-width:0.75px ``` Hier nochmal die Möglichkeiten aus der @fig-mermaid-simple-reg-01 mit mehr Details zu der obigen Flowchart. Für mehr Informationen musst du dann die entsprechenden Kapitel besuchen, hier passt dann alles wirklich nicht hin für ein so komplexen Thema wie die Modellierung von Daten. - Wir rechnen mit dem Fit des Modells einen Gruppenvergleich oder eine Inferenzstatistik (siehe @sec-anova zu der ANOVA oder @sec-posthoc zu multiplen Vergleichen) -- wir wollen also einen Vergleich zwischen Gruppen rechnen. Wir haben damit dann kein kontinuierliches $x$ vorliegen sondern einen Faktor mit verschiedenen Leveln als Gruppen. - Wir rechnen ein kausales Modell, uns interessieren die Effekte (siehe @sec-simple-kausal) -- wir sind an den Effekten von verschiedenen Einflussvariablen interessiert. Dabei geht es dann weniger um eine Behandlungsgruppe sondern um viele verschiedene Einflussvariablen, die verstanden werden wollen in ihrem Einfluss auf das Outcome $y$. Wir können uns dann die Ergebnisse durch die Funktionen `augment()` und `glance()` aus dem R Paket `{broom}` näher anschauen. - Wir rechnen ein prädiktives Modell, uns interessiert der Wert *neuer* Werte (siehe @sec-simple-pred) -- wir wollen dann also eine Vorhersage rechnen. Wir bauen uns ein komplexes Modell und wollen mit diesem Modell zukünftige Werte vom Outcome $y$ vorhersagen. Wenn wir also neue Einflussvariablen $x$ messen, können wir dann anhand des Modells unbekannte Outcomes $y$ bestimmen. Wir können uns dann mehr Informationen zu der Güte der Vorhersage mit der Funktion `conf_mat()` aus dem R Paket `{tidymodels}` anschauen. ### Inferenzmodell In diesem Kapitel betrachten wir nicht die Analyse von Gruppenvergleichen oder der Inferenzstatistik. Das machen wir dann in eigenen Kapiteln. Bitte besuche das @sec-anova zu mehr Informationen zu der einfaktoriellen und zweifaktoriellen ANOVA oder @sec-posthoc zu multiplen Vergleichen. Es würde dieses Kapitel sprengen, wenn wir uns dann hier die Sachlage nochmal hier aufrollen würde. ### Kausales Modell {#sec-simple-kausal} Im Folgenden rechnen wir ein kausales Modell, da wir an dem Effekt des $x$ interessiert sind. Wenn also das $x_1$ um eine Einheit ansteigt, um wie viel verändert sich dann das $y$? Der Schätzer $\beta_1$ gibt uns also den Einfluss oder den kausalen Zusammenhang zwischen $y$ und $x_1$ wieder. Im ersten Schritt schauen wir uns die Ausgabe der Funktion `lm()` in der Funktion `summary()` an. Daher pipen wir das Objekt `fit_1` in die Funktion `summary()`. ```{r} #| eval: false fit_1 |> summary ``` Wir erhalten folgende Ausgabe dargestellt in @fig-lin-reg-3. ![Die `summary()` Ausgabe des Modells `fit_1`. Hier nicht erschrecken, wir kriegen hier sehr viele Informationen aufeinmal, die wir teilweise nicht alle benötigen. Es hängt dann eben sehr stark von der Fragestelung ab.](images/statistical_modeling_3.png){#fig-lin-reg-3 fig-align="center" width="100%"} Was sehen wir in der Ausgabe der `summary()` Funktion? Als erstes werden uns die Residuen wiedergegeben. Wenn wir nur wenige Beobachtungen haben, dann werden uns die Residuen direkt wiedergegeben, sonst die Verteilung der Residuen. Mit der Funktion `augment()` aus dem R Paket `{broom}` können wir uns die Residuen wiedergeben lassen. Die Residuen schauen wir uns aber nochmal im @sec-lin-reg-quality genauer an. ```{r} fit_1 |> augment() ``` Im zweiten Block erhalten wir die Koeffizienten (eng. *coefficients*) der linearen Regression. Das heißt, wir kriegen dort $\beta_0$ als y-Achsenabschnitt sowie die Steigung $\beta_1$ für das Gewicht. Dabei ist wichtig zu wissen, dass immer als erstes der y-Achsenabschnitt `(Intercept)` auftaucht. Dann die Steigungen der einzelnen $x$ in dem Modell. Wir haben nur *ein* kontinuierliches $x$, daher ist die Interpretation der Ausgabe einfach. Wir können die Gradengleichung wie folgt formulieren. $$ jump\_length \sim 0.97 + 0.51 \cdot weight $$ Was heißt die Gleichung nun? Wenn wir das $x$ um eine Einheit erhöhen dann verändert sich das $y$ um den Wert von $\beta_1$. Wir haben hier eine Steigung von $0.51$ vorliegen. Ohne Einheit keine Interpretation! Wir wissen, dass das Gewicht in \[mg\] gemessen wurde und die Sprungweite in \[cm\]. Damit können wir aussagen, dass wenn ein Floh 1 mg mehr wiegt der Floh 0.51 cm weiter springen würde. Schauen wir nochmal in die *saubere* Ausgabe der `tidy()` Funktion. Wir sehen nämlich noch einen $p$-Wert für den Intercept und die Steigung von `weight`. ```{r} fit_1 |> tidy() ``` Wenn wir einen $p$-Wert sehen, dann brauchen wir eine Nullhypothese, die wir dann eventuell mit der Entscheidung am Signifikanzniveau $\alpha$ von 5% ablehnen können. Die Nullhypothese ist die Gleichheitshypothese. Wenn es also keinen Effekt von dem Gewicht auf die Sprungweite gebe, wie groß wäre dann $\beta_1$? Wir hätten dann keine Steigung und die Grade würde parallel zur x-Achse laufen. Das $\beta_1$ wäre dann gleich null. $$ \begin{align*} H_0: \beta_i &= 0\\ H_A: \beta_i &\neq 0 \\ \end{align*} $$ Wir haben für jedes $\beta_i$ ein eigenes Hypothesenpaar. Meistens interessiert uns der Intercept nicht. Ob der Intercept nun durch die Null geht oder nicht ist eher von geringem Interessen. Spannder ist aber wie sich der $p$-Wert berechnet. Der $p$-Wert basiert auf einer t-Statistik, also auf dem t-Test. Wir rechnen für jeden Koeffizienten $\beta_i$ einen t-Test. Das machen wir in dem wir den Koeffizienten `estimate` durch den Fehler des Koeffizienten `std.error` teilen. $$ \begin{align*} T_{(Intercept)} &= \cfrac{\mbox{estimate}}{\mbox{std.error}} = \cfrac{0.969}{0.445} = 2.18\\ T_{weight} &= \cfrac{\mbox{estimate}}{\mbox{std.error}} = \cfrac{0.510}{0.232} = 2.20\\ \end{align*} $$ Wir sehen in diesem Fall, dass weder der Intercept noch die Steigung von `weight` signifikant ist, da die $p$-Werte mit $0.081$ und $0.079$ leicht über dem Signifikanzniveau von $\alpha$ gleich 5% liegen. Wir haben aber einen starkes Indiz gegen die Nullhypothese, da die Wahrscheinlichkeit die Daten zu beobachten sehr gering ist unter der Annahme das die Nullhypothese gilt. Zun Abschluß noch die Funktion `glance()` ebenfalls aus dem R Paket `{broom}`, die uns erlaubt noch die Qualitätsmaße der linearen Regression zu erhalten. Wir müssen nämlich noch schauen, ob die Regression auch funktioniert hat. Die Überprüfung geht mit einem $x$ sehr einfach. Wir können uns die Grade ja anschauen. Das geht dann mit einem Model mit mehreren $x$ nicht mehr und wir brauchen andere statistische Maßzahlen. ```{r} fit_1 |> glance() ``` ### Prädiktives Modell {#sec-simple-pred} Neben dem kausalen Modell gibt es auch die Möglichkeit ein prädiktives Modell zu rechnen. Im Prinzip ist die Sprache hier etwas ungenau. Wir verwenden das gefittete Modell nur anders. Anstatt das Modell `fit_1` in die Funktion `summary()` zu pipen, pipen wir die das Modell in die Funktion `predict()`. Die Funktion `predict()` kann dann für neue Daten über die Option `newdata =` das $y$ vorhersagen. In unserem Fall müssen wir uns deshalb ein `tibble` mit einer Spalte bauen. Wir haben ja oben im Modell auch nur ein $x_1$ mit aufgenommen. Später können wir natürlich auch für multiple Modelle die Vorhersage machen. Wichtig ist, dass die Namen gleich sind. Das heißt in dem neuen Datensatz müssen die Spalten *exakt* so heißen wir in dem alten Datensatz in dem das Modell gefittet wurde. ```{r} simple_new_tbl <- tibble(weight = c(1.7, 1.4, 2.1, 3.0)) predict(fit_1, newdata = simple_new_tbl) |> round(2) ``` Wie wir sehen ist die Anwendung recht einfach. Wir haben die vier `jump_length` Werte vorhergesagt bekommen, die sich mit dem Fit des Modells mit den neuen `weight` Werten ergeben. In @fig-scatter-lin-pred sehen wir die Visualisierung der vier vorhergesagten Werte. Die Werte müssen auf der Geraden liegen. ```{r} #| echo: false #| message: false #| label: fig-scatter-lin-pred #| fig-align: center #| fig-height: 5 #| fig-width: 6 #| fig-cap: "Scatterplot der *alten* Beobachtungen der Sprungweite in \\[cm\\] und dem Gewicht in \\[mg\\]. Sowie der *neuen* vorhergesagten Beobachtungen auf der Geraden." pred_tbl <- bind_rows(mutate(simple_tbl, status = "beobachtet"), tibble(weight = c(1.7, 1.4, 2.1, 3.0), jump_length = predict(fit_1, newdata = tibble(weight)), status = "vorhergesagt")) ggplot(pred_tbl, aes(weight, jump_length, color = status, shape = status)) + stat_smooth(method = "lm", se = FALSE, fullrange = TRUE, color = "black") + labs(color = "", shape = "") + geom_point(size = 4) + theme_minimal() + scale_color_okabeito() + xlim(0, 3.5) + ylim(0, 3.5) ``` Wir werden später in der [Klassifikation](#sec-class-basic), der Vorhersage von $0/1$-Werten, sowie in der multiplen Regression noch andere Prädktionen und deren Maßzahlen kennen lernen. Im Rahmen der simplen Regression soll dies aber erstmal hier genügen.