21  Repetitionstutorat 14

# Wichtig: Ladet eure packages über den folgenden Befehl. Ladet keine weiteren packages.
pacman::p_load("sjPlot", "texreg", "haven", "readr", "readxl", "fixest", "MASS","glm.predict", "scales", "tidyverse")

22 Programm für heute: Eure Fragen aus dem Forum und den Tutoraten

  • Teil 1: Allgemeine Fragen
  • Teil 2: Repetition herausfordernde Aufgaben

23 Teil 1: Allgemeine Fragen

23.1 1.1 Wann muss man einen Seed setzen?

Einen Seed brauchen wir immer dann, wenn Zufallsprozesse involviert sind. Für die Tutorate trifft das auf folgende Funktionen zu:

?rbinom
?rnorm
?rbernoulli
?mvrnorm
?sample_n
?sample
?glm.predict::predicts

Setzt vor diese Funktionen jeweils einen beliebigen Seed, damit wir eure Resultate auf die Kommazahlen genau reproduzieren könnten.

  • Signalwörter:
  • Simulation/simuliere
  • Werte aus der Normalverteilung/aus einer anderem Verteilung ziehen
  • Stichprobe ziehen
  • Predicted probabilities (Vorhergesagte Wahrscheinlichkeiten berechnen)

23.2 1.2 Predicts() allgemein? Was sind die Unterschiede zwischen predict und predicts? Wann nutzt man den position-Parameter etc.?

Aufgabe: Stelle den Zusammenhang zwischen der Anzahl Opfer und der Wahrscheinlichkeit für die Todesstrafe graphisch dar.

df_penalty <- read_dta("death_penalty.dta")
df_penalty <- df_penalty %>%
  mutate(death = ifelse(death == 99, NA, death),
         ac = ifelse(ac == 99, NA, ac))


# Unser Regressionsmodell
mb <- glm(death ~ v2 + fv + ac + yv, 
              family = binomial("logit"),
              data = df_penalty)

# Predicted probabilities mit predict
## Werte für Kontrollvariable berechnen
table(df_penalty$fv) # Nominal -> Modus
median(df_penalty$ac) # Metrisch -> Median oder arithmetisches Mittel (mean)
table(df_penalty$yv) # Nominal -> Modus
# Ordinal -> Modus oder Median

p1 <- predict(mb,
newdata = data.frame(v2 = c(0,1,2,3,4,5,6), fv = c(0,0,0,0,0,0,0), 
                     ac = c(rep(0,7)), yv = c(rep(1,7))), type = "response")

p1

Nachteile von predict aus base-r

  • Einzusetzende Werte für Kontrollvariablen müssen bekannt sein
  • Mehr Schreibarbeit
  • Keine Konfidenzintervalle
  • Gibt vektor zurück statt dataframe

Aufgrund dieser Nachteile nutzen wir predicts aus dem package glm.predict.

# Predicted probabilities mit predicts()
set.seed(666)
preds2 <- predicts(mb, "0-6, 1; mode; median; mode")

preds2 %>%
  ggplot(aes(x = v2, y = mean,
             ymin = lower, ymax = upper)) +
  geom_ribbon(alpha = 0.3) +
  geom_line() +
  scale_y_continuous(labels = scales::percent)

## Niveau der Kontrollvariablen wichtig, da marginaler Effekt davon abhängt
set.seed(666)
preds <- predicts(mb, "0-6, 1; 1; 2; 1")

preds %>%
  ggplot(aes(x = v2, y = mean,
             ymin = lower, ymax = upper)) +
  geom_ribbon(alpha = 0.3) +
  geom_line() +
  scale_y_continuous(labels = scales::percent)

## Weiterer Vorteil von predicts: Möglichkeit Discrete Change zu berechnen
preds_dc <- predicts(mb, "0-6, 1; mode; median; mode", position = 1)
preds_dc

Predicts berechnet discrete changes als Wahrscheinlichkeit für Wert 1 - Wahrscheinlichkeit für Wert 2. Hier also Wechsel der Wahrscheinlichkeit, wenn wir statt 1, 0 Opfer haben.

Wir können die Interpretation auch umkehren:

preds_dc[1,7]*(-1)

Ausnahme: Falls ihr die Modellqualität eines binären Modells überprüfen sollt, nutzt ihr predict bzw. den Code von den Vorlesungsfolien

23.3 1.3 Recoding, group_by und summarize

# Vorbereitung
df_selects19 <- haven::read_dta("Selects2019.dta")
df_selects19[df_selects19 == ""] <- NA # NA's korrekt kodieren
df_selects19 <- df_selects19 %>% 
  mutate(
    leftright = as.numeric(case_when(
      leftright == "0 - Left" ~ "0",
      leftright == "10 - Right" ~ "10", 
      T ~ leftright))
  )

# Recoding
df_selects19 <- df_selects19 %>% 
  mutate(alter_cat = case_when(alter <= 34 ~ "18-34", # Wenn alter kleinergleich 34, dann speichere den Wert "18-34"
                             alter > 34 & alter < 55 ~ "35-54", # Wenn alter grösser 34 und kleiner 55, dann "35-54"
                             alter > 54 & alter < 75 ~ "55-74", # Wenn alter grösser 54 und kleiner 75, dann "55-74
                             alter >= 75 ~ "75+", # Wenn alter grössergleich 75 dann "75+"
                             TRUE ~ NA_character_)) # Wenn keine der Bedingungen zutrifft, speichere als NA. Hier NA_character_, weil die Variable als character gespeichert werden soll. Wenn numerisch, dann NA_real_, wenn logisch, dann NA_logical_


# Ziel: Durchschnitt der linksrechts-Einstellung pro Altersgruppen und Geschlechtskombination

# MUTATE: Neue Variable mu erstellen
agemeans <- df_selects19 %>% 
  group_by(alter_cat, gender) %>% 
  mutate(mu = mean(leftright, na.rm = TRUE))
# Wir haben so viele Zeilen wie Beobachtungen


# SUMMARIZE: Fasse jede Gruppe in einer Zeile zusammen
agemeans <- df_selects19 %>% 
  group_by(alter_cat, gender) %>% 
  summarize(mu = mean(leftright, na.rm = TRUE))
# Jetzt haben wir nur noch so viele Zeilen, wie Gruppen

23.4 1.4 Falls möglich, eine kurze Übersicht über die verschiedenen Modelle und wann man welches verwenden sollte + Did/ feols/ before-after

df_postage <- read_rds("postage.rds")

# Abhängige Variable ist intervallskaliert/kontinuierlich
# Hier: Abweichung Gemeindestimmbeteiligung von Schweizer Stimmbeteiligung

# Wir haben hier Panel-Daten, beobachten also Gemeinden über Zeit, vor und nach der Einführung von vorfrankierten Couverts.

# QUERSCHNITTSDESIGN: 
# Wir haben nur Daten zu einem Zeitpunkt, von Gemeinden mit und ohne vorfrankierte Couverts
quer <- df_postage |> 
  filter(vote_date == max(vote_date))

mod_quer <- lm(turnout_diff ~ postage,
                        data = quer)

# => Effekt ist kausal identifiziert, wenn es keine Confounder gibt


# BEFORE AND AFTER:
# Wir haben nur Daten von Gemeinden welche zu irgendeinem Zeitpunkt vorfrankierte Couverts eingeführt haben
ba_communes <- df_postage |>
  filter(postage == 1) |> 
  select(municipality) |>
  distinct()

ba <- df_postage |>
  filter(municipality %in% ba_communes$municipality)

mod_ba <- lm(turnout_diff ~ postage,
                        data = ba)

# => Effekt ist kausal identifiziert, wenn kein zeitlich variierendes Confounding vorliegt, z.B. keine zeitlichen Trends


# DIFFERENCE IN DIFFERENCES
# Wir haben Daten von zwei Zeitpunkten von Couvert-Gemeinden und Nicht-Couvert-Gemeinden
did <- df_postage %>%
  filter(time %in% c(-1,1))

# DiD-Modell schätzen.
mod_did <- lm(turnout_diff ~ factor(time)*group,
                data = did) 
## => Effekt ist kausal identifiziert, wenn parallel trends Annahme plausibel ist


# FIXED EFFECTS
# Wir haben Panel-Daten, mehrere Zeitpunkte, Couvert- und Nicht-Couvert-Gemeinden
mod_twoway <- feols(turnout_diff ~ postage | muninr + votedateid, 
                   #muninr = unit fixed effect, votedateid = time fixed effect
   data = df_postage)
## =>  Effekt ist kausal identifiziert, wenn es keine Confounder gibt, welche über Zeit und über Gemeinden variieren



# Alle Modelle im Vergleich
screenreg(list(mod_quer, mod_ba, mod_did, mod_twoway))







# Abhängige Variable ist binär/dichotom
# Hier: Todesstrafe verhängt (1) oder nicht (0)
logit <- glm(death ~ v2 + fv + ac + yv, 
              family = binomial,
              data = df_penalty)


probit <- glm(death ~ v2 + fv + ac + yv, 
              family = binomial("probit"),
              data = df_penalty)

# Wahl des Modells substantiell meistens egal
# => Effekt ist kausal, wenn es keine Confounder gibt

23.5 1.5 Wie ausführlich müssen die Interpretationen von Modellen sein?

Orientiert euch an den Beispielen in den Tutoratsunterlagen: - Besprecht sämtliche Koeffizienten. Intercept muss nicht unbedingt interpretiert werden (bei binärer Regression sowieso nicht!) - Bei linearer Regression: Interpretiert die signifikanten Koeffizienten inhaltlich: Wenn x um eine Einheit steigt, dann verändert sich y um XY Einheiten. Wenn ein Koeffizient nicht signifikant ist, könnt ihr ihn nicht auf diese Weise interpretieren, dann reicht “Variable X hat keinen statistisch signifikanten Effekt auf Y” - Bei binärer Regression: Interpretiert bei signifikanten Koeffizienten die Richtung des Effekts (positiv/negativ) - Interpretiert bei Interaktionen sämtliche Bestandteile, so wie in den Übungsblättern. - Achtet darauf die inhaltliche Frage zu beantworten, falls es eine gibt (bspw. wie beurteilst du Hypothese XY)

# Beispiel binärer Fall
mb <- glm(death ~ v2 + fv + ac + yv, 
              family = binomial("logit"),
              data = df_penalty)
summary(mb)

Die Variable v2 ist positiv und statistisch signifikant. Die Todesstrafe wird demnach mit der steigenden Anzahl von Opfern wahrscheinlicher. Dies entspricht auch unserer Intuition, wonach schwerere Verbrechen härter bestraft werden. Alle anderen Variablen weisen keinen statistisch signifikanten Zusammenhang mit der Wahrscheinlichkeit der Verhängung der Todesstrafe auf

# Geben wir nur den Interaktionsterm alter * gender an, werden die beiden
# Hauptterme automatisch hinzugefügt
selects_fit2 <- lm(leftright ~ alter * gender, data = df_selects19)
summary(selects_fit2)

# Alternativ könnten wir schreiben
#lm(lire ~ alter + gender + alter * gender, data = df_selects19)

  • (Intercept) gibt den durchschnittlichen Skalenwert für einen 0-jährigen Mann an. Die Interpretation dieses Wertes ist bedeutungslos, da im Datensatz keine 0-jährigen Personen enthalten sind (Extrapolation).

  • alter: Steigt das Alter um ein Jahr, steigt der durchschnittliche Skalenwert für Männer um ca. r round(coef(selects_fit2)[2],2) Einheiten, ceteris paribus.

  • genderFemale: Für eine 0-jährige Frau ist der durchschnittliche Skalenwert um ca. r abs(round(coef(selects_fit2)[3],2)) Skalenpunkte tiefer als für einen 0-jährigen Mann, ceteris paribus. Beachte, dass die Interpretation dieses Wertes aufgrund der Extrapolation bedeutungslos ist.

  • alter:genderFemale: Der Effekt von alter auf die Selbsteinschätzung ist für Frauen um rund r round(coef(selects_fit2)[4],3) Skalenpunkte stärker als für Männer.

24 Teil 2: Repetition herausfordernde Aufgaben

Abstimmung: Welche der folgenden herausfordernden Aufgaben sollen wir nochmal gemeinsam anschauen? - Übungsblatt 2: Simulation eines Experiments - Übungsblatt 3: Kontrafaktisches Szenario - Übungsblatt 5: R2 mit einer simulierten Variable - Übungsblatt 6: Confounding und das Problem von überschätzten/unterschätzten Koeffizienten - Übungsblatt 7: Demeaning - Übungsblatt 9: Visuelle Darstellung parallel trends - Übungsblatt 11: Konfidenzintervall mit For-Loop

24.1 ÜBUNGSBLATT 2: Aufgabe 5

Sei \(D_i\in\{0,1\}\) ein Indikator für den Treatment-Status einer Untersuchungseinheit \(i\). Für alle Einheiten in der Treatment-Gruppe gilt \(D_i=1\), während \(D_i=0\) für alle Kontrolleinheiten. Ferner sei \(Y_i(0)\) das potentielle Ergebnis unter Kontrollstatus, während \(Y_i(1)\) das potentielle Ergebnis unter Treatmentstatus darstellt.

  1. Simuliere die beobachteten Werte \(Y\) einer randomisierten kontrollierten Studie (engl. randomized controlled trial RCT) mit 1000 Teilnehmenden. Der Treatment-Effekt betrage 2.55 für sämtliche Individuen. Nehme an, dass die potentiellen Ergebnisse aus folgendem Datengenerierungsprozess stammen:

\[\begin{align*} Y(0) &= 1 + \varepsilon, & \varepsilon\sim\mathcal{N}(0,1) \\ Y(1) &= 1 + \text{TE} + \varepsilon, & \varepsilon\sim\mathcal{N}(0,1) \end{align*}\]

Lösungshilfe: Gehe Schritt für Schritt vor. (1) Setze zunächst einen seed = 32. (2) Erstelle mit tibble() einen Datensatz mit 1000 Reihen, welcher für jede Observation einen Fehlerterm, den Treatment-Effekt und den Treatment-Status enthält. Den Treatment-Status kannst du mit dem Code rbinom(n=1000, size=1, prob =0.5) randomisieren. (3) Generiere zwei Variablen für die potentiellen Ergebnisse, indem du die obenstehenden Formeln des Datengenerierungsprozesses in R nachbaust. (4) Nutze die Ergebnisse aus den zwei vorangegangenen Schritten um zum Schluss die beobachteten Werte \(Y\) zu generieren (4).

# TIDYVERSE VARIANTE
# Seed setzen für Nachvollziehbarkeit
set.seed(32)

# Treatment randomisieren, Fehlerterm und Treatment-Effekt generieren
exp <- tibble(
  treatment = rbinom(n = 1000, size = 1, prob = 0.5), # Randomisierung, wer Treatment erhält
  e = rnorm(n= 1000, mean = 0, sd = 1), # Individueller Fehler, Basisabweichung
  effect = 2.55 # Treatment-Effekt, hier für alle gleich, in der Realität aber häufig unterschiedlich (für diese Übung irrelevant)
)

# Alles in Formel Datengenerierungsprozesse einsetzen und potentielle Ergebnisse Y(1) und Y(0) generieren
exp <- exp %>%
  mutate(
    y_0 = 1 + e,
    y_1 = 1 + effect + e
  )

# Beobachtete Werte Y generieren
# Wenn treatment = 1, dann y_1, sonst y_0
exp <- exp %>%
  mutate(
    y = ifelse(treatment == 1, y_1, y_0)
  )
  1. Berechne nun sowohl den Mittelwert der beobachteten Werte für die Kontroll- als auch die Treatment-Gruppe und bilde die Differenz beider Werte. Wiederhole dies für beide potentiellen Ergebnisse. Vergleiche die beiden erhaltenen Differenzen.
# Beobachtete Werte (y)
# Tidyverse
mean_treated <- exp %>%
  filter(treatment == 1) %>%
  summarize(mean_y = mean(y))

mean_control <- exp %>%
  filter(treatment == 0) %>%
  summarize(mean_y = mean(y))

mean_treated - mean_control

# Mit base r
mean(exp$y[exp$treatment == 1]) - mean(exp$y[exp$treatment == 0])


# Potential outcomes (y_0 und y_1)
mean(exp$y_1 - exp$y_0)

Wir können erkennen, dass \(\mathbb{E}[Y|D=1]-\mathbb{E}[Y|D=0]\approx\mathbb{E}[Y(1)-Y(0)]\). In Worten ausgedrückt, sehen wir, dass der Unterschied in den beobachteten Mittelwerten approximativ dem durchschnittlichen Treatment-Effekt (engl. average treatment effect ATE) entspricht. Dies gilt, da das Treatment zufällig zugewiesen wurde, wodurch sich die potenziellen Outcomes in beiden Gruppen im Erwartungswert nicht unterscheiden.

Der beobachtete Mittelwert entspricht nicht genau dem Treatment-Effekt, da es aufgrund der randomisierten Treatment-Zuweisung in dieser spezifischen Iteration der Zuweisung mal diese, mal jene Observation ‘trifft’. Über viele Iterationen und viele Zuweisungen würde der Mittelwert der beobachteten Mittelwerte genau dem Treatment-Effekt entsprechen.

Übrigens: Den Treatment-Effekt, also in diesem Fall die Differenz zwischen den Mittelwerten der Treatment- und der Kontrollgruppe, können wir auch mittels einer linearen Regression berechnen. Dafür müssen wir den tatsächlich beobachteten Outcome auf den Treatment-Indikator regressieren. Das Resultat dieser einfachen Regression (ohne Kontrollvariablen) stellt jedoch ebenfalls nur dann den tatsächlichen Treatment-Effekt dar, wenn die Zuteilung zur Treatment-Gruppe randomisiert erfolgt ist.

m0 <- lm(y ~ treatment, data = exp)
coef(m0)

24.2 ÜBUNGSBLATT 3: Aufgabe 3c

  1. Importiere den Datensatz mean_turnout.csv und führe ihn mit dem bestehenden Datensatz CO2_Gesetz.csv zusammen.
#| message: false

# Importiere Datensätze
df_turn <- read_csv("mean_turnout.csv")
df_co2 <- read_csv("CO2_Gesetz.csv")

df_co2 <- df_co2 %>% 
  left_join(df_turn, by = c("canton_id", "canton_name", "mun_id", "mun_name"))

# Verbinden von Datensätzen über Schlüsselspalten
# left_join(): Alle Zeilen von links, passende von rechts
# right_join(): Alle Zeilen von rechts, passende von links
# inner_join(): Nur Zeilen, die in beiden vorkommen
# full_join(): Alle Zeilen aus beiden Tabellen

# Datensätze stumpf zusammenkleben
# bind_rows(): Tabellen untereinander stapeln (gleiche Spalten)
# bind_cols(): Tabellen nebeneinander kleben (gleiche Zeilenanzahl)
  1. (herausfordernd) Möglicherweise fragst du dich, ob das \(\text{CO}_2\)-Gesetz angenommen worden wäre wenn die Stimmbeteiligung in den ländlichen Gemeinden am 13. Juni nicht so aussergewöhnlich hoch gewesen wäre. Berechne die absolute Anzahl an Ja-Stimmen für jede Gemeinde in diesem kontrafaktischen Szenario anhand der durchschnittlichen Stimmbeteiligung der sechs vorangehenden Abstimmungsterminen (mean_turnout_pct). Nehme an, dass alle anderen Faktoren gleich bleiben (Anzahl Stimmberechtigte, Verhältnis von Ja- zu Nein-Stimmen etc.).
# Erstelle neuen Data Frame für kontrafaktisches Szenario
df_co2_kf <- df_co2 %>%
  select(canton_id, canton_name, mun_id, mun_name, anzahlStimmberechtigte, mean_turnout_pct, jaStimmenInProzent) %>% 
  mutate(stimmbeteiligung = mean_turnout_pct / 100,
         AnteilJaStimmen = jaStimmenInProzent / 100, 
         
         # Berechne gültige Stimmen und absolute Anzahl Ja-Stimmen
         gueltigeStimmen = stimmbeteiligung * anzahlStimmberechtigte,
         jaStimmenAbsolut = round(AnteilJaStimmen * gueltigeStimmen)
         )
  1. Wäre das \(\text{CO}_2\)-Gesetz in obigem Szenario angenommen worden?
# Tatsächlicher Anteil Ja-Stimmen
sum(df_co2$jaStimmenAbsolut) / sum(df_co2$gueltigeStimmen) * 100

# Anteil Ja-Stimmen in kontrafaktischem Szenario
sum(df_co2_kf$jaStimmenAbsolut) / sum(df_co2_kf$gueltigeStimmen) * 100

24.3 ÜBUNGSBLATT 6: Aufgabe 4

=> siehe PDF in Ordner

24.4 ÜBUNGSBLATT 5: Aufgabe 3

  1. Schätze nun eine multiple lineare Regression ohne Interaktion, welche leftright durch die beiden Prädiktoren alter und gender erklärt. Interpretiere den Koeffizienten für die Altersvariable.
selects_fit2 <- lm(leftright ~ alter + gender, data = df_selects19)
coef(selects_fit2)

Eine Erhöhung des Alters um ein Jahr ist, ceteris paribus, mit einer Erhöhung des durchschnittlichen Links-Rechts Skalenwerts um ca. 0.019 Einheiten assoziiert. Je älter eine Person ist, desto rechter ist sie.

  1. Was ist der \(R^2\)-Wert der in 3 a) geschätzten Regression? Wie interpretierst du diesen Wert? Tipp: Verwende dazu summary(fit), wobei fit dein geschätztes Modell darstellt.
summary(selects_fit2)
# mit $r.squared erhalten wir nur den R^2-Wert des summary-Outputs
summary(selects_fit2)$r.squared

Durch das geschätzte Modell können wir nur ca. 3 % der Varianz in der politischen Selbsteinschätzung erklären.

  1. (herausfordernd) Füge dem Selects-Datensatz nun eine simulierte Variable var_sim hinzu, welche normalverteilt ist mit Mittelwert 3 und Varianz 9. Verwende dazu einen seed von 1. Füge diese neue Variable anschliessend dem Regressionsmodell aus 3 a) hinzu. Verwende die Funktion screenreg() aus dem Package texreg und vergleiche die \(R^2\)- sowie die angepassten \(R^2\)-Werte zwischen dem hier und in 3 a) geschätzten Modellen. Beachte, dass wir die geschätzten Modelle mit list(fit1, fit2) zusammenfügen müssen, da screenreg() als erstes Argument eine Liste erwartet. Verwende das zusätzliche Argument digits = 5, um die Werte besser vergleichen zu können. Was fällt dir auf?
set.seed(1)
df_selects19 <- df_selects19 %>%
  mutate(var_sim = rnorm(n = nrow(df_selects19), mean = 3, sd = 3))


# beachte, dass die Länge des generierten Vektors der Anzahl Zeilen in
# df_selects19 entsprechen muss

selects_fit3 <- lm(leftright ~ alter + gender + var_sim, data = df_selects19)


screenreg(list(selects_fit2, selects_fit3), digits = 5)

# Hinweis für Forschungsarbeiten:
wordreg(list(selects_fit2, selects_fit3), digits = 5, file = "regtable.docx")

Wir können erkennen, dass der \(R^2\)-Wert für die Regression, welche die zufällig generierte Variable var_sim enthält, einen leicht höheren \(R^2\)-Wert aufweist, als die Regression ohne var_sim. Fügen wir einem Regressionsmodell mit \(p\) Prädiktoren \(k\) weitere Prädiktoren hinzu, wird die erklärte Varianz \(SSR=\sum_{i=1}^n(\hat{Y}_i-\hat{Y})^2\) des Modells mit \(p+k\) Prädiktoren niemals kleiner sein als \(SSR\) für das Modell mit \(p\) Prädiktoren. Daher wird der Determinationskoeffizient \(R^2\) immer ansteigen (oder gleich bleiben), wenn wir einem Modell weitere Prädiktoren hinzufügen, selbst wenn diese wie var_sim völlig zufällige Werte enthalten.

Im Gegensatz zum \(R^2\) hängt das angepasste \(R^2\) (engl. adjusted \(R^2\)) von der Anzahl Prädiktoren in einem Modell ab und bestraft das Hinzufügen von Prädiktoren, welche nicht mit der abhängigen Variable zusammenhängen. Dementsprechend sehen wir, dass das angepasste \(R^2\) für das Modell, welches var_sim als zusätzlichen Prädiktor enthält etwas geringer ist wie für das Modell ohne var_sim.

24.5 ÜBUNGSBLATT 7: Aufgabe 3

  1. Lies den resourcecurse.dta Datensatz ein und benutze wie in der Vorlesung institutions als abhängige und oil als unabhängige Variable.

Der unit fixed effects Schätzer wird manchmal auch within Schätzer genannt, da er nur Variation innerhalb einer Einheit berücksichtigt. Man kann diesen Schätzer finden indem man ein unit fixed effects Modell schätzt, aber auch durch Datenmanipulation vor dem Schätzen des Modells (anders gesagt: wir sollten auf zwei unterschiedliche Herangehensweisen zum selben Resultat kommen). Wir wollen uns nun anschauen, wie die zweite Variante funktioniert, um die Mechanik hinter fixed effects Modellen besser zu verstehen.

Die Grundidee ist einfach. Wir wollen nur Variation innerhalb der Einheiten analysieren und nicht zwischen Einheiten. Eine Möglichkeit besteht darin, die \(Y\) und die \(X\) Variable pro Land jeweils um den Mittelwert zu bereinigen. Wir können also bspw. für jedes Land \(i\) den Mittelwert von \(X_{it}\) bestimmen, den wir als \(\bar{X_i}\) bezeichnen können. Nun generieren wir ein neues \(X^*\), welches um diesen Mittelwert bereinigt ist, also \(X^*_{it} = X_{it} - \bar{X_i}\). Mit der \(Y\) Variable machen wir dasselbe. Führe diese Schritte für die Variablen institutions und oil in R aus und speichere das Resultat in einem neuen Datensatz df_within.

Tipp: Beachte das Datenformat mancher Variabeln und verwende group_by() und ungroup().

df_resource <- foreign::read.dta("resourcecurse.dta")



df_within <- df_resource %>%
  group_by(country) %>%
  mutate(oil = as.numeric(oil),
         Y_WithinAverage = mean(institutions, na.rm = T), 
         X_WithinAverage = mean(oil, na.rm = T)) %>%
  ungroup() %>% # Gruppierung auflösen
  mutate(Y_demeaned_within = institutions - Y_WithinAverage,
         X_demeaned_within = oil - X_WithinAverage)
df_within

Wir können nun die Varianz in der erklärenden Variable anschauen und sehen, dass die Varianz in der de-meaned \(X\) Variable kleiner ist, da wir nur noch die Varianz über Zeit innerhalb eines Landes behalten.

var(df_within$oil, na.rm = TRUE)
var(df_within$X_demeaned_within, na.rm = TRUE)
  1. Schätze nun ein unit-fixed-effects Modell mit den nicht-manipulierten Daten. Verwende institutions als Outcome und oil als erklärende Variable. Wie gross ist der Effekt von oil? Und ist das konsistent mit der Resource Curse Theorie?
unit_fe <- feols(institutions ~ oil | country, data = df_within)
coef(unit_fe)

Die Resultate dieses Modelles sind konsistent mit der Theorie, da mehr Öl zu schlechteren Institutionen führt. In der Vorlesung haben wir noch weitere Modelle gesehen und auch Argumente, weshalb obiges Modell wahrscheinlich nicht frei von Confounding ist.

  1. Schätze nun ein normales Regressionsmodell, aber verwende die neuen Variablen, die du in Aufgabe 3a erstellt hast. Was fällt dir auf, wenn du den Koeffizientenschätzer von Öl in diesem Modell mit dem Schätzer in 3b vergleichst?
result <- lm(Y_demeaned_within ~ X_demeaned_within, data = df_within)
coef(result)

Die beiden Schätzer sind identisch! Es zeigt sich also, dass wir entweder fixed effects Modelle schätzen können oder aber eine normale Regression, wenn wir zuerst die Daten bereinigen (de-meaning).

24.6 ÜBUNGSBLATT 9: Aufgabe 3b

Die Orginaldaten wurden so aufbereitet, dass sie für die Difference-in-Differences Methode angewendet werden können. Die Variable group klassifiziert die Gemeinden als Teil der Treatment oder Kontroll-Gruppe. Die Variable time zeigt den zeitlichen Verlauf an. Die Variable turnout_diff beschreibt die Abweichung der Stimmbeteiligung zur durchschnittlichen schweizweiten Stimmbeteiligung an diesem Tag. Diese nutzen wir als abhängige Variable. Zeige graphisch, ob die Modellannahme für die Difference-in-Differences Methode zutrifft.

#Möglichkeit 1: Werte händisch berechnen
df_postage %>%
  #gruppieren nach Zeit und Gruppe
  group_by(group, time) %>%
  #durchschnittliche Abweichung für jeden Zeitpunkt pro Gruppe berechnen
  summarize(mean_diff = mean(turnout_diff, na.rm = T)) %>%
  #gruppieren um beide linien darstellen zu können
  ggplot(aes(x = time, y = mean_diff,
             color = group)) +
  #vertikale Linie für Zeitpunkt 0
  geom_vline(xintercept = 0, linetype= "dashed") +
  #trend darstellen
  geom_point() +
  geom_smooth(se = F) +
  theme_bw()


#Möglichkeit 2: Werte ggplot berechnen lassen
df_postage %>%
  ggplot(aes(x = time, y = turnout_diff,
             color = group))+
  geom_vline(xintercept = 0, linetype= "dashed") +
  #Daten nach Gruppe und Zeit mit mean Funktion summieren
  #Als Linie darstellen
  stat_summary(fun = "mean", geom = "line") +
  theme_bw()

24.7 ÜBUNGSBLATT 11: Aufgabe 4 (herausfordernd)

Generiere mithilfe des untenstehenden Codes einen Datensatz.

Der Code erzeugt einen Datensatz mit 2000 Beobachtungen, welcher für diese Übung die Population/Grundgesamtheit bildet. Die Variable V1 agiert als unsere abhängige Variable und V2 als unsere unabhängige Variable. Mithilfe der Funktion sample_n() sollst du nun 100 Individuen zufällig aus dieser Population ziehen und damit eine lineare Regression rechnen. Speichere das obere und untere Ende des Konfidenzintervalles dieser Regression. Das Konfidenzintervall erhältst du mit der Funktion confint(). Nutze einen Loop und wiederhole den gesamten Prozess (Sample ziehen und Konfidenzintervall abspeichern) 1000 Mal.

Rechne anschliessend eine Regression mit der gesamten Population. Der resultierende Koeffizient dieser Regression ist somit unser unverzerrter wahrer Wert. Stelle deine Konfidenzintervalle diesem Wert gegenüber. Wie häufig war der wahre Wert nicht im Konfidenzintervall?

#Varianz-Kovarianz-Matrix
sigma <- rbind(c(1,-0.8), c(-0.8,1))

#Durchschnitt der Verteilung
mu <- c(10, 2) 

set.seed(666)
df <- as.data.frame(mvrnorm(n=2000, mu=mu, Sigma=sigma))
# Leere Vektoren erzeugen
coef <- rep(NA, 1000)
upper <- rep(NA, 1000)
lower <- rep(NA, 1000)

# Loop, der Sample zieht, Regression schätzt und Konfidenzintervall abspeichert.
for (i in 1:1000) {
  
  df_sample <- sample_n(df, size = 100)
  reg <- lm(V1 ~ V2, df_sample)
  coef[i] <- reg$coefficients[2]
  lower[i] <- confint(reg)[2,1]
  upper[i] <- confint(reg)[2,2]
  
}



# Wahren Wert berechnen
full_model <- lm(V1 ~ V2, df)

# Wie häufig war wahrer Wert nicht im Konfidenzintervall?
n_not <- tibble(upper = upper,
                  lower = lower,
                  true_val = full_model$coefficients[2]) %>%
  filter(upper < true_val | lower > true_val) %>% 
  nrow()
n_not

# Alternativer Code ohne dataframe
sum(upper < full_model$coefficients[2] | lower > full_model$coefficients[2])




# Illustration mittels plot
df_regs <- tibble(coef = coef,
                  lower = lower,
                  upper = upper,
                  ID = 1:1000) %>%
  mutate(interval = case_when( # Die Variable "interval" benötigen wir für die 
                               # Einfärbung der Intervalle in ggplot
    upper < full_model$coefficients[2] ~ "Out",
    lower > full_model$coefficients[2] ~ "Out",
    T ~ "In"
  ))

ggplot(df_regs, aes(y = reorder(ID, coef),
                    x = coef,
                    xmin = lower,
                    xmax = upper,
                    color = interval))+
  geom_vline(xintercept = full_model$coefficients[2],
             linetype = "longdash",
             color = "grey")+
  geom_pointrange(alpha= 0.1)+
  labs(y=NULL)+
  scale_color_manual(values=c("dodgerblue", "darkred"))+
  theme_classic()+
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        legend.position = "none")

In gut 95% aller Regressionen (r 1000 - n_not von 1000) war der wahre Wert innerhalb des Konfidenzintervalls zu finden. Dies entspricht der durch den p-Wert = 0.05 festgelegten akzeptablen Fehlerquote. Somit liegen wir bei 1000 Schätzungen lediglich in rund 50 Fällen (in diesem Beispiel sind es genau 41) ausserhalb unseres angegebenen Unsicherheitsbereiches.

?

Lernuhr
Uhrzeit
Lerntimer