---
title: "Tutorial 8"
author: "Lars Abt"
format: pdf
editor: visual
---

## Differences-In-Differences

```{r}
# Pakete laden
pacman::p_load(data.table, ggplot2)

# Data.Table anlegen
did_data <- data.table(
    Gruppe = c("Treatment", "Treatment", "Control", "Control"),
    Zeit = c("Vorher", "Nacher", "Vorher", "Nacher"),
    Wert = c(70, 85, 75, 80)
    )

# Kopie der Daten erstellen, als Sicherungskopie (und für spätere Replizierbarkeit, da kein Originaldatensatz eingelesen werden kann)
did_data_backup <- copy(did_data)

# Daten visualisieren
did_plot <- ggplot(did_data, aes(x = Zeit, y = Wert, fill = Gruppe)) +
    geom_bar(stat = "identity", position = "dodge") +
    labs(title = "Differences-In-Differences Visualisierung", x = "Zeitpunkt", y = "Wert") +
    theme_minimal()
did_plot
```

Wir sehen, dass die Werte für beide Gruppen nach dem Treatment ansteigen, aber der Anstieg bei der Treatment-Gruppe deutlich stärker ist. Dies ist der Effekt, den wir mit einem Differences-In-Differences Schätzer messen wollen.

## Differences-In-Differences Schätzer

```{r}
# Wir nutzen unsere Sicherungskopie der Daten, sodass wir diesen Code-Chunk immer wieder ausführen können, ohne dass wir unsere Daten dauerhaft verändern und Replizierbarkeit verloren geht.
did_data <- copy(did_data_backup)

# Zeitwerte in numerische Werte übersetzen
did_data[, Zeit := ifelse(Zeit == "Vorher", 0, 1)]

# Gruppenwerte in numerische Werte übersetzen
did_data[, Gruppe := ifelse(Gruppe == "Treatment", 1, 0)]

# Berechnung des Effekts mit unserer "normalen" Formel über die Differenzen
diff1 <- did_data[Gruppe == 1 & Zeit == 1]$Wert - did_data[Gruppe == 1 & Zeit == 0]$Wert
diff2 <- did_data[Gruppe == 0 & Zeit == 1]$Wert - did_data[Gruppe == 0 & Zeit == 0]$Wert

diff_effect <- diff1 - diff2
diff_effect
```

Diff-In-Diff Effekt beträgt 10.


## Jetzt als Lineare Regression

```{r}
# Erstellen unseres Interaktionsterms
did_data[, Interaktion := Gruppe * Zeit]
# Wird nur zu 1, wenn wir in Treatment-Gruppe und nach Verabreichung des Treatments

# Lineares Regressionsmodell schätzen
fit_model <- lm(Wert ~ Gruppe + Zeit + Interaktion, data = did_data)
summary(fit_model)

# alternativ zu Interaktion: Gruppe * Zeit
fit_model <- lm(Wert ~ Gruppe + Zeit + Gruppe*Zeit, data = did_data)
summary(fit_model)
```

Wir sehen, dass unsere Diff-in-Diff Schätzung mit Hilfe eines linearen Regressionsmodells den gleichen Effekt schätzt, wie unsere manuelle Berechnung. Dieser ist in unserem Interaktionsterm wiederzufinden, da dieser genau den Fall des isolierten Treatment-Effekts darstellt. Zudem schätzen wir auch zeitliche Effekte in unserer `post`-Variable und Gruppeneffekte in unserer `group`Variable, welch hier aber keine Bedeutung besitzen.


### Bootstrapping

```{r}
# Working Directory überprüfen
# getwd()

# Daten herunterladen
download.file("https://statistik.julianhinz.com/tutorials/tutorial7/avocado.csv", destfile = "../input/avocado.csv")

# Daten einlesen
avocado_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(avocado_data)
str(avocado_data)

# Variablen auswählen und vernünftig umbenennen
data = avocado_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)

# Besonders relevant!
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)
    # type und region nicht zwingend notwendig. aber vlt. schöner
    bootstrap_elastizität[i] <- coef(bootstrap_model)[2]
    print(paste0(i, ": ", bootstrap_elastizität[i]))
    # nur für die Übersichtlichkeit
}
```

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
ci_normal

comparison_plot <- result_plot +
    geom_vline(xintercept = ci_normal,
                linetype = "dashed", color = "orange")
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.