Robuste Statistik

# Robuste Statistik
### Statistik mit R <a href='https://therbootcamp.github.io'> The R Bootcamp </a> <a href='https://therbootcamp.github.io/SmR_2021Oct/'> </a>  <a href='https://therbootcamp.github.io'> </a>  <a href='mailto:therbootcamp@gmail.com'> </a>  <a href='https://www.linkedin.com/company/basel-r-bootcamp/'> </a>
### Oktober 2021

---

<div class="my-footer">
 
 
 <img src="https://raw.githubusercontent.com/therbootcamp/therbootcamp.github.io/master/_sessions/_image/by-sa.png" height=14 style="vertical-align: middle"/>
 
 <a href="https://therbootcamp.github.io/">
 
 
 www.therbootcamp.com
 
 
 </a>
 <a href="https://therbootcamp.github.io/">
 
 Statistik mit R | Oktober 2021
 
 </a>
 
 </div>

---

# Robuste Statistik

### Parametrische Annahmen

<ul>
 <li class="m1">Parametrische Modelle / Tests kommen mit <high>starken Annahmen</high>.</li>
 <li class="m2">Für <high>Effizienz und Verzerrungsfreiheit</high> sollten diese Annahmen erfüllt sein.</li>
</ul>

### Nicht-parametrische Statistik

<ul>
 <li class="m3">Nicht-parametrische Modelle kommen <high>ohne Annahmen über die Verteilung</high> in der Population aus.</li>
 <li class="m4">Involviert <high>Rangbildung</high>, <high>Kombinatorik</high> und <high>Bootstrap</high>.</li>
</ul>

]

]

---

# Annahmen der Regression

Auf Deutsch...

<ul>
 <li class="m1">Passen die Daten?</li>
 <ul>
 <li>Quantitativ oder dichotom (A1)</li>
 <li>Varianz grösser Null (A2)</li>
 <li>Interkorrelation moderat (A3)</li>
 </ul>
 <li class="m2">Passt das Modell?</li>
 <ul>
 <li>Linearer Zusammenhang (A4, A5, A7)</li>
 <li>Konstante Fehlervarianz (A6)</li>
 <li>Normalverteilte Fehler (A8)</li>
 </ul> 
</ul>

]

<img src="image/assumptions_berry.png" height=440px> 
from <a href="">Berry (1993)</a>

]

---

# Passen die Daten?

<ul>
 <li class="m1">Ein typisches Daten-Problem sind <high>zu hohe Prädiktorkorrelationen</high>.</li>
 
 <ul>
 <li><high>Singularität</high> - Mindestens ein Prädiktor kann <high>vollständig</high> durch andere erklärt werden.</li>
 
 <li><high>Multikollinearität</high> - Mindestens ein Prädiktor kann <high> fast vollständig</high> durch andere erklärt werden.</li>
 </ul>
 
 <li class="m2">Der <high>Variance Inflation Factor</high> (VIF) schätzt das Ausmass an Multikollinearität. Literatur empfiehlt VIF<10 <a href="https://therbootcamp.github.io/SwR_2019Apr/_sessions/RobustStats/literature/Stine1995.pdf">(Stine, 1995)</a>.</li>
</ul>

]

]
---

# Passt das Modell?

<ul>
 <li class="m1">Die Passung des Modells wird bestenfalls über <high>graphische Illustrationen</high> bestimmt.</li>
 
 <li class="m2"><high>Statistische Kennwerte</high> können zusätzliche Informationen liefern.</li>
</ul>

<table style="cellspacing:0; cellpadding:0; border:none; padding-top:10px">
 <col width="50%">
 <col width="50%">
<tr>
 <td bgcolor="white">
 Annahme
 </td>
 <td bgcolor="white">
 Statistik
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Linearität
 </td>
 <td bgcolor="white">
 <mono>lm</mono>, <mono>glm</mono> (curve fitting)
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Homoskedastizität
 </td>
 <td bgcolor="white">
 <mono>bartlett.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Normalität
 </td>
 <td bgcolor="white">
 <mono>skewness</mono>, <mono>kurtosis</mono>, <mono>shapiro.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Ausreisser
 </td>
 <td bgcolor="white">
 <mono>cooks.distance</mono>, <mono>mahalanobis</mono>, etc.
 </td> 
</tr>
</table>

]

Einfache Regression 
<img src="image/assumptions.png" height=420px style="padding-top:10px">

]

---

# Passt das Modell?

<table style="cellspacing:0; cellpadding:0; border:none; padding-top:10px">
 <col width="50%">
 <col width="50%">
<tr>
 <td bgcolor="white">
 Annahme
 </td>
 <td bgcolor="white">
 Statistik
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Linearity
 </td>
 <td bgcolor="white">
 <mono>lm</mono>, <mono>glm</mono> (curve fitting)
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Homoskedastizität
 </td>
 <td bgcolor="white">
 <mono>bartlett.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Normalität
 </td>
 <td bgcolor="white">
 <mono>skewness</mono>, <mono>kurtosis</mono>, <mono>shapiro.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Ausreisser
 </td>
 <td bgcolor="white">
 <mono>cooks.distance</mono>, <mono>mahalanobis</mono>, etc.
 </td> 
</tr>
</table>

]

Residualanalyse 
<img src="image/assumptions_res.png" height=420px style="padding-top:10px">

]

---

# Fehlende Prädiktoren

<br2>

<ul>
 <li class="m1"><high>Verletzungen der Linearitätsannahme</high> sind am gravierendsten.</li>
 
 <li class="m2"><high>Fehlende Prädiktoren</high> oft verantwortlich.</li>
 
 <li class="m3">Ob <high>zusätzliche Prädiktoren</high> die Passung verbessern kann getestet werden: <mono>anova(model_1, model_2)</mono>.</li>
</ul>

]

]

---

# Fehlende Prädiktoren

<ul>
 <li class="m1"><high>Test auf quadratischen Trend</high> in Szenario C.</li>
</ul>

```
## Analysis of Variance Table
## 
## Model 1: error_rate ~ stress
## Model 2: error_rate ~ stress + stress2
## Res.Df RSS Df Sum of Sq F Pr(>F) 
## 1 298 265 
## 2 297 112 1 153 408 <2e-16 ***
## ---
## Signif. codes: 
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.'
## 0.1 ' ' 1
```

]

]

---

# Nicht-parametrische Statistik

<ul>
 <li class="m1">Wenn Annahmen nicht zutreffen, kann auf nicht-parametrische Methoden zurückgegriffen werden, die höhrere <high>Robustheit</high> aufweisen</li>
</ul>

<table style="cellspacing:0; cellpadding:0; border:none; padding-top:10px">
 <col width="50%">
 <col width="50%">
<tr>
 <td bgcolor="white">
 Ansatz
 </td>
 <td bgcolor="white">
 Methode
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Rang-Kombinatorik
 </td>
 <td bgcolor="white">
 <mono>wilcox.test</mono>, <mono>friedman.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Häufigkeit-Kombinatorik
 </td>
 <td bgcolor="white">
 <mono>sign.test</mono>, <mono>chisq.test</mono>
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 M-Schätzung
 </td>
 <td bgcolor="white">
 <mono>rq</mono> (quantile regr.), <mono>rfit</mono> (rank regr.)
 </td> 
</tr>
<tr>
 <td bgcolor="white">
 Bootstrap
 </td>
 <td bgcolor="white">
 <mono>boot</mono> (anything)
 </td> 
</tr>
</table>

]

]

---

# Rangsummentest

<ul>
 <li class="m1">Gruppenvergleiche können mit dem <high>Rangsummentest</high> von Wilcoxon bzw. Mann-Whitney durchgeführt werden</li>
 <li class="m2">Der Rangsummentest testet, ob die <high>Summen gmeinsam gebildeter Ränge</high> in den Gruppen identisch sind.</li>
</ul>

]

]

---

# Wilcoxon test

]

```r
# data
group_1 <- c(18, 24, 29, 12, 11, 31)
group_2 <- c(27, 16, 23, 8, 15, 21)

# wilcoxon test
wilcox.test(group_1, group_2)
```

```
## 
## 	Wilcoxon rank sum exact test
## 
## data:  group_1 and group_2
## W = 22, p-value = 0.6
## alternative hypothesis: true location shift is not equal to 0
```

]

---

# Vorzeichentest

<ul>
 <li class="m1">Einfacher Test für abhängige Daten basierend auf der <high>Kombinatorik von Häufigkeiten</high>.</li>
 <li class="m2">Der Vorzeichentest testet, ob die <high>Vorzeichen von Differenzen (+,-)</high> gleich wahrscheinlich sind.</li>
</ul>

]

]

---

# Vorzeichentest

]

```r
# data
time_1 <- c(18, 24, 29, 12, 11, 31)
time_2 <- c(27, 16, 23, 8, 15, 21)

# sign test
N <- length(time_1)
n_plus <- sum(time_1 - time_2 > 0)
dbinom(x = n_plus : N, 
 size = N, 
 prob = .5) %>% # H0
 sum()
```

```
## [1] 0.3438
```
]

---

# Robuste Regression

<ul>
 <li class="m1">Die Robustheit von Regression kann verbessert durch alternative Fehlerfunktionen &rho;(e) verbessert werden.</li>
</ul>

<table style="cellspacing:0; cellpadding:0; border:none; padding-top:10px">
 <col width="13%">
 <col width="13%">
 <col width="35%">
 <col width="30%">
<tr>
 <td bgcolor="white">
 Funktion
 </td>
 <td bgcolor="white">
 Paket
 </td> 
 <td bgcolor="white">
 Fehlerfunktion
 </td>
 <td bgcolor="white">
 Beschreibung
 </td>
</tr>
<tr>
 <td bgcolor="white">
 <mono>lm</mono>
 </td>
 <td bgcolor="white">
 <mono>stats</mono>
 </td> 
 <td bgcolor="white">
 <img src="image/robust_lm.png" height=23px>
 </td>
 <td bgcolor="white">
 Kleinste-Quadrate
 </td>
</tr>
<tr>
 <td bgcolor="white">
 <mono>rq</mono>
 </td>
 <td bgcolor="white">
 <mono>quantreg</mono>
 </td> 
 <td bgcolor="white">
 <img src="image/robust_quantile.png" height=20px>
 </td>
 <td bgcolor="white">
 Quantil-Regression
 </td>
</tr>
<tr>
 <td bgcolor="white">
 <mono>rfit</mono>
 </td>
 <td bgcolor="white">
 <mono>Rfit</mono>
 </td> 
 <td bgcolor="white">
 <img src="image/robust_rank.png" height=40px>
 </td>
 <td bgcolor="white">
 Rangbasierte Regression
 </td>
</tr>
</table>
]

]

---

# Robuste Regression

<ul>
 <li class="m1">Die Robustheit von Regression kann verbessert durch alternative Fehlerfunktionen &rho;(e) verbessert werden.</li>
</ul>

```r
# Quantile regression
library(quantreg)
m <- rq(formula = y ~ x, 
 data = outlier_df)
summary(m)
```

```
## 
## Call: rq(formula = y ~ x, data = outlier_df)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd
## (Intercept)  0.01626     -0.04042
## x            0.47497      0.30777
##             upper bd
## (Intercept)  0.13365
## x            0.54799
```

]

---

# Robuste Regression

<ul>
 <li class="m1">Die Robustheit von Regression kann verbessert durch alternative Fehlerfunktionen &rho;(e) verbessert werden.</li>
</ul>

```r
# Rank-based regression
library(Rfit)
m <- rfit(formula = y ~ x, 
 data = outlier_df)
summary(m)
```

```
## Call:
## rfit.default(formula = y ~ x, data = outlier_df)
## 
## Coefficients:
##             Estimate Std. Error
## (Intercept)   0.0328     0.0551
## x             0.3048     0.0417
##             t.value p.value    
## (Intercept)    0.60    0.55    
## x              7.31 2.2e-12 ***
## ---
## Signif. codes:  
##   0 '***' 0.001 '**' 0.01 '*'
##   0.05 '.' 0.1 ' ' 1
## 
## Multiple R-squared (Robust): 0.05718 
## Reduction in Dispersion Test: 19.29 p-value: 2e-05
```

]

---

# Bootstraptests

<ul>
 <li class="m1">Bootstrap tests basieren auf wiederholten <high>Zufallsziehungen der Daten</high> (mit Zurücklegen).</li>
 
 <li class="m2">Für jedes Bootstrap sample wird die relevante Statistik berechnet, wodurch eine <high>verteilungsunabhängige, simulierte Stichprobenverteilung</high> entsteht.</li>
</ul>
]

<img src="image/muench.png" height=380px style="border-radius:0px"> 
 Baron Münchhausen, from <a href="https://de.wikipedia.org/wiki/Datei:M%C3%BCnchhausen-Sumpf-Hosemann.png">wikipedia.org</a>

]

---

# Bootstraptests

Schritt 1 Ziehe `\(R\)` <high>Bootstrap-Stichproben</high> `\(B_i\)` von Grösse `\(n\)` Beobachtungen (= eine Zeile im Datensatz).

Schritt 2 <high>Berechne Statistik</high> `\(f(B_i)\)`, z.B., einen Gruppenunterschied oder ein Regressionsgewicht, für jede der `\(R\)` Bootstrap-Stichproben.

Schritt 3 <high>Evaluiere die Simulierte Statistik</high>, z.B., mit einem Konfidenzinterval.

]

]

---

# Bootstraptests

Schritt 1 Ziehe `\(R\)` <high>Bootstrap-Stichproben</high> `\(B_i\)` von Grösse `\(n\)` Beobachtungen (= eine Zeile im Datensatz).

Schritt 2 <high>Berechne Statistik</high> `\(f(B_i)\)`, z.B., einen Gruppenunterschied oder ein Regressionsgewicht, für jede der `\(R\)` Bootstrap-Stichproben.

Schritt 3 <high>Evaluiere die Simulierte Statistik</high>, z.B., mit einem Konfidenzinterval.

]

```r
# Bootstrap
library(boot)

# bootstrap function
stat_fun <- function(data, indices){
 data <- data[indices,] # bootstrap
 m <- lm(error_rate ~ stress + stress2, 
 data = data)
 coefficients(m)
 }

# bootstrap samples
B <- boot(stress_df, 
 statistic = stat_fun, 
 R = 1000)
```

]

---

# Bootstraptests

Schritt 1 Ziehe `\(R\)` <high>Bootstrap-Stichproben</high> `\(B_i\)` von Grösse `\(n\)` Beobachtungen (= eine Zeile im Datensatz).

Schritt 2 <high>Berechne Statistik</high> `\(f(B_i)\)`, z.B., einen Gruppenunterschied oder ein Regressionsgewicht, für jede der `\(R\)` Bootstrap-Stichproben.

Schritt 3 <high>Evaluiere die Simulierte Statistik</high>, z.B., mit einem Konfidenzinterval.

]

```r
# Bootstrap CI for stress
boot.ci(B, index = 2)
```

```
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = B, index = 2)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.1534,  0.2883 )   ( 0.1478,  0.2806 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.1694,  0.3021 )   ( 0.1589,  0.2916 )  
## Calculations and Intervals on Original Scale
```

]

---

# Bootstraptests

Schritt 1 Ziehe `\(R\)` <high>Bootstrap-Stichproben</high> `\(B_i\)` von Grösse `\(n\)` Beobachtungen (= eine Zeile im Datensatz).

Schritt 2 <high>Berechne Statistik</high> `\(f(B_i)\)`, z.B., einen Gruppenunterschied oder ein Regressionsgewicht, für jede der `\(R\)` Bootstrap-Stichproben.

Schritt 3 <high>Evaluiere die Simulierte Statistik</high>, z.B., mit einem Konfidenzinterval.

]

```r
# Bootstrap CI for stress2
boot.ci(B, index = 3)
```

```
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = B, index = 3)
## 
## Intervals : 
## Level      Normal              Basic         
## 95%   ( 0.6626,  0.7902 )   ( 0.6653,  0.7920 )  
## 
## Level     Percentile            BCa          
## 95%   ( 0.6541,  0.7807 )   ( 0.6668,  0.7931 )  
## Calculations and Intervals on Original Scale
```

]

---

<h1><a href="https://therbootcamp.github.io/SwR_2021Oct/_sessions/RobustStats/RobustStats_practical.html">Practical</a></h1>