---
title: "Tutorium 5"
format: html
editor: visual
---

## Setup

```{r}
if (!require("pacman")) install.packages("pacman"); library(pacman)
p_load(ggplot2)
```
Das übliche Setup. In diesem Tutorium benötigen wir nur das ggplot2 package.

## Loops

# Zahlen von 1 bis 10:

```{r}
for (i in 1:10) {
  print(i)
}
```

Wir nutzen eine einfache "for"-Schleife, um die Zahlen von 1 bis 10 auszugeben. Die Funktion "print" gibt die jeweilige Zahl aus.

# Fibonacci-Zahlen:

```{r}
# Fibonacci-Zahlen bis 100
fibonacci <- c(1, 1)
while (TRUE) {
  next_fib <- tail(fibonacci, 1) + tail(fibonacci, 2)[1]
  if (next_fib > 100) break
  fibonacci <- c(fibonacci, next_fib)
}
print(fibonacci)
```

Anschließend berechnen wir die Fibonacci-Zahlen bis 100. Dazu initialisieren wir die ersten beiden Zahlen der Fibonacci-Reihe und berechnen dann die nächsten Zahlen, indem wir die letzten beiden Zahlen addieren. Die Schleife wird beendet, wenn die nächste Zahl größer als 100 ist.

## Cars Datensatz

```{r}
data(mtcars)
head(mtcars)
str(mtcars)
summary(mtcars)
```

## Plotten des Cars Datensatzes

```{r}
# mpg gegen cyl
ggplot(mtcars, aes(x = cyl, y = mpg)) + geom_point() + labs(title = "mpg gegen cyl")

# mpg gegen disp
ggplot(mtcars, aes(x = disp, y = mpg)) + geom_point() + labs(title = "mpg gegen disp")

# mpg gegen hp
ggplot(mtcars, aes(x = hp, y = mpg)) + geom_point() + labs(title = "mpg gegen hp")

# mpg gegen drat
ggplot(mtcars, aes(x = drat, y = mpg)) + geom_point() + labs(title = "mpg gegen drat")

# mpg gegen wt
ggplot(mtcars, aes(x = wt, y = mpg)) + geom_point() + labs(title = "mpg gegen wt")
```

Mit Hilfe dieser Plots können wir uns einen ersten Eindruck über die vorliegenden Zusammenhänge machen. Dies können wir nutzen, um anschließend passende Modelle zu erstellen.


## Forward und Backward Selection

# Forward Selection "von Hand"

```{r}
# Wir starten mit einem leeren Modell:
model_0 <- lm(mpg ~ 1, data = mtcars)
model_1 <- lm(mpg ~ cyl, data = mtcars)
model_2 <- lm(mpg ~ disp, data = mtcars)
model_3 <- lm(mpg ~ hp, data = mtcars)
model_4 <- lm(mpg ~ drat, data = mtcars)
model_5 <- lm(mpg ~ wt, data = mtcars)
model_6 <- lm(mpg ~ qsec, data = mtcars)
model_7 <- lm(mpg ~ vs, data = mtcars)
model_8 <- lm(mpg ~ am, data = mtcars)
model_9 <- lm(mpg ~ gear, data = mtcars)
model_10 <- lm(mpg ~ carb, data = mtcars)
summary(model_0)
summary(model_1)
summary(model_2)
summary(model_3)
summary(model_4)
summary(model_5)
summary(model_6)
summary(model_7)
summary(model_8)
summary(model_9)
summary(model_10)
```

Wir suchen nach der Variable mit dem kleinsten p-Wert (sofern dieser kleiner als 0.05 ist). In unserem Fall ist dies die Variable "cyl". Wir fügen diese Variable in unser Modell ein und suchen dann die Variable mit dem kleinsten p-Wert, die noch nicht im Modell enthalten ist. Dies wiederholen wir, bis kein p-Wert mehr kleiner als 0.05 ist.

```{r}
model_1_1 <- lm(mpg ~ cyl + disp, data = mtcars)
model_1_2 <- lm(mpg ~ cyl + hp, data = mtcars)
model_1_3 <- lm(mpg ~ cyl + drat, data = mtcars)
model_1_4 <- lm(mpg ~ cyl + wt, data = mtcars)
model_1_5 <- lm(mpg ~ cyl + qsec, data = mtcars)
model_1_6 <- lm(mpg ~ cyl + vs, data = mtcars)
model_1_7 <- lm(mpg ~ cyl + am, data = mtcars)
model_1_8 <- lm(mpg ~ cyl + gear, data = mtcars)
model_1_9 <- lm(mpg ~ cyl + carb, data = mtcars)
summary(model_1_1)
summary(model_1_2)
summary(model_1_3)
summary(model_1_4)
summary(model_1_5)
summary(model_1_6)
summary(model_1_7)
summary(model_1_8)
summary(model_1_9)
```

Die Variable "wt" hat den kleinsten p-Wert (0.000222) und wird daher in das Modell aufgenommen.

```{r}
model_2_1 <- lm(mpg ~ cyl + wt + disp, data = mtcars)
model_2_2 <- lm(mpg ~ cyl + wt + hp, data = mtcars)
model_2_3 <- lm(mpg ~ cyl + wt + drat, data = mtcars)
model_2_4 <- lm(mpg ~ cyl + wt + qsec, data = mtcars)
model_2_5 <- lm(mpg ~ cyl + wt + vs, data = mtcars)
model_2_6 <- lm(mpg ~ cyl + wt + am, data = mtcars)
model_2_7 <- lm(mpg ~ cyl + wt + gear, data = mtcars)
model_2_8 <- lm(mpg ~ cyl + wt + carb, data = mtcars)
summary(model_2_1)
summary(model_2_2)
summary(model_2_3)
summary(model_2_4)
summary(model_2_5)
summary(model_2_6)
summary(model_2_7)
summary(model_2_8)
```

Nun erreicht keine weitere Variable unsere Schwelle von 0.05, sodass wir das Modell mit diesen beiden Variablen "abschließen".  

# Backward Selection "von Hand"

```{r}
# Wir starten mit einem kompletten Modell (10 Variablen):
model_10 <- lm(mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb, data = mtcars)
summary(model_10)
```

Wir überprüfen die Werte in der Splate "Pr(>|t|)":
Die Variable "cyl"  hat den höchsten Wert mit 0.9161 und wird daher ausgeschlossen.

```{r}
model_9 <- lm(mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb, data = mtcars)
summary(model_9)
```

Die Variable "vs" mit 0.8433 wird ausgeschlossen.

```{r}
model_8 <- lm(mpg ~ disp + hp + drat + wt + qsec + am + gear + carb, data = mtcars)
summary(model_8)
```

Die Variable "carv" mit 0.7470 wird ausgeschlossen.

```{r}
model_7 <- lm(mpg ~ disp + hp + drat + wt + qsec + am + gear, data = mtcars)
summary(model_7)
```

Die Variable "gear" mit 0.61964 wird ausgeschlossen.

```{r}
model_6 <- lm(mpg ~ disp + hp + wt + qsec + am + drat, data = mtcars)
summary(model_6)
```

Die Variable "drat" mit 0.46240 wird ausgeschlossen.

```{r}
model_5 <- lm(mpg ~ disp + hp + wt + qsec + am, data = mtcars)
summary(model_5)
```

Die Variable "disp" mit 0.29897 wird ausgeschlossen.

```{r}
model_4 <- lm(mpg ~ hp + wt + qsec + am, data = mtcars)
summary(model_4)
```

Die Variable "hp" mit 0.22309 wird ausgeschlossen.

```{r}
model_3 <- lm(mpg ~ wt + qsec + am, data = mtcars)
summary(model_3)
```

Nun sind alle signifikant "genug", sodass wir keine weitere Variable ausschließen.


## Forward Selection (mit einer Schleife)

```{r}
# Welche Eigenschaft wollen wir modellieren?
response <- "mpg"

# Alle anderen Spalten werden zu potenziell erklärenden Variablen
predictors <- setdiff(names(mtcars), response)

# Hilfsvariablen für eine erfolgreiche Selektion
current_predictors_fwd <- c() # Wir starten ohne vorausgewählte Prädiktoren
best_model_fwd <- lm(mpg ~ 1, data = mtcars) # Wir starten mit einem "leeren" Modell
best_pvalue_fwd <- 1 # Wir haben noch keinen besten p-Wert, also setzen wir diesen auf das Maximum "1"

# Wir definieren unsere Grenze, ab wann wir Variablen nicht mehr hinzufügen
threshold <- 0.05

# Unsere For-Schleife für die Vorwärtsselektion
for (i in 1:length(predictors)) {
  pvalues_fwd <- c()
  
  for (predictor in setdiff(predictors, current_predictors_fwd)) {
    formula_fwd <- as.formula(paste("mpg ~", paste(c(current_predictors_fwd, predictor), collapse = " + ")))
    model_fwd <- lm(formula_fwd, data = mtcars)
    pvalue <- summary(model_fwd)$coefficients[predictor, "Pr(>|t|)"]
    pvalues_fwd[predictor] <- pvalue
  }
  
  best_predictor_fwd <- names(which.min(pvalues_fwd))
  best_pvalue_fwd <- min(pvalues_fwd)
  
  if (best_pvalue_fwd < threshold) {
    current_predictors_fwd <- c(current_predictors_fwd, best_predictor_fwd)
    best_model_fwd <- lm(as.formula(paste("mpg ~", paste(current_predictors_fwd, collapse = " + "))), data = mtcars)
  } else {
    break # Wir brechen die Schleife ab, wenn der beste p-Wert über dem Schwellenwert liegt
  }
}

# Wir geben uns anschließend das (vermeintlich) beste Modell aus
cat("Forward Selection Model:\n")
summary(best_model_fwd)
```

## Backward Selection (mit einer Schleife)

```{r}
# Welche Eigenschaft wollen wir modellieren?
response <- "mpg"

# Alle anderen Spalten werden zu potenziell erklärenden Variablen
predictors <- setdiff(names(mtcars), response)

# Hilfsvariablen für eine erfolgreiche Selektion
current_predictors_bwd <- predictors
best_model_bwd <- lm(as.formula(paste("mpg ~", paste(current_predictors_bwd, collapse = " + "))), data = mtcars)

# Wir definieren unsere Grenze, ab wann wir Variablen nicht mehr hinzufügen
threshold <- 0.05

# Unsere While-Schleife für die Rückwärtsselektion
while (length(current_predictors_bwd) > 1) {
  pvalues_bwd <- summary(best_model_bwd)$coefficients[-1, "Pr(>|t|)"] # Wir schauen nicht auf den Intercept
  worst_predictor_bwd <- names(which.max(pvalues_bwd))
  worst_pvalue_bwd <- max(pvalues_bwd)
  
  if (worst_pvalue_bwd > threshold) {
    current_predictors_bwd <- setdiff(current_predictors_bwd, worst_predictor_bwd)
    best_model_bwd <- lm(as.formula(paste("mpg ~", paste(current_predictors_bwd, collapse = " + "))), data = mtcars)
  } else {
    break # Wir brechen die Schleife ab, wenn der beste p-Wert unter dem Schwellenwert liegt
  }
}

# Wir geben uns anschließend das (vermeintlich) beste Modell aus
cat("Backward Selection Model:\n")
summary(best_model_bwd)
```

Wir sehen, dass wir etwas unterschiedliche Modelle, mit unterschiedlichen Variablen erhalten. Für unsere Forward Selection werden uns nur die Variablen "wt" und "cyl" als signifikant und Bestandteil des Modells ausgegeben. Für die Backward Selection hingegen erhalten wir die Variablen "wt", "qsec" und "am" als signifikant und Bestandteil des Modells. Dies zeigt, dass die Auswahl der Variablen von der Methode abhängt, die wir verwenden. Zum Teil werden sogar völlig verschiedene Variablen ausgewählt. Es geht also nicht nur um die Frage "Wie viele?", sondern auch "Welche überhaupt?".


## AIC und BIC

# Berechnung "per Hand"

```{r}
# Wir betrachten hier unsere Modelle aus der Forward und Backward Selection:
model_fwd <- lm(mpg ~ wt + cyl, data = mtcars)
model_bwd <- lm(mpg ~ wt + qsec + am, data = mtcars)

# Wir schauen uns die logLikelihood an
loglik_fwd <- logLik(model_fwd)
loglik_bwd <- logLik(model_bwd)

# Wir verwenden die Formeln für AIC und BIC (aus der Vorlesung)
aic_hand_fwd <- -2 * as.numeric(loglik_fwd) + 2 * length(coef(model_fwd))
aic_hand_bwd <- -2 * as.numeric(loglik_bwd) + 2 * length(coef(model_bwd))

bic_hand_fwd <- -2 * as.numeric(loglik_fwd) + log(nrow(mtcars)) * length(coef(model_fwd))
bic_hand_bwd <- -2 * as.numeric(loglik_bwd) + log(nrow(mtcars)) * length(coef(model_bwd))

# Ausgabe im nächsten Block
```

# Zum Vergleich: Nutzen von R-Funktionen

```{r}
# Wir wenden die Funktionen AIC und BIC auf unsere Modelle an
aic_r_fwd <- AIC(model_fwd)
aic_r_bwd <- AIC(model_bwd)

bic_r_fwd <- BIC(model_fwd)
bic_r_bwd <- BIC(model_bwd)

# Wir geben uns die Ergebnisse aus um sie zu vergleichen
cat("Für unser Forward Selection Modell:\n")
cat("AIC per Hand:", aic_hand_fwd, "\n")
cat("AIC mit AIC Funktion:", aic_r_fwd, "\n")
cat("BIC per Hand:", bic_hand_fwd, "\n")
cat("BIC mit BIC Funktion:", bic_r_fwd, "\n")

cat("\nFür unser Backward Selection Modell:\n")
cat("AIC per Hand:", aic_hand_bwd, "\n")
cat("AIC mit AIC Funktion:", aic_r_bwd, "\n")
cat("BIC per Hand:", bic_hand_bwd, "\n")
cat("BIC mit BIC Funktion:", bic_r_bwd, "\n")
```

Und "von Hand" berechneten Werte sind etwas geringer als die Werte, welche R uns ausgibt. Allerdings ist die Struktur unverändert. AIC gibt uns etwas geringere Werte für unser Modell aus der Backward Selection aus, ebenso wie BIC. Dies deutet darauf hin, dass für unsere vorliegenden Daten die Backward Selection etwas "erfolgreicher" war, da wir niedrigere Werte anstreben.


## Vor- und Nachteile von Forward und Backward Selection

# Forward Selection
*Vorteile*: Beginnt mit einem einfachen Modell, daher risikoärmer für Overfitting; weniger rechenintensiv bei vielen Prädiktoren.
*Nachteile*: Kann wichtige Variablen übersehen, die in Kombination mit anderen Variablen signifikant wären.

# Backward Selection
*Vorteile*: Beginnt mit einem vollständigen Modell, daher umfassender; erkennt Kombinationseffekte von Variablen besser.
*Nachteile*: Höheres Risiko für Overfitting; rechenintensiver bei vielen Prädiktoren.

## AIC und BIC

# AIC (Akaike Information Criterion):
*Vorteile*: Flexibler, da es Modelle mit höherer Komplexität bevorzugt, wenn sie die Daten besser erklären.
*Nachteile*: Neigt dazu, komplexere Modelle zu wählen, was zu Overfitting führen kann.

# BIC (Bayesian Information Criterion):
*Vorteile*: Bestraft komplexere Modelle stärker als AIC, daher konservativer; bevorzugt einfachere Modelle, die die Daten gut erklären.
*Nachteile*: Kann zu einfachere Modelle bevorzugen, die wichtige Variablen ausschließen.

# Anwendung:
AIC und BIC sind sinnvoll, um Modelle zu vergleichen und Overfitting zu vermeiden.
AIC wird oft verwendet, wenn Vorhersagegenauigkeit im Vordergrund steht.
BIC wird bevorzugt, wenn das Ziel die Identifikation der wahren Modellstruktur ist.


