# Übung 4 {.unnumbered}

```{r}
#| echo: false
source("files/_common.R")
```

In der Übung beschäftigen wir uns mit Daten des Untergangs der Titanic. Schauen Sie sich jeweils die Befehle an und finden Sie heraus, was die Befehle machen. Machen Sie sich Notizen, wenn Sie etwas nicht verstehen.

## Daten einlesen {.unnumbered .unlisted .smaller .scrollable}

Legen Sie für diese Übung einen Ordner an, erstellen Sie dort eine qmd und kopieren Sie die Folgenden Daten in einen dortigen Unterordner «data»:

Suchen und nehmen Sie den Datensatz "train.csv": Downlaod der Daten: https://www.kaggle.com/competitions/titanic/data

```{r}
#| echo: true
#| include: true
#| code-fold: true
#| output: html

# DATEN_titanic <- read_csv("data/titanic/train.csv") # lese beim ersten Mal die Daten ein

DATEN_titanic <- readRDS("data/titanic/train.RDS") # Lese nach dem ersten Mal so die Daten ein.
saveRDS(DATEN_titanic, "data/titanic/train.RDS") # speichere die Daten
```

## Datenaufbereitung

```{r}
#| echo: true
#| include: true
#| code-fold: true
#| output: html
DATEN_titanic <- DATEN_titanic |> 
    mutate(Kinder = Age < 14) |> # Kinder werden als unter 14 Jährige definiert
    mutate(Age_z = Age - mean(Age, na.rm = TRUE)) |> # Das alter um den Mittelwert verschieben (zentrieren)
    mutate(Pclass_f = factor(Pclass, levels = c(3, 2, 1)), 
           Cabin_D = ifelse(is.na(Cabin), 0, 1))
```


... und mal die Daten angucken:

```{r}
#| echo: true
#| include: true
#| code-fold: true
#| output: html

DATEN_titanic |> # mache eine Zusammenfassung der Daten
  summary() 
```

Die PassengerId ist einefach eine Identifikationsnummer. 

- Es gibt dann eine Variable, die "Survived" heisst, die ein Minimum von 0 hat und ein Maximum von 1. Das deutet sehr auf eine Dummy hin. Da der Durchschnitt ("Mean") = 0.38 ist, wissen wir jetzt schon, dass 38 Prozent der Passagiere überlebt haben (der Mittelwert einer Dummy ist immer der Prozentsatz der 1er-Gruppe). 
- Dann kommt noch der Name als Zeichenvariable,
- das Alter, das von 0.42 bis 80 geht. Von 177 Personen fehlen die Altersangaben. 

Informieren Sie sich über die übrigen Variablen auf (https://www.kaggle.com/competitions/titanic/data)[«Kaggle»]. 

## Daten in Trainings- und Testdaten aufteilen {.unnumbered .unlisted .smaller}

Was passiert hier?

```{r}
#| echo: true
#| include: true
#| code-fold: true

# Setze eine Zufallszahl, damit die Ergebnisse replizierbar sind, also nicht jedes Mal eine neue Zufallszahl gesetzt wird und die Ergebnisse (bisschen) abweichen
set.seed(12345)

# Ziehe eine Zufallsstichprobe aus dem Filmdatensatz und bezeichne ihn als "train", also Trainingsdatensatz, mit dem die "Maschine" trainiert wird
train <- DATEN_titanic |> 
  sample_frac(.75)

# Bilde aus dem Rest der nicht für "train" gezogenen Fälle einen Test-Datensatz, indem nach 'id' die Fälle aus "Filme" das Gegenteil (anti) von zusammengetan (join) werden.
test  <- anti_join(DATEN_titanic, train, by = 'PassengerId')

```

Sehen kann man nach diesem r-Chunk übrigens nichts, weil nur Datensätze im Hintergrund aufgeteilt wurden. Also suchen wir mal nach guten Datenvisualisierungen.

**Übung**: Experimentieren Sie mit verschiedenen Zahlen in set.seed() und grösseren und kleineren Werten in sample_frac().

### Datenvisualisierung {.scrollable .smaller .unnumbered .unlisted }

Eine einfache Darstellungsmöglichkeit ist ein sogenannter "Scatterplot" der die Lage von Fällen in ein Koordinatensystem einteilt, das durch zwei Variablen gebildet wird. Im Beispiel (Abbildung \@ref(titanic-Datenvisualisierungen)\vpageref{titanic-Datenvisualisierungen}) ist es das Alter der Passagiere "Age" und der Fahrpreis "Fare". Als Zweites haben wir ein Balkendiagramm für die Passagierklassen. Im letzten sehr aufwendigen Plotgrafik werden die Zweierbeziehungen aller Variablen dargestellt, also wie sie miteinander korrelieren (obere Nebendiagonalen), wie ihre Verteilung ist (auf der Diagonale mit Namen) und wie ihre gemeinsame Streuung ist, also ein Scatterplot in der unteren Nebendiagonalen. Mehr zu diesen SPLOM finden Sie hier: https://cran.r-project.org/web/packages/psych/vignettes/intro.pdf.

```{r}
#| echo: true
#| code-fold: true
#| code-summary: "R-Code"

# Erstelle einen Scatterplot (Punktewolke)

DATEN_titanic |> 
  ggplot(aes(x = Age, y = Fare)) + 
  geom_point() 

# Erstelle ein Balkendiagramm
DATEN_titanic |> 
  ggplot(aes(x=Pclass)) +
  geom_bar()

# Alle auf einmal
psych::pairs.panels(DATEN_titanic) 
```


### Modellbildung für den Fahrpreis {.unnumbered .unlisted .scrollable .smaller}

**Übung:** Experimentieren Sie mit dem Syntax! Kopieren Sie sich die Zeile für das Modell, löschen von "Age_z" bis "Kinder" alles heraus und schätzen Sie mal. Schauen Sie sich das Ergebnis gut an und achten Sie darauf, was passiert, wenn Sie die Summanden für I(Age_z^2) usw. wieder in das Modell tun. Am Ende können Sie versuchen das Modell durch weitere Variablen ergänzen und verbessern oder andere Teile wieder herausnehmen. Manche Variablen wurden erst noch erstellt (zB "Kinder" oder "Age_z"). Die entsprechende Datenaufbereitung.qmd können Sie <a href="./data/Datenaufbereitung.qmd" downloaden > hier </a> abrufen .

**Übung:** Interpretieren Sie den Output!

```{r}
#| echo: true
#| code-fold: true

fit_titanic <- lm(log(Fare + 1) ~ Sex + Age_z + I(Age_z^2) + Survived + Pclass_f + Kinder, data = train)

broom::tidy(fit_titanic) |>
  mutate(across(where(is.numeric), ~round(.,3)))|>
  gt::gt()

broom::glance(fit_titanic)|> 
  round(2)|>
  gt::gt()
```


### Regressionsoutput mit `sjPlot::tab_model` {.scrollable .smaller .unnumbered .unlisted }

```{r}
#| echo: true
#| code-fold: true

sjPlot::tab_model(fit_titanic,
                  show.std = TRUE, # zeige die standardisierten Koeffizienten
                  show.est = TRUE, # zeige die unstandardisierten estimates
                  show.r2 = TRUE, # zeige R^2
                  show.fstat = FALSE, # zeige die F-Statistik
                  string.est = "b", # beschrifte die Estimator mit b
                  string.std = "std. b" # beschrifte die standardisierten b mit std. b
                  )
```

Als auch nicht ganz schlechte Alternative kann gtsummary::tbl_regression benutzt werden, was durch weitere gepipte Befehle angepasst und ergänzt werden kann. 

### Regressionsoutput mit `gtsummary::tbl_regression` {.scrollable .smaller .unnumbered .unlisted }

**Übung:** Was für ein Problem sehen Sie im folgenden Output? Bei welcher Befehlszeile müssten Sie das # wegnehmen, um das Problem zu lösen?

```{r}
#| echo: true
#| code-fold: true

## Alternative zum sjPlot-Ouptut. Sie können entscheiden, was Sie besser finden. 
fit_titanic |> 
  gtsummary::tbl_regression() |> 
  gtsummary::add_vif() |> 
#  gtsummary::add_global_p() |> # Gebe die Signifikanzen je Variable raus (Varianzaufklärung)
#  gtsummary::modify_header(estimate = "**b**")   |> # Sorgt dafür, dass die b's nicht Beta heissen, sondern b
  gtsummary::add_glance_source_note()
```

### Voraussetzungschecks {.unnumbered .unlisted}

**Übung:** Suchen Sie die Befehle und führen Sie sie aus:

Die möglichen Verletzungen der Voraussetzungen sind:

- Multikollinearität
- Deutliche Abweichung von der Normalverteilung in den Fehlern
- Heteroskedastizität
- Outlier

### Algorithmus für den Fahrpreis {.unnumbered .unlisted .scrollable .smaller}

Prognosewerte erstellen:

```{r}
#| echo: true
#| code-fold: true

preds <- predict(fit_titanic, test) 

modelEval <- cbind(log(test$Fare), preds) |> 
  as_tibble() 

# Beschrifte die Spalten
colnames(modelEval) <- c('actual', 'predicted')

plot <- ggplot(modelEval, aes(x = predicted, y = actual)) + 
  # Create a diagonal line:
  geom_abline(lty = 2) + 
  geom_point(alpha = 0.5) + 
  labs(y = "Tatsächliche Werte", x = "Vorhergesagte Werte") 

plot
```

### Überlebensprognose {.unnumbered .unlisted .scrollable .smaller}

Zweite Analyse: Überlebensprognose (Dummy) mit linearer Regression. Was denken Sie?

```{r}
#| echo: true
#| code-fold: true

DATEN_titanic <- readRDS("data/titanic/train.RDS")

model1 <- lm(Survived ~ Pclass_f + Sex + Age_z + I(Age_z^2) + Kinder, data=train)

broom::tidy(model1) |>
 mutate(p.value = scales::pvalue(p.value)) |> # damit der p-Wert nicht mit endlosen Nullen dargestellt wird
 mutate(across(where(is.numeric), ~round(.x, 3))) |> # über alle numerischen Variablen hinweg, runde sie auf 3 Stellen
 gt::gt()

broom::glance(model1)|>
 mutate(p.value = scales::pvalue(p.value)) |> # damit der p-Wert nicht mit endlosen Nullen dargestellt wird
 mutate(across(where(is.numeric), ~round(.x, 3))) |> # über alle numerischen Variablen hinweg, runde sie auf 3 Stellen
 gt::gt()
```

## Logistische Regression {.unnumbered .unlisted .scrollable .smaller}

```{r}
#| echo: true
#| code-fold: true

model2 <- glm(Survived ~ Pclass_f + Sex + Kinder,  family=binomial, data=train)

# Berechnung der Odds-Ratios (OR) für Pclass_f2 und Pclass_f1 um sie unten im Text verwenden zu können
OR_Pclass_f2 <- round(exp(summary(model2)$coefficients["Pclass_f2", 1]), 2) 
OR_Pclass_f1 <- round(exp(summary(model2)$coefficients["Pclass_f1", 1]), 2)

broom::tidy(model2)|>
 mutate(p.value = scales::pvalue(p.value)) |> # damit der p-Wert nicht mit endlosen Nullen dargestellt wird
 mutate(across(where(is.numeric), ~round(.x, 3))) |> # über alle numerischen Variablen hinweg, runde sie auf 3 Stellen
 gt::gt()

broom::glance(model2)|> 
 mutate(across(where(is.numeric), ~round(.x, 3))) |> # damit der p-Wert nicht mit endlosen Nullen dargestellt wird
  gt::gt()
```


### Kennwerte {.unnumbered .unlisted .scrollable .smaller}

**Übung:** Interpretieren Sie den Output

```{r}
#| echo: true
#| code-fold: true

model2 <- glm(Survived ~ Pclass_f + Sex +  Kinder,  family=binomial, data=train)

sjPlot::tab_model(model2,title = "Überlebensanalyse zum Titanicunglück mit sjPlot",
                  show.est = TRUE, # zeige die estimates
                  show.r2 = TRUE, # zeige R^2
                  show.fstat = TRUE # zeige die F-Statistik
                  )
```

### Wieder auch mit `gtsummary::` {.unnumbered .unlisted .scrollable .smaller}

```{r}
#| echo: true
#| code-fold: true

model2 |> 
  gtsummary::tbl_regression(exponentiate = TRUE) |> 
  gtsummary::add_vif() |> 
#  gtsummary::add_global_p() |> # damit könnte man die Signifikanz ganzer Faktoren testen, statt jede Ausprägung gegen die Referenz
  gtsummary::add_glance_source_note() |> 
  gtsummary::modify_caption("**Überlebensanalyse zum Titanicunglück mit gtsummary**")
```

## Voraussetzungschecks {.unnumbered .unlisted}

**Übung:** Führen Sie die Checks für die Voraussetzungen aus!

## Vorhersagetest {.unnumbered .unlisted .scrollable .smaller}

Was sagt dieser Test?

```{r}
#| echo: true
#| code-fold: true
#| eval: false
#| output: html

preds <- predict(model2, test) 

modelEval <- cbind(test$Survived, preds) |> 
  as.tibble() |> 
  mutate(preds = ifelse(preds > 0.5,1,0))

# Beschrifte die Spalten
colnames(modelEval) <- c('actual', 'predicted')

modelEval |> 
sjmisc::flat_table() 

modelEval |> 
  filter(!is.na(predicted)) |> 
  mutate(test = actual == predicted, na.rm = TRUE) |> 
  summarise(Accuracy = mean(test))
```