071_Regression.Rmd

\part{Geleitetes Modellieren}




# Lineare Regression


```{r echo = FALSE, out.width = "30%", fig.align = "center"}
knitr::include_graphics("images/FOM.jpg")
```

```{r echo = FALSE, out.width = "10%", fig.align = "center"}
knitr::include_graphics("images/licence.png")
```


```{block2, ziele-regression, type='rmdcaution', echo = TRUE} 
Lernziele:


- Wissen, was man unter Regression versteht.
- Die Annahmen der Regression überprüfen können.
- Regression mit kategorialen Prädiktoren durchführen können.
- Die Modellgüte bei der Regression bestimmen können.
- Interaktionen erkennen und ihre Stärke einschätzen können.



```


Für dieses Kapitel benötigen Sie folgende Pakete:
```{r libs-regr}
library(caret)  # Modellieren
library(tidyverse)  # Datenjudo, Visualisierung,...
library(gridExtra)  # Mehrere Plots kombinieren
library(modelr)  # Residuen und Schätzwerte zum Datensatz hinzufügen
library(broom)  # Regressionswerte geordnet ausgeben lassen
```


## Die Idee der klassischen Regression

Regression\index{Regression} ist eine bestimmte Art der *Modellierung* von Daten. Wir legen eine Gerade 'schön mittig' in die Daten; damit haben wir ein einfaches Modell der Daten (vgl. Abb. \@ref(fig:bsp-regression)). Die Gerade 'erklärt' die Daten: Für jeden X-Wert liefert sie einen Y-Wert als Vorhersage zurück.

```{r bsp-regression, fig.cap = "Beispiel für eine Regression"}
stats_test <- read.csv("data/stats_test.csv")

stats_test %>% 
  ggplot +
  aes(x = study_time, y = score) +
  geom_jitter() +
  geom_abline(intercept = 24, 
              slope = 2.3, 
              color = "red")

```

Wie wir genau die Regressionsgerade berechnet haben, dazu gleich mehr. Fürs Erste begnügen wir uns mit der etwas groberen Beobachtung, dass die Gerade 'schön mittig' in der Punktewolke liegt. 

Schauen wir uns zunächst die Syntax genauer an.

```{block2, pseudo-regression, type='rmdpseudocode', echo = TRUE} 

Lade die CSV-Datei mit den Daten als `stats_test`.
  
Nehme `stats_test` UND DANN...  
starte ein neues Diagramm mit ggplot UND      
definiere das Diagramm (X-Achse, Y-Achse) UND DANN  
zeichne das Geom "Jitter" (verwackeltes Punktediagramm) UND DANN  
und zeichne danach eine Gerade ("abline" in rot). 


```

Eine Regression zeigt anhand einer Regressionsgeraden einen "Trend" in den Daten an (s. weitere Beispiele in Abb. \@ref(fig:bsp-regression2)).

```{r bsp-regression2, echo = FALSE, fig.cap = "Zwei weitere Beispiele für Regressionen"}

stats_test %>% 
   ggplot +
  aes(x = self_eval, y = score) +
  geom_jitter() +
  geom_smooth(method = "lm", color = "red", se = FALSE) -> p1


stats_test %>% 
   ggplot +
  aes(x = interest, y = score) +
  geom_jitter() +
  geom_smooth(method = "lm", color = "red", se = FALSE) -> p2


gridExtra::grid.arrange(p1, p2, ncol = 2)


```


Eine Regression lädt förmlich dazu ein, Vorhersagen zu treffen: Hat man erstmal eine Gerade, so kann man für jeden X-Wert ("Prädiktor") eine Vorhersage für den Y-Wert ("Kriterium") treffen. Anhand des Diagramms kann man also für jede Person (d.h. jeden Wert innerhalb des Wertebereichs von `study_time` oder einem anderen Prädiktor) einen Wert für `score` vorhersagen. Wie gut die Vorhersage ist, steht erstmal auf einen anderen Blatt.

Man beachte, dass eine Gerade über ihre *Steigung* und ihren *Achsenabschnitt* festgelegt ist; in Abb. \@ref(fig:bsp-regression) ist die Steigung 2.3 und der Achsenabschnitt 24. Der Achsenabschnitt zeigt also an, wie viele Klausurpunkte man "bekommt", wenn man gar nicht lernt (Gott bewahre); die Steigung gibt eine Art "Wechselkurs" an: Wie viele Klausurpunkte bekomme ich pro Stunde, die ich lerne. 

Unser Modell ist übrigens einfach gehalten: Man könnte argumentieren, dass der Zusatznutzen der 393. Stunde lernen geringer ist als der Zusatznutzen der ersten paar Stunden. Aber dann müssten wir anstelle der Gerade eine andere Funktion nutzen, um die Daten zu modellieren. Lassen wir es erst einmal einfach hier.

Als "Pseudo-R-Formel" ausgedrückt: 
```
score = achsenabschnitt + steigung*study_time
```

Die Vorhersage für die Klausurpunkte (`score`) einer Person sind der Wert des Achsenabschnitts plus das Produkt aus der Anzahl der gelernten Stunden mal den Zusatznutzen pro gelernter Stunde.


Aber wie erkannt man, ob eine Regression "gut" ist - die Vorhersagen also präzise?

In R kann man eine Regression so berechnen:

```{r lm1-stats-test}
lm(score ~ study_time, data = stats_test)
```

`lm` steht dabei für "lineares Modell"; allgemeiner gesprochen lautet die Rechtschreibung für diesen Befehl:

```
lm(kriterium ~ praediktor, data = meine_datentabelle)
```

Um ausführlichere Informationen über das Regressionsmodell zu bekommen, kann man die Funktion `broom::tidy` nutzen:


```{r eval = FALSE}
mein_lm <- lm(kriterium ~ praediktor, data = meine_datentabelle)
tidy(mein_lm)
```


Natürlich kann das auch ~~in der Pfeife rauchen~~ mit der Pfeife darstellen:

```
lm(kriterium ~ praediktor, data = meine_datentabelle) %>% 
  summary
```


## Vorhersagegüte

Der einfache Grundsatz lautet: Je geringer die Vorhersagefehler, desto besser; Abb. \@ref(fig:resids-plot) zeigt ein Regressionsmodell mit wenig Vorhersagefehler (links) und ein Regressionsmodell mit viel Vorhersagefehler (rechts).

```{r resids-plot, echo = FALSE, results = "hold", fig.cap = "Geringer (links) vs. hoher (rechts) Vorhersagefehler"}

set.seed(42)  
N      <- 100
beta   <- 0.4
intercept <- 1


sim <- data_frame(
  x = rnorm(N),
  error1 = rnorm(N, mean = 0, sd = .5),
  error2 = rnorm(N, mean = 0, sd = 2),
  y1 = intercept + x*beta + error1,
  y2 = intercept + x*beta + error2,
  pred = 1 + x*beta
)



p1 <- ggplot(sim, aes(x, y1)) + 
  geom_abline(intercept = intercept, slope = beta, colour = "red") +
  geom_point(colour = "#00998a") +
  geom_linerange(aes(ymin = y1, ymax = pred), colour = "grey40") +
  ylim(-6,+6)


p2 <- ggplot(sim, aes(x, y2)) + 
  geom_abline(intercept = intercept, slope = beta, colour = "red") +
  geom_point(colour = "#00998a") +
  geom_linerange(aes(ymin = y2, ymax = pred), colour = "grey40") +
  ylim(-6,+6)


grid.arrange(p1, p2, ncol = 2)

```


In einem Regressionsmodell lautet die grundlegenden Überlegung zur Modellgüte damit:

>    Wie groß ist der Unterschied zwischen Vorhersage und Wirklichkeit?

Die Größe des Unterschieds (Differenz, "Delta") zwischen vorhergesagten (geschätzten) Wert und Wirklichkeit, bezeichnet man als *Fehler*, *Residuum* oder Vohersagefehler, häufig mit $\epsilon$ (griechisches e wie "error") abgekürzt.


Betrachten Sie die beiden Plots in Abb. \@ref(fig:resids-plot). Die rote Linie gibt die *vorhergesagten* (geschätzten) Werte wieder; die Punkte die *beobachteten* ("echten") Werte. Je länger die blauen Linien, desto größer die Vorhersagefehler. Je größer der Vorhersagefehler, desto schlechter. Und umgekehrt.

>   Je kürzer die typische "Abweichungslinie", desto besser die Vohersage.


Sagt mein Modell voraus, dass Ihre Schuhgröße 49 ist, aber in Wahrheit liegt sie bei 39, so werden Sie dieses Modell als schlecht beurteilen, wahrscheinlich.

Leider ist es nicht immer einfach zu sagen, wie groß der Fehler sein muss, damit das Modell als "gut" bzw. "schlecht" gilt. Man kann argumentieren, dass es keine wissenschaftliche Frage sei, wie viel "viel" oder "genug" ist [@uncertainty]. Das ist zwar plausibel, hilft aber nicht, wenn ich eine Entscheidung treffen muss. Stellen Sie sich vor: Ich zwinge Sie mit der Pistole auf der Brust, meine Schuhgröße zu schätzen.

Eine einfache Lösung ist, das beste Modell unter mehreren Kandidaten zu wählen.

Ein anderer Ansatz ist, die Vorhersage in Bezug zu einem Kriterium zu setzen. Dieses "andere Kriterium" könnte sein "einfach die Schuhgröße raten". Oder, etwas intelligenter, Sie schätzen meine Schuhgröße auf einen Wert, der eine gewisse Plausibilität hat, also z.B. die durchschnittliche Schuhgröße des deutschen Mannes. Auf dieser Basis kann man dann quantifizieren, ob und wie viel besser man als dieses Referenzkriterium ist.

### Mittlere Quadratfehler
Eine der häufigsten Gütekennzahlen ist der *mittlere quadrierte Fehler* (engl. "mean squared error", MSE), wobei Fehler wieder als Differenz zwischen Vorhersage (`pred`) und beobachtete Wirklichkeit (`obs`, `y`) definiert ist. Dieser berechnet für jede Beobachtung den Fehler, quadriert diesen Fehler und bilden dann den Mittelwert dieser "Quadratfehler", also einen *mittleren Quadratfehler*. Die englische Abkürzung *MSE* ist auch im Deutschen gebräuchlich.

$$ MSE = \frac{1}{n} \sum{(pred - obs)^2} $$

Konzeptionell ist dieses Maß an die Varianz angelehnt. Zieht man aus diesem Maß die Wurzel, so erhält man den sog. *root mean square error* (RMSE), welchen man sich als die Standardabweichung der Vorhersagefehler vorstellen kann. In Pseudo-R-Syntax:

```
RMSE <- sqrt(mean((df$pred - df$obs)^2))
```

Der RMSE hat die selben Einheiten wie die zu schätzende Variable, also z.B. Schuhgrößen-Nummern.




### R-Quadrat ($R^2$)
$R^2$, auch *Bestimmtheitsmaß*\index{Bestimmtheitsmaß} oder *Determinationskoeffizient*\index{Determinationskoeffizient} genannt, setzt die Höhe unseres Vorhersagefehlers\index{Vorhersagefehler} im Verhältnis zum Vorhersagefehler eines "Nullmodell". Das Nullmodell hier würde sagen, wenn es sprechen könnte: "Keine Ahnung, was ich schätzen soll, mich interessieren auch keine Prädiktoren, ich schätzen einfach immer den Mittelwert der Grundgesamtheit!".

Analog zum Nullmodell-Fehler spricht auch von der Gesamtvarianz oder $SS_T$ (sum of squares total); beim Vorhersagefehler des eigentlichen Modells spricht man auch von $SS_M$ (sum of squares model).

Damit gibt $R^2$ an, wie gut unsere Vorhersagen im Verhältnis zu den Vorhersagen des Nullmodells sind. Ein $R^2$ von 25% (0.25) hieße, dass unser Vorhersagefehler 25% *kleiner* ist als der der Nullmodells. Ein $R^2$ von 100% (1) heißt also, dass wir den kompletten Fehler reduziert haben (Null Fehler übrig) - eine perfekte Vorhersage. Etwas formaler, kann man $R^2$ so definieren:

$$ R^2 = 1 - \left( \frac{SS_T - SS_M}{SS_T} \right)$$

Präziser, in R-Syntax:

```
R2 <- 1 - sum((df$pred - df$obs)^2) / sum((mean(df$obs) - df$obs)^2)
```


Praktischerweise gibt es einige R-Pakete, z.B. `caret`, die diese Berechnung für uns besorgen:

```{r R2-caret, eval = FALSE}
postResample(obs = obs, pred = pred)
```

Hier steht `obs` für beobachtete Werte und `pred` für die vorhergesagten Werte (beides numerische Vektoren). Dieser Befehl gibt sowohl RMSE als auch $R^2$ wieder. Wir betrachten gleich ein Beispiel an echten Daten.




```{block2, r-nicht-als-guete, type='rmdcaution', echo = TRUE}

Verwendet man die Korrelation (r) oder $R^2$ als Gütekriterium, so sollte man sich über folgenden Punkt klar sein. Bei Skalierung der Variablen ändert sich die Korrelation nicht; das gilt auch für $R^2$. Beide Koeffizienten ziehen allein auf das *Muster* der Zusammenhänge ab - nicht die Größe der Abstände. Aber häufig ist die Größe der Abstände zwischen beobachteten und vorhergesagten Werten das, was uns interessiert. In dem Fall wäre der MSE vorzuziehen.

```



## Die Regression an einem Beispiel erläutert

Schauen wir uns den Datensatz zur Statistikklausur noch einmal an. Welchen Einfluss hat die Lernzeit auf den Klausurerfolg? Wie viel bringt es also zu lernen? Wenn das Lernen keinen Einfluss auf den Klausurerfolg hat, dann kann man es ja gleich sein lassen... Aber umgekehrt, wenn es viel bringt, ok gut, dann könnte man sich die Sache (vielleicht) noch mal überlegen. Aber was heißt "viel bringen" eigentlich?

>   Wenn für jede Stunde Lernen viele zusätzliche Punkte herausspringen, dann bringt Lernen viel. Allgemeiner: Je größer der Zuwachs im Kriterium ist pro zusätzliche Einheit des Prädiktors, desto größer ist der Einfluss des Prädiktors.

Natürlich könnte jetzt jemand argumentieren, dass die ersten paar Stunden lernen viel bringen, aber dann flacht der Nutzen ab, weil es ja schnell einfach und trivial wird. Aber wir argumentieren (erstmal) so nicht. Wir gehen davon aus, dass jede Stunde Lernen gleich viel (oder wenig) Nutzen bringt.

>   Geht man davon aus, dass jede Einheit des Prädiktors gleich viel Zuwachs bringt, unabhängig von dem Wert des Prädiktors, so geht man von einem linearen Einfluss aus.


Versuchen wir im ersten Schritt die Stärke des Einfluss an einem Streudiagramm abzuschätzen (s. Abb. \@ref(fig:bsp-regression)).


Hey R - berechne uns die "Trendlinie"! Dazu nimmt man den Befehl `lm`:


```{r}
mein_lm <- lm(score ~ study_time, data = stats_test)
tidy(mein_lm)

```

`lm` steht für 'lineares Modell', eben weil eine *Linie* als Modell in die Daten gelegt wird. Aha. Die Steigung der Geraden beträgt 2.3 - das ist der Einfluss des Prädiktors Lernzeit auf das Kriterium Klausurerfolg! Man könnte sagen: Der "Wechselkurs" von Lernzeit auf Klausurpunkte. Für jede Stunde Lernzeit bekommt man offenbar 2.3 Klausurpunkte (natürlich viel zu leicht). Wenn man nichts lernt (`study_time == 0`) hat man 24 Punkte.

>    Der Einfluss des Prädiktors steht unter 'estimate'. Der Kriteriumswert wenn der Prädiktor Null ist steht unter '(Intercept)'.



Malen wir diese Gerade in unser Streudiagramm (Abbildung \@ref(fig:stats-test-scatter2).


```{r stats-test-scatter2, fig.cap = "Streudiagramm von Lernzeit und Klausurerfolg"}
ggplot(data = stats_test) +
  aes(y = score, x = study_time) +
  geom_jitter() +
  geom_abline(slope = 2.3, intercept = 24, color = "red")
```

Jetzt kennen wir die Stärke (und Richtung) des Einflusses der Lernzeit. Ob das viel oder wenig ist, ist am besten im Verhältnis zu einem Referenzwert zu sagen.

Die Gerade wird übrigens so in die Punktewolke gelegt, dass die (quadrierten) Abstände der Punkte zur Geraden minimal sind. Dies wird auch als *Kriterium der Kleinsten Quadrate*\index{Kriterium der Kleinsten Quadrate} (*Ordinary Least Squares*, *OLS*)\index{Ordinary Least Squares} bezeichnet.

Jetzt können wir auch einfach Vorhersagen machen. Sagt uns jemand, ich habe "viel" gelernt (Lernzeit = 4), so können wir den Klausurerfolg grob im Diagramm ablesen.

Genauer geht es natürlich mit dieser Rechnung:

$y = 4*2.3 + 24$

Oder mit diesem R-Befehl:



```{r}
predict(mein_lm, data.frame(study_time = 4))
```



Berechnen wir noch die Vorversagegüte des Modells. Dazu kann man den Befehl `summary` nehmen, oder auch `broom::glance`. `glance` gibt Informationen zur Modellgüte zurück und das in Form eines Dateframes. Summary liefert eine Menge Informationen mit einem infomrationen Ausdruck, aber nicht in Form eines Dataframes (sondern in Form einer Liste).

```{r}
glance(mein_lm)
```


Das Bestimmtheitsmaß $R^2$ ist mit `r round(summary(mein_lm)$r.squared,2)` "ok": `r round(summary(mein_lm)$r.squared*100)`-\% der Varianz des Klausurerfolg wird im Modell 'erklärt'. 'Erklärt' meint hier, dass wenn die Lernzeit konstant wäre, würde die Varianz von Klausurerfolg um diesen Prozentwert sinken.


## Überprüfung der Annahmen der linearen Regression

Aber wie sieht es mit den Annahmen aus?

- Die *Linearität des Zusammenhangs* haben wir zu Beginn mit Hilfe des Scatterplots überprüft. Es schien einigermaßen zu passen.

- Zur Überprüfung der *Normalverteilung der Residuen* zeichnen wir ein Histogramm (s. Abbildung \@ref(fig:resid-distrib)). Die *Residuen*\index{Residuen} können über den Befehl `add_residuals` (Paket `modelr`) zum Datensatz hinzugefügt werden. Dann wird eine Spalte mit dem Namen `resid` zum Datensatz hinzugefügt.  

Hier scheint es zu passen:

```{r resid-distrib, fig.cap = "Die Residuen verteilen sich hinreichend normal."}
stats_test %>% 
  add_residuals(mein_lm) %>% 
  ggplot +
  aes(x = resid) +
  geom_histogram()
```


Sieht passabel aus. Übrigens kann man das Paket `modelr` auch nutzen, um sich komfortabel die vorhergesagten Werte zum Datensatz hinzufügen zu lassen (Spalte `pred`):

```{r}
stats_test %>% 
  add_predictions(mein_lm) %>% 
  select(pred) %>% 
  head
```


- *Konstante Varianz*: Dies kann z. B. mit einem Scatterplot der Residuen auf der Y-Achse und den vorhergesagten Werten auf der X-Achse überprüft werden. Bei jedem X-Wert sollte die Varianz der Y-Werte (etwa) gleich sein (s. Abbildung \@ref(fig:tips-preds-resid)).

Die geschätzten (angepassten) Werte kann man über den Befehl `add_predictions()` aus dem Paket `modelr` bekommen. Die Fehlerwerte entsprechend mit dem Befehl `add_residuals`. 

```{r tips-preds-resid, fig.cap = "Vorhergesagte Werte vs. Residualwerte im Datensatz tips"}

stats_test %>% 
  add_predictions(mein_lm) %>% 
  add_residuals(mein_lm) %>% 
  ggplot() +
  aes(y = resid, x = pred) +
  geom_point()

```


Die Annahme der konstanten Varianz scheint verletzt zu sein: Die sehr großen vorhersagten Werte können recht genau geschätzt werden; aber die mittleren Werte nur ungenau.
Die Verletzung dieser Annahme beeinflusst *nicht* die Schätzung der Steigung, sondern die Schätzung des Standardfehlers, also des p-Wertes der Einflusswerte.

- *Extreme Ausreißer*: Extreme Ausreißer scheint es nicht zu geben.

- *Unabhängigkeit der Beobachtungen*: Wenn die Studenten in Lerngruppen lernen, kann es sein, dass die Beobachtungen nicht unabhängig voneinander sind: Wenn ein Mitglied der Lerngruppe gute Noten hat, ist die Wahrscheinlichkeit für ebenfalls gute Noten bei den anderen Mitgliedern der Lerngruppe erhöht. Böse Zungen behaupten, dass 'Abschreiben' eine Gefahr für die Unabhängigkeit der Beobachtungen sei.



```{block2, tips-uebung, type='rmdexercises', echo = TRUE}
1.  Wie groß ist der Einfluss des Interesss?

2.  Für wie aussagekräftig halten Sie Ihr Ergebnis aus 1.?

3.  Welcher Einflussfaktor (in unseren Daten) ist am stärksten?

```



## Regression mit kategorialen Prädiktoren
Vergleichen wir interessierte und nicht interessierte Studenten. Dazu teilen wir die Variable `interest` in zwei Gruppen (1-3 vs. 4-6) auf:

```{r}
stats_test$interessiert <- stats_test$interest > 3
```


Vergleichen wir die Mittelwerte des Klausurerfolgs zwischen den Interessierten und Nicht-Interessierten:

```{r}
stats_test %>% 
  group_by(interessiert) %>% 
  summarise(score = mean(score)) -> score_interesse

score_interesse
```

Aha, die Interessierten haben im Schnitt mehr Punkte; aber nicht viel.

```{r}
stats_test %>% 
  na.omit %>% 
  ggplot() +
  aes(x = interessiert, y = score) +
  geom_jitter(width = .1) +
  geom_point(data = score_interesse, color = "red", size = 5) +
  geom_line(data = score_interesse, group = 1, color = "red")
```

Mit `group=1` bekommt man eine Linie, die alle Punkte verbindet (im Datensatz `score_interesse` sind es dieser zwei). Wir haben in dem Fall nur zwei Punkte, die entsprechend verbunden werden.


###  Aufgaben
1. Visualisieren Sie den Gruppenunterschied auch mit einem Boxplot!


2. Berechnen Sie ein lineares Modell dazu!


*Lösung:*

1. Boxplot: `qplot(x = interessiert, y = score, data = stats_test, geom = "boxplot")`


2. Lineares Modell:


```{r}
lm2 <- lm(score ~ interessiert, data = stats_test)
summary(lm2)
```

Der Einfluss von `interessiert` ist statistisch signifikant (p = .03). Der Stärke des Einflusses ist im Schnitt 1.6 Klausurpunkte (zugunsten `interessiertTRUE`). Das ist genau, was wir oben herausgefunden haben.


```{block2, tips-uebung2, type='rmdexercises', echo = TRUE}
3.  Vie ist der Einfluss von `study_time`, auch in zwei Gruppen geteilt?

4.  Wie viel \% der Variation des Klausurerfolgs können Sie durch das Interesse modellieren?

```




## Multiple Regression
Aber wie wirken sich mehrere Einflussgrößen *zusammen* auf den Klausurerfolg aus?


```{r}
lm3 <- lm(score ~ study_time + interessiert, data = stats_test)
summary(lm3)
```

Interessant ist das *negative* Vorzeichen vor dem Einfluss von `interessiertTRUE`! Die multiple Regression untersucht den 'Nettoeinfluss' jedes Prädiktors. Den Einfluss also, wenn der andere Prädiktor *konstant* gehalten wird. Anders gesagt: Betrachten wir jeden Wert von `study_time` separat, so haben die Interessierten jeweils im Schnitt etwas *weniger* Punkte (jesses). Allerdings ist dieser Unterschied nicht statistisch signifikant.

>   Die multiple Regression zeigt den 'Nettoeinfluss' jedes Prädiktor: Den Einfluss dieses Prädiktor, wenn der andere Prädiktor oder die anderne Prädiktoren konstant gehalten werden.


Hier haben wir übrigens dem Modell aufgezwungen, dass der Einfluss von Lernzeit auf Klausurerfolg bei den beiden Gruppen gleich groß sein soll  (d.h. bei Interessierten und Nicht-Interessierten ist die Steigung der Regressionsgeraden gleich). Das illustriert sich am einfachsten in einem Diagramm (s. Abbildung \@ref(fig:no-interakt)).


```{r no-interakt, echo = FALSE, fig.cap = "Eine multivariate Analyse fördert Einsichten zu Tage, die bei einfacheren Analysen verborgen bleiben", fig.height = 7}

library(viridis)

df1 <- data_frame(
  interessiert = c(T, F),
  slope = c(2.3, 2.3),
  intercept = c(23.7, 24)
)


p1 <- stats_test %>% 
  na.omit %>% 
  ggplot +
  aes(x=study_time, y = score) +
  geom_jitter(width = .1) +
  #facet_wrap(~interessiert) +
  geom_abline(data = df1, mapping = aes(slope = slope, intercept = intercept, color = interessiert)) +
  guides(color = FALSE) +
  scale_color_viridis(discrete = TRUE) +
  labs(title = "A") +
  coord_cartesian(ylim = c(25,35))


df2 <- data_frame(
  slope = rep(-.33, 5),
  intercept = rep(c(24), 5),
  study_time = 1:5)


p2 <- stats_test %>% 
  na.omit %>% 
  ggplot +
  aes(x=interessiert, y = score) +
  geom_jitter(width = .1) +
  facet_wrap(~study_time) +
  geom_abline(data = df2, mapping = aes(slope = slope, intercept = intercept, color = study_time)) +
  scale_color_viridis() +
  guides(color = FALSE) +
  labs(title = "B")


gridExtra::grid.arrange(p1, p2, nrow = 2)
```


Diese *multivariate*\index{multivariat} Analyse (mehr als 2 Variablen sind beteiligt) zeigt uns, dass die Regressionsgerade nicht gleich ist in den beiden Gruppen (Interessierte vs. Nicht-Interessierte; s. Abbildung \@ref(fig:no-interakt)): Im Teildiagramm A sind die Geraden (leicht) versetzt. Analog zeigt Teildiagramm B, dass die Interessierten (`interessiert == TRUE`) geringe Punktewerte haben als die Nicht-Interessierten, wenn man die Werte von `study_time` getrennt betrachtet.

>    Die multivariate Analyse zeigt ein anderes Bild, ein genaueres Bild als die einfachere Analyse. Ein Sachverhalt, der für den ganzen Datensatz gilt, kann in Subgruppen anders sein.

Ohne multivariate Analyse hätten wir dies nicht entdeckt. Daher sind multivariate Analysen sinnvoll und sollten gegenüber einfacheren Analysen bevorzugt werden.

Man könnte sich jetzt noch fragen, ob die Regressionssgerade in Abbildung \@ref(fig:no-interakt) parallel sein müssen. Gerade hat unser R-Befehl sie noch gezwungen, parallel zu sein. Gleich lassen wir hier die Zügel locker. Wenn die Regressionsgerade nicht mehr parallel sind, spricht man von *Interaktionseffekten*.


Das Ergebnis des zugrunde-liegenden F-Tests (vgl. Varianzanalyse) wird in der letzten Zeile angegeben (`F-Statistic`). Hier wird $H_0$ also verworfen.




## Interaktionen 

Es könnte ja sein, dass die Stärke des Einflusses von Lernzeit auf Klausurerfolg in der Gruppe der Interessierten anders ist als in der Gruppe der Nicht-Interessierten. Wenn man nicht interessiert ist, so könnte man argumentieren, dann bringt eine Stunden Lernen weniger als wenn man interessiert ist. Darum müssten die Steigungen der Regressionsgeraden in den beiden Gruppen unterschiedlich sein. Schauen wir uns es an. Um R dazu zu bringen, die Regressionsgeraden frei variieren zu lassen, so dass sie nicht mehr parallel sind, nutzen wir das Symbol `*`, dass wir zwischen die betreffenden Prädiktoren schreiben:

```{r}
lm4 <- lm(score ~ interessiert*study_time, data = stats_test)
summary(lm4)
```

Interessanterweise zeigen die Interessierten nun wiederum - betrachtet man jede Stufe von `study_time` einzeln - bessere Klausurergebnisse als die Nicht-Interessierten. Ansonsten ist noch die Zeile `interessiertTRUE:study_time` neu. Diese Zeile zeigt die Höhe des *Interaktionseffekts*\index{Interaktionseffekt}. Bei den Interessierten ist die Steigung der Geraden um 0.32 Punkte geringer als bei den Interessierten. Der Effekt ist klein und nicht statistisch signifikant, so dass wir wahrscheinlich Zufallsrauschen überinterpretieren. Aber die reine Zahl sagt, dass bei den Interessierten jede Lernstunde weniger Klausurerfolg bringt als bei den Nicht-Interessierten. Auch hier ist eine Visualisierung wieder hilfreich. 



```{r interakt-stats-test, echo = FALSE, fig.cap = "Eine Regressionsanalyse mit Interaktionseffekten"}



df1 <- data_frame(
  interessiert = c(F, T),
  slope = c(2.44, 2.44-.32),
  intercept = c(23.6, 23.6+.66)
)


p1 <- stats_test %>% 
  na.omit %>% 
  ggplot +
  aes(x=study_time, y = score) +
  geom_jitter(width = .1) +
  #facet_wrap(~interessiert) +
  geom_abline(data = df1, mapping = aes(slope = slope, intercept = intercept, color = interessiert)) +
  guides(color = FALSE) +
  scale_color_viridis(discrete = TRUE) +
  labs(title = "A") +
  coord_cartesian(ylim = c(25,35))



df2 <- stats_test %>% 
  data_grid(
    study_time = 1:5,
    interessiert = c(F, T)
  ) %>% 
  add_predictions(lm4)
  


p2 <- stats_test %>% 
  na.omit %>% 
  ggplot +
  aes(x=interessiert, group = study_time) +
  geom_jitter(aes(y = score), width = .1) +
  geom_line(data = df2, aes(y = pred, color = factor(study_time))) +
  scale_color_viridis(discrete = TRUE) +
  theme(legend.position = "bottom") +
  ggtitle("B")
  


gridExtra::grid.arrange(p1, p2, nrow = 2)
```


Wir wir in Abbildung \@ref(fig:interakt-stats-test) sehen, ist der Einfluss von `study_time' je nach Gruppe (Wert von `interessiert`) unterschiedlich (Teildiagramm A). Analog ist der Einfluss des Interesses (leicht) unterschiedlich, wenn man die fünf Stufen von `study_time` getrennt betrachtet. 

>   Sind die Regressionsgerade nicht parallel, so liegt ein Interaktionseffekt vor. Andernfalls nicht.






## Fallstudie zu Overfitting {#overfitting-casestudy}

Vergleichen wir im ersten Schritt eine Regression, die die Modellgüte anhand der *Trainingsstichprobe* schätzt mit einer Regression, bei der die Modellgüte in einer *Test-Stichprobe* überprüft wird.


Betrachten wir nochmal die einfache Regression von oben. Wie lautet das $R^2$?

```{r lm-overfitting1}

lm1 <- lm(score ~ study_time, data = stats_test)

```

Es lautet `round(summary(lm1)$r.squared, 2)`.

Im zweiten Schritt teilen wir die Stichprobe in eine Trainings- und eine Test-Stichprobe auf. Wir "trainieren" das Modell anhand der Daten aus der Trainings-Stichprobe:

```{r lm-overfitting2}
train <- stats_test %>% 
  sample_frac(.8, replace = FALSE)  # Stichprobe von 80%, ohne Zurücklegen

test <- stats_test %>% 
  anti_join(train)  # Alle Zeilen von "df", die nicht in "train" vorkommen

lm_train <- lm(score ~ study_time, data = train)
```


Dann testen wir (die Modellgüte) anhand der *Test*-Stichprobe. Also los, `lm_train`, mach Deine Vorhersage:

```{r lm-overfitting-predict}
lm2_predict <- predict(lm_train, newdata = test)
```

Diese Syntax sagt:

```{block2, lm2-predict-block, type='rmdpseudocode', echo = TRUE}
Speichere unter dem Namen "lm2_predict" das Ergebnis folgender Berechnung:  
Mache eine Vorhersage ("to predict") anhand des Modells "lm2",   
wobei frische Daten ("newdata = test") verwendet werden sollen. 

```

Als Ergebnis bekommen wir einen Vektor, der für jede Beobachtung des Test-Samples den geschätzten (vorhergesagten) Klausurpunktewert speichert.

```{r R2-postresample}
caret::postResample(pred = lm2_predict, obs = test$score)

```

Die Funktion `postResample` aus dem Paket `caret` liefert uns zentrale Gütekennzahlen unser Modell. Wir sehen, dass die Modellgüte im Test-Sample deutlich *schlechter* ist als im Trainings-Sample. Ein typischer Fall, der uns warnt, nicht vorschnell optimistisch zu sein!


>   Die Modellgüte im in der Test-Stichprobe ist meist schlechter als in der Trainings-Stichprobe. Das warnt uns vor Befunden, die naiv nur die Werte aus der Trainings-Stichprobe berichten.


## Aufgaben^[F, R, R, F, F, F, F, F, R]



```{block2, exercises-regr, type='rmdexercises', echo = TRUE} 
Richtig oder Falsch!?

1. X-Wert: Kriterium; Y-Wert: Prädiktor.

1. Der Y-Wert in der einfachen Regression wird berechnet als Achsenabschnitt plus *x* mal die Geradensteigung.

1. $R^2$ liefert einen *relativen* Vorhersagefehler und MSE einen *absoluten* (relativ im Sinne eines Anteils).

1. Unter 'Ordinary Least Squares' versteht man eine abschätzige Haltung gegenüber Statistik.

5. Zu den Annahmen der Regression gehört Normalverteilung der *Kriteriumswerte*.

1. Die Regression darf nicht bei kategorialen Prädiktoren verwendet werden.

1. Mehrere bivariate Regressionsanalysen (1 Prädiktor, 1 Kriterium) sind einer multivariaten Regression i.d.R. vorzuziehen.

1. Interaktionen erkennt man daran, dass die Regressionsgeraden *nicht* parallel sind.

```


## Befehlsübersicht

Tabelle \@ref(tab:befehle-regression) stellt die Befehle dieses Kapitels dar. 


```{r befehle-regression, echo = FALSE}

df <- readr::read_csv("includes/Befehle_Regression.csv")


knitr::kable(data.frame(df), caption = "Befehle des Kapitels 'Regression'")




```