outlier_vision.Rmd

---
title: "Ausreisser"
output: 
  html_document:
    includes:
      after_body: footer.html 
---

# Datensatz
In diesem Beispiel wurden Name, Alter, Geschlecht und Sehvermögen (auf einer Skala von 1 bis 10) von ingsesamt 29 Personen aufgenommen. Wir wollen so tun als hätten wir selbst den Datensatz erhoben.

```{r, include=FALSE}
library(data.table)
library(dplyr)
dat <- fread("D:/RKurse/MyDatasets/vision fixed (other) LM.txt", 
             dec=",", na.strings = "", stringsAsFactors=TRUE)
dat <- dat %>% select(Person, Ages, Gender, Vision)
```

<div class = "row"> <div class = "col-md-6">
```{r}
head(dat) # Erste Zeilen des Datensatzes
```
</div> <div class = "col-md-6">
```{r, eval=FALSE}
str(dat) # Struktur des Datensatzes
```
```{r, echo=FALSE, warning=FALSE, message=FALSE, error=FALSE}
str(dat, width=40, strict.width="cut")
```
</div> </div>

# Fragestellung
Genau wie im [Alkohol-Beispiel](1n_drinks.html) könnten wir nun unterschiedliche Fragen an den Datensatz stellen, wie z.B. ob in dieser Stichprobe die Frauen oder Männer ein höheres Durchschnittsalter haben. Doch wieder wollen wir uns hier nur auf die *zwei numerischen Variablen* `Ages` und `Vision` beschränken und demnach fragen ob das Alter einer Person etwas mit dessen Sehvermögen zu tun hat. Ebenfalls analog zum [Alkohol-Beispiel](1n_drinks.html) soll die Korrelation und eine Regressionsgerade geschätzt werden. Der Unterschied zum vorangegangenen Beispiel wird sich während der Analyse weiter unten zeigen.

# Deskriptive Statistik
Um ein Gefühl für die Daten zu bekommen, betrachten wir einige Kennzahlen zu den Daten, sowie einen Plot.

<div class = "row"> <div class = "col-md-6">
```{r, fig.height = 3, fig.width = 4, fig.align = "center"}
summary(dat[, c("Ages", "Vision")])
```
</div> <div class = "col-md-6">
```{r}
plot(x=dat$Ages, y=dat$Vision)
```
</div> </div>

### Fehlwerte?
Beim Betrachten des `summary()` Outputs fällt auf, dass es Fehlwerte (=`NA`) im Datensatz gibt. Wir prüfen diese indem wir uns die Datenzeile näher anschauen.

```{r}
subset(dat, is.na(Ages)) # Zeige nur die Daten mit Fehlwert in der Ages-Splate
```

Für Enrique haben wir demnach keine weitere Angabe außer seinen Namen. In solch einem Fall würden wir nun versuchen uns zu erinnern warum das so ist und ob wir die Daten nicht doch noch nachträglich beschaffen könnten. Wenn uns aber beispielsweise einfällt, dass Enrique während der Datenaufnahme unerwartet abbrechen musste und wir ihn so schnell auch nicht wieder erreichen können, dann ist dieser Datenpunkt zumindest für diese Analyse komplett wertlos. Tatsächlich kann man sogar argumentieren, dass er verwirrt, da der Datensatz 29 Zeilen hat, in den Analysen aber nur 28 davon brauchbar sein werden. Wir entscheiden uns also in diesem Fall die Datenzeile deshalb sicherheitshalber zu löschen.

An dieser Stelle zeige ich gleich drei verschiedene Möglichkeiten alle Datenzeilen zu löschen, die Fehlwerte in der `Ages` Spalte haben (= die eine Zeile mit Enrqiue). Für mehr Infos zur zweiten und dritten Methode siehe [das Kapitel mit weiteren Tipps zu R](datr_moreadvanced.html):

```{r, include=FALSE}
dat <- dat[!is.na(Ages)]
```

<div class = "row"> <div class = "col-md-6">
```{r, eval=FALSE}
dat <- subset(dat, !is.na(Ages)) # Basis R
dat <- dat[!is.na(Ages)] # data.table
dat <- dat %>% subset(!is.na(Ages)) # dplyr
```
</div> <div class = "col-md-6">
```{r}
nrow(dat) # Anzahl Zeilen nach der Löschung
```
</div> </div>

### Negativer Trend
Nun können wir uns die Deskriptive genauer anschauen und stellen beim Betrachten des Plots fest, dass es scheinbar einen negativen Trend gibt: Mit zunehmendem Alter sinkt das Sehvermögen. Was man sich an diesem Punkt auch vergegenwärtigen sollte ist, dass die jüngste Person 22 und die älteste 55 Jahre alt ist. Das bedeutet für uns, dass all unsere Ergebnisse am Ende streng genommen nur für Personen in dieser Altersspanne gelten. Nähmen wir z.B. das später geschätzte Regressionsmodell um das Sehvermögen einer 4-jährigen Person zu vorherzusagen, dann würden wir *"extrapolieren"*, da dieses Alter klar außerhalb der Daten liegt, welche zur Bildung des Modells vorhanden waren. Solche Vorhersagen sind natürlich stets mit Vorsicht zu genießen.

# Korrelation & Regression
Wie im vorangegangenen Kapitel wollen wir also Korrelation und Regression berechnen. Beim Betrachten des Plots oben können wir schon jetzt sagen, dass wir eine negative, wenn auch nicht ganz so starke Korrelation und eine negative Steigung erwarten:

```{r}
cor <- cor.test(dat$Ages, dat$Vision) # Korrelation
cor
reg <- lm(data=dat, formula=Vision ~ Ages) # Regressionsmodell y = a + bx
summary(reg)
```

Und wir stellen beides wieder in einem Plot mit `ggplot` dar:

```{r, message=FALSE, warning=FALSE}
library(ggplot2) # Erzeugt Plot
library(ggpubr) # Zusätzlicher Befehl "stat_cor" und "stat_regline_equation" (siehe unten)
```
```{r}
ggplot(data=dat, aes(x=Ages, y=Vision)) + # Definiere Daten
  geom_point(size=3) + # Scatter plot mit Punkten der Größe 3
  scale_x_continuous(name="Alter",       limits=c(20, 60), breaks=seq(20, 60, by=10)) + # x-Achse
  scale_y_continuous(name="Sehvermögen", limits=c( 0, 10), breaks=seq( 0, 10, by=1)) + # y-Achse
  theme_classic() + # Simple, klassische Formatierung
  geom_smooth(method='lm', formula=y~x, se=FALSE) + # Füge Gerade ein
  stat_regline_equation(label.x=50, label.y=9) + # Füge Gleichung ein
  stat_cor(method="pearson", label.x=40, label.y=3) # Füge Korrelation mit p-Wert ein
```

# Der Ausreißer
### Schritt 1: Finde eventuelle Ausreißer
An diesem Punkt könnte man das Beispiel als abgehakt ansehen. Es fällt aber auf, dass ein Punkt deutlich aus der Reihe tanzt - nämlich die Person mit Sehvermögen von 3. Diese Person hat das mit Abstand schlechteste Sehvermögen (das zweitschlechteste liegt bei 6) und hinzukommt, dass die Person auch noch zu den jungen Personen gehört, sodass dieser Wert komplett gegen den Trend der restlichen Daten geht. 

### Schritt 2: Mehr über Ausreißer erfahren
In solch einem Fall sollte immer erst mehr über den Ausreißer in Erfahrung gebracht werden:
```{r}
subset(dat, Vision<5) # Zeige nur Datenzeilen mit Vision < 5
```
Es handelt sich also um den 26-jährigen Mann namens Rolando. Hier eine Möglichkeit um ihn in `ggplot` mittels des packages `ggrepel` hervorzuheben. Im Code sieht man, dass nur die Punkte gekennzeichnet werden sollen, die eine Vision<5 haben, sodass nur Rolando gekennzeichnet wird.

```{r, message=FALSE, warning=FALSE}
library(ggrepel) # Zusätzlicher Befehl "geom_text_repel" (siehe unten)
```
```{r}
ggplot(data=dat, aes(x=Ages, y=Vision)) + # Definiere Daten
  geom_point(size=3) + # Scatter plot mit Punkten der Größe 3
  scale_x_continuous(name="Alter",       limits=c(20, 60), breaks=seq(20, 60, by=10)) + # x-Achse
  scale_y_continuous(name="Sehvermögen", limits=c( 0, 10), breaks=seq( 0, 10, by=1)) + # y-Achse
  theme_classic() + # Simple, klassische Formatierung
  geom_smooth(method='lm', formula=y~x, se=FALSE) + # Füge Gerade ein
  geom_text_repel(data=dat[Vision<5], aes(label=Person), nudge_y=1, nudge_x=1) # Kennzeichne Rolando
```

### Schritt 3: Entscheide über Ausreißer
Hätten wir die Daten selbst erhoben, würde wir nun versuchen uns an Rolando zu erinnern. 

**Option 1: Klingt komisch, ist aber so:** Würden wir direkt jemanden mit dicken Brillengläsern vor Augen haben und uns erinnern, dass er überraschend schlecht gesehen hat, hätten wir den auffälligen Wert bestätigt. Wir würden nichts am Wert ändern und er würde weiterhin dem generellen Trend widersprechen, aber das ändert nichts daran, dass er die Wahrheit abbildet.

**Option 2: Korrigiere den Fehler:** Könnten wir uns allerdings an Rolando erinnern und wären uns ziemlich sicher, dass er eigentlich ein normales Sehrvermögen hatte, dann liegt ein Problem vor. Nicht selten kommt es vor, dass man sich z.B. bei der Dateneingabe vertippt und deshalb ein falscher Wert im Datensatz steht. Könnten wir nun zurück zu unseren schriftlichen Aufzeichnungen vom Tag der Datenerhebung gehen und dort Rolandos korrekten Wert für das Sehvermögen einsehen, so könnten wir nachträglich den Wert im Datensatz korrigieren und fortfahren.

**Option 3: Lösche den Fehler:** Es kann aber auch sein, dass wir uns wie bei Option 2 sicher sind, dass der Wert nicht bei 3 lag, wir aber keine Möglichkeit haben den korrekten Wert aufzutreiben. In solch einem Fall kann man sich dafür entscheiden den Wert zu löschen. Dann haben wir zwar weniger Daten, aber zumindest ist das Ergebnis nicht durch den offensichtlich falschen Wert verfälscht. 

> Das Löschen von Ausreißern aus einem Datensatz ist eine sensible Angelegenheit und sollte nicht leichtfertig und schon gar nicht heimlich durchgeführt werden. Leider kommt es vor, dass selbst Wissenschaftler mehr oder weniger subjektiv entscheiden welche Werte "schlichtweg nicht wahr sein können" und diese dann aus dem Datensatz entfernen ohne es in der Publikation zu vermerken. In manchen Fällen führt dies sogar dazu, dass das Ergebnis erst danach signifikant wird. Mehr dazu im [Kapitel zum p(roblem)-Wert](stat_pwert.html).

# Korrelation & Regression - nochmal
In diesem Beispiel wollen wir uns dazu entscheiden den Ausreißer zu löschen. Demnach entfernen wir Rolando aus dem Datensatz und führen beide Analysen nochmals durch. Wir speichern die Ergebnisse diesmal aber in neue R-Objekte, sodass wir die alten nicht überschreiben:

```{r}
dat_noRolando <- subset(dat, Person!="Rolando") # Entferne Rolando

cor_noRolando <- cor.test(dat_noRolando$Ages, dat_noRolando$Vision) # Korrelation
reg_noRolando <- lm(data=dat_noRolando, formula=Vision ~ Ages) # Regressionsmodell y = a + bx
```

### Vergleich Korrelationen
Zuerst vergleichen wir die Korrelationen mit und ohne Rolando:
<div class = "row"> <div class = "col-md-6">
```{r}
cor # mit Rolando
```
</div> <div class = "col-md-6">
```{r}
cor_noRolando # ohne Rolando
```
</div> </div>

Die Korrelation ist natürlich auch ohne Rolando negativ, sie hat aber einen erheblichen Sprung gemacht. Vorher lag sie bei `r round(cor$estimate, 2)` und ohne Rolando bei `r round(cor_noRolando$estimate, 2)`. Der p-Wert war in diesem Beispiel zwar schon vorher signifikant, doch auch dieser wurde im rechten Output deutlich kleiner. Das ist erstaunlich, da wir nur eine einzige Beobachtung von 28 entfernt haben. Es wird also klar welchen Einfluss einzelne Ausreißer haben können.

### Vergleich Regressionen
<div class = "row"> <div class = "col-md-6">
```{r}
summary(reg) # mit Rolando
```
</div> <div class = "col-md-6">
```{r}
summary(reg_noRolando) # ohne Rolando
```
</div> </div>

Beim Vergleich der geschätzten Regressionsmodelle zeigt sich auf den ersten Blick kein so deutliches Bild wie bei der Korrelation: Der Schätzwert für die Achsenabschnitte änderte sich von `r round(summary(reg)$coefficients[1,1],2)` zu  `r round(summary(reg_noRolando)$coefficients[1,1],2)` und der der Steigung von `r round(summary(reg)$coefficients[2,1],2)` zu  `r round(summary(reg_noRolando)$coefficients[2,1],2)`. Und tatsächlich muss man etwas genauer hinschauen wenn man die Geraden in Plots vergleicht:

<div class = "row"> <div class = "col-md-6">
```{r}
# Mit Rolando
ggplot(data=dat, 
       aes(x=Ages, y=Vision)) +
geom_point(size=3) + 
scale_x_continuous(name="Alter",       
                   limits=c(20, 60), 
                   breaks=seq(20, 60, by=10)) +
scale_y_continuous(name="Sehvermögen", 
                   limits=c( 0, 10), 
                   breaks=seq( 0, 10, by=1)) +
theme_classic() +
geom_smooth(method='lm', formula=y~x, se=F)
```
</div> <div class = "col-md-6">
```{r}
# Ohne Rolando
ggplot(data=dat_noRolando, 
       aes(x=Ages, y=Vision)) +
geom_point(size=3) + 
scale_x_continuous(name="Alter",       
                   limits=c(20, 60), 
                   breaks=seq(20, 60, by=10)) +
scale_y_continuous(name="Sehvermögen", 
                   limits=c( 0, 10), 
                   breaks=seq( 0, 10, by=1)) +
theme_classic() +
geom_smooth(method='lm', formula=y~x, se=F)
```
</div> </div>

Sicherlich ergibt es Sinn, dass die Gerade im Plot mit Rolando links etwas tiefer liegt - der Punkt von Rolando zieht sie quasi ein Stück weit zu sich. Insgesamt ist dieser Unterschied jedoch minimal. Es bleibt also die Frage wie jetzt zu bemessen ist, dass doch die Gerade ohne Rolando sehr viel besser zu der Gesamtheit aller Punkte passt, als die Gerade mit Rolando. 

### Bestimmtheitsmaß
Dafür gibt es das **Bestimmtheitsmaß** $R^2$ (auch Determinationskoeffizient), welches auch in der vorletzten Zeile des `summary(reg)` Outputs angegeben ist. Es gibt an, wie viel Streuung in den Daten durch ein vorliegendes lineares Regressionsmodell „erklärt“ werden kann. Demnach wird es in Prozent gemessen, kann also Werte zwischen 0% und 100% annehmen. Hier werden wir nun also fündig, dass auch bei der Regression der Rauswurf von Rolando zu einer erheblichen Verbesserung des Bestimmtheißtsmaßes geführt hat, da es mit Rolando nur bei `r round(summary(reg)$r.squared,3)*100` % lag und ohne Rolando auf `r round(summary(reg_noRolando)$r.squared,3)*100` % anstieg.

# Endergebnis
Schließlich wollen wir den Plot zum Abschlus ohne Rolando und mit allen Ergebnissen erstellen:

```{r}
ggplot(data=dat_noRolando, aes(x=Ages, y=Vision)) + # Definiere Daten
  geom_point(size=3, color="darkgrey") + # Scatter plot mit grauen Punkten der Größe 3
  scale_x_continuous(name="Alter",       limits=c(20, 60), breaks=seq(20, 60, by=10)) + # x-Achse
  scale_y_continuous(name="Sehvermögen", limits=c( 0, 10), breaks=seq( 0, 10, by=1)) + # y-Achse
  theme_classic() + # Simple, klassische Formatierung
  geom_smooth(method='lm', formula=y~x, se=FALSE, color="orange") + # Füge orange Gerade ein
  stat_regline_equation(aes(label=paste(..eq.label.., ..rr.label.., sep = "~`,`~~~")),
                        label.x=50, label.y=9) + # Füge Gleichung inkl. R2 ein
  stat_cor(method="pearson", label.x=40, label.y=3) + # Füge Korrelation mit p-Wert ein
  labs(caption="Zwei Personen wurden afgrund von fehlenden oder fehlerhaften Werten aus dem Datensatz entfernt")
```