---
title: "Tutorial 7"
format: html
editor: visual
---

## Setup

```{r}
# Pakete laden, mit Hilfe von pacman
if (!require("pacman")) install.packages("pacman"); library("pacman")
p_load(data.table)
p_load(ggplot2)

# Daten laden
download.file("https://statistik.julianhinz.com/tutorials/tutorial7/avocado.csv",
              destfile = "../input/avocado.csv")
original_data <- fread("../input/avocado.csv")
```

Wir laden unseren Paketmanager `pacman`, und laden unseren Datensatz programmatisch herunter, und speichern ihn in dem dafür vorgesehen Ordner. Anschließend laden wir die Daten mit Hilfe von `fread` in R.


# Mit dem Datensatz arbeiten

```{r}
# Daten ansehen
head(original_data)
str(original_data)

# Variablen auswählen und vernünftig umbenennen
data = original_data[, .(date = Date,
                average_price = AveragePrice,
                total_volume = `Total Volume`,
                region = region,
                type = type)]
```

Das Umbennen der Variablen und die Selektion ist nicht zwingend nötig, erlaubt uns aber hier, mit unseren Namenskonventionen konstant zu bleiben, und unnötigen Ballast zu entfernen. Dies ist nur ratsam, wenn wir sicher sind, dass wir die anderen Variablen nicht mehr benötigen.


## Plotten

```{r}
# Variation über die Zeit
plot_data_time = data[, .(average_price = mean(average_price)), by = .(date, type)]
ggplot(plot_data_time) +
    theme_minimal() +
    geom_line(aes(x = date, y = average_price, color = type)) +
    scale_x_date(NULL) +
    labs(title = "Durchschnittspreis von Avocados über die Zeit",
        y = "Durchschnittspreis",
        color = "Typ") +
    theme(legend.position = "bottom")

## Variation über die Regionen
plot_data_region = data[, .(average_price = mean(average_price),
                     total_volume = mean(total_volume)), by = .(region, type)]
ggplot(plot_data_region) +
    theme_minimal() +
    geom_point(aes(x = total_volume,
                   y = average_price,
                   color = region)) +
    scale_x_log10("Total Volume") +
    scale_y_log10("Average Price") +
    labs(title = "Durchschnittspreis von Avocados nach Region") +
    theme(legend.position = "right")
```

Wir sehen klare Unterschiede bei den Durchschnittspreisen wenn wir nach Typ unterscheiden, und erkennen, dass Bio-Avocados immer teuer sind als klassische Avocados. Zudem sind die Verläufe insgesamt, absteigend und ansteigend, nahezu identisch. Auch die Regionen haben einen großen Einfluss auf den Preis und zeigen deutliche Unterschiede an, allerdings gibt es so viele Regionen, dass die weitere Interpretation des Plots schwierig ist.


## Lineare Regression

```{r}
# Verschiedene mögliche Gleichungen der linearen Regression
regression0 = lm(log(total_volume) ~ log(average_price), data = data)
summary(regression0)

regression1 = lm(log(total_volume) ~ log(average_price) + type, data = data)
summary(regression1)

regression2 = lm(log(total_volume) ~ log(average_price) + region, data = data)
summary(regression2)

regression3 = lm(log(total_volume) ~ log(average_price) + type + region, data = data)
summary(regression3)

regression4 = lm(log(total_volume) ~ log(average_price) * type, data = data)
summary(regression4)

```

Wie interpretieren wir das ganze?
Wir schauen uns `regression4`einmal genauer an:

Der Koeffizient für `log(average_price)` ist -1.43. Dies können wir so deuten, dass eine 1%ige Erhöhung des Durchschnittspreises mit einer Abnahme des Gesamtvolumens um etwa 1.43% einhergeht (da wir hier ein log-log-Modell haben).
Der Koeffizient für `typeorganic` ist -3.28. Hier liegt ein log-level-Modell vor, dies können wir so deuten, dass, wenn wir `log(average_price)` konstant halten, der Wert für `log(total_volume)` (also der Gesamtabsatz) für Bio-Produkte um 3.28*100% niedriger ist als für konventionelle Produkte.
Und der Interaktionsterm zwischen `log(average_price)` und `typeorganic` ist 0.28. Dies können wir so deuten, dass der negative Einfluss von steigenden Preisen auf den Absatz für Bio-Produkte um 0.28 Einheiten geringer ist als für konventionelle Produkte (also -1.43 + 0.28 = -1.15).


## Bootstrap

```{r}
# Seed setzen
set.seed(12345)

n <- 250
bootstrap_elastizität = numeric(n)

for (i in 1:n) {
    bootstrap_data <- data[sample(1:nrow(data), size = nrow(data), replace = TRUE)]
    bootstrap_model <- lm(log(total_volume) ~ log(average_price) + type + region,
                          data = bootstrap_data)
    bootstrap_elastizität[i] <- coef(bootstrap_model)[2]
    print(paste0(i, ": ", bootstrap_elastizität[i]))
}
```

Wir erhalten variierende, immer etwas abweichende Werte, welche jedoch alle um den vermeintlich wahren Wert, welchen wir in unserer Stichprobe finden, streuen. 


```{r}  
# Ergebnisse visualisieren
bootstrap_plot <- ggplot() +
    theme_minimal() +
    geom_histogram(aes(x = bootstrap_elastizität),
    binwidth = 0.0025, fill = "skyblue", color = "white") +
    labs(title = "Bootstrap-Schätzungen der Preiselastizität",
         x = "Preiselastizität", y = "Häufigkeit")
bootstrap_plot
```


## Bootstrap Konfidenzintervall

````{r}
# Konfidenzintervall berechnen
quantile(bootstrap_elastizität, probs = c(0.025, 0.975))

# Plot mit Konfidenzintervall und Schätzung
result_plot <- bootstrap_plot +
    geom_vline(xintercept = mean(bootstrap_elastizität),
                linetype = "dotted", color = "purple") +
    geom_vline(xintercept = quantile(bootstrap_elastizität, probs = c(0.025, 0.975)),
                linetype = "dashed", color = "red")
result_plot
```

Wir sehen, dass unser Konfidenzintervall um den wahren Wert streut, und dass die Schätzung in der Mitte des Konfidenzintervalls liegt. Dies ist ein gutes Zeichen, dass unser Bootstrap-Verfahren funktioniert hat.


## Vergleich mit traditionellem Konfidenzintervall

```{r}
# Konfidenzintervall mit normaler Verteilung
summary(regression3)
estimate <- summary(regression3)$coefficients[2, 1]
estimate_sd <- summary(regression3)$coefficients[2, 2]
ci_normal <- estimate + c(-1, 1) * qnorm(0.975) * estimate_sd

comparison_plot <- result_plot +
    geom_vline(xintercept = ci_normal,
                linetype = "dashed", color = "green")
comparison_plot
```

Wir sehen, dass wir mit Hilfe unseres Bootstraps eine gute Schätzung für das Konfidenzintervall erhalten haben, welche dem traditionellen Konfidenzintervall ähnlich ist. Dies ist ein gutes Zeichen, dass unser Bootstrap-Verfahren funktioniert hat. Das positive daran ist, dass unser Bootstrap völlig ohne Verteilungsannahme funktioniert, und lediglich auf unseren Daten basiert, was es zu einem sehr nützlichen Werkzeug für uns macht.