--- title: "Empirische Studien in R" author: "Stefanie Radnik" date: "`r format(Sys.time(), '%d %B, %Y')`" output: html_document: theme: cosmo toc: yes toc_float: yes keep_md: yes pdf_document: toc: yes params: viridis_palette: viridis --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) knitr::opts_chunk$set(warning = FALSE, message = FALSE) options(dplyr.summarise.inform = FALSE) ``` In diesem R-Code wird das Paper „The Long-Term Impact of Miliary Service on Health: Evidence from World War II and Korean War Veterans" von Kelly Bedard und Olivier Deschênes behandelt. Dabei wird der langfristige gesundheitliche Einfluss des Militärdienstes anhand der Vorsterblichkeit der Veteranen betrachtet. Dafür werden auch die verschiedenen Todesursachen untersucht und speziell der Zusammenhang zwischen diesen und des heute noch im Militär üblichen Rauchens ermittelt. Die im Paper verwendeten Daten wurden zur Verfügung gestellt und hier genutzt. Die Kommentierung im Code dient lediglich der Leserlichkeit des Codes, daher wird am Anfang jeder Darstellung einmal kurz erläutert, was das Ziel der Darstellung ist und am Ende die Ergebnisse interpretiert. Hinweise: - Teilweise benötigt die Fehlerberechnung etwas mehr Laufzeit, also nicht wundern, wenn einzelne Code Zeilen etwas länger brauchen bis diese fertig durchgelaufen sind. - Einige Variablen werden hier mehrfach verwendet, also werden für die jeweiligen Tabellen überschrieben, daher ist es notwendig für die Ausgabe der Tabellen oder Figures immer den gesamten Code Chunk laufen zu lassen. Für die Tabelle 6 werden die Ergebnisse der Tabellen 2,3 und 5 verwendet, weshalb diese vorher generiert werden sollten. Bevor die einzelnen Tabellen generiert werden können, müssen natürlich aber erstmal die Packages und Daten geladen werden. Ich empfehle auch am Anfang erstmal alles laufen zu lassen und bei Bedarf die verschiedenen Tabellen/Abbildungen einzeln. - Jeder Code Chunk (außer die ersten beiden) enthalten jeweils eine Tabelle/Abbildung ```{r} #notwendige packages laden library (haven) library(ggplot2) library(dplyr) library(fastDummies) library(lmtest) library(sandwich) library(AER) library(gtsummary) library(kableExtra) library(lfe) library(tibble) ``` ```{r} #bereitgestellte Daten laden und in Variablen speichern data_67<- read_dta("Daten/cps_tobacco67.dta") data_90<- read_dta("Daten/cps_tobacco90s.dta") data_mcod_f<- read_dta("Daten/mcod_f.dta") data_mcod_m<- read_dta("Daten/mcod_m.dta") ``` # Tabelle 1 Die Tabelle 1 des Papers beinhaltet einen ersten Überblick der Todesursachen für Männer und Frauen, die zwischen 1920 und 1939 geboren sind und zwischen 40 und 75 Jahre alt geworden sind. Aufgrund der Einfachheit der Struktur der Tabelle, habe ich mich dazu entschlossen diese, als Balkendiagramm darzustellen, um so die Unterschiede zwischen Männern und Frauen auch optisch besser darzustellen. Die Angaben beziehen sich dabei immer auf Tote/1000 Menschen. ```{r Table 1, fig.dim=c(10,8)} # Todesursachenrate für Tote/1000 Maenner berechnen total_m = mean(1000*data_mcod_m$total) ischem_m = mean(1000*data_mcod_m$ischem) lungcan_m = mean(1000*data_mcod_m$lungcan) colon_m = mean(1000*data_mcod_m$colon) cerebro_m = mean(1000*data_mcod_m$cerebro) respir_m = mean(1000*data_mcod_m$chronic)+ mean(1000*data_mcod_m$pneum) acc_suic_m = mean(1000*data_mcod_m$accid)+mean(1000*data_mcod_m$suicide) # Todesursachenrate für Tote/1000 Frauen berechnen total_f = mean(1000*data_mcod_f$total) ischem_f = mean(1000*data_mcod_f$ischem) lungcan_f = mean(1000*data_mcod_f$lungcan) colon_f = mean(1000*data_mcod_f$colon) cerebro_f = mean(1000*data_mcod_f$cerebro) respir_f = mean(1000*data_mcod_f$chronic)+ mean(1000*data_mcod_f$pneum) acc_suic_f = mean(1000*data_mcod_f$accid)+mean(1000*data_mcod_f$suicide) #berechnen der sonstigen Todesursachen/1000 für Maenner und Frauen other_m = total_m - ischem_m - lungcan_m - colon_m - cerebro_m - respir_m - acc_suic_m other_f = total_f - ischem_f - lungcan_f - colon_f - cerebro_f - respir_f - acc_suic_f #Vektor mit allen Parametern erstellen c_tot = c(other_f, ischem_f, lungcan_f, colon_f, cerebro_f, respir_f, acc_suic_f, other_m, ischem_m, lungcan_m, colon_m, cerebro_m, respir_m, acc_suic_m) #Tabelle erstellen mit allen ermittelten Werten table1 = data.frame(c('other','Ischemic heart disease', 'Lung cancer', 'Colon cancer', 'Cerebrovascular disease', 'Respiratory diseases', 'Accidents and suicides'), c(rep("female", times = 7), rep("male", times = 7)), c_tot) #Spalten der Tabelle umbenennen colnames(table1) <- c('type', 'sex', 'value') #'fill'- Funktion bei Erstellung der Grafik soll nach Größe von 'value' sortiert sein, damit die häufigste Todesursache ganz unten ist und die seltenste Todesursache ganz oben, dazu wird die Variable order erstellt, um die gewuenschte Reihnfolge zu erhalten order <- table1 %>% group_by(type) %>% summarize(total = sum(value)) %>% arrange(total) %>% pull(type) #Grafik erstellen mit zwei Balken (male/female) auf der x-Achse und den akkumulierten Werten fuer die Todesursachen auf der y-Achse, farblich aufgeteilt nach den verschiedenen Todesursachen table1 %>% ggplot(mapping=aes(x = sex, y = value, fill = factor(type, levels=order), label = round(value,2)))+ geom_col(width = 0.5, color = "black", size = 0.2)+ geom_text(position = position_stack(vjust = 0.5))+ ylab('Deaths per 1000')+ scale_fill_discrete(name = "Types")+ theme_grey(base_size = 22) ``` Es ist deutlich zu erkennen, dass die Todesrate bei Männern fast doppelt so hoch ist, wie bei Frauen. Besonders markant ist der Unterschied bei den Todesursachen "Ischemic heart disease" und "Lung cancer", wo die Todesrate bei den Männern jeweils mehr als doppelt so hoch ist, wie bei den Frauen. # Abbildung 1 Vor Beginn der eigentlichen Analyse soll ein Eindruck des Zusammenhangs der Veteranenrate einer Kohorte und der mortality rate erlangt werden. Dazu wird ein Plot mit zwei Graphen und zwei y-Achsen erstellt. Der eine Graph enthält die Informationen zu der Veteranenrate für die Kohorten 1920–1939 und der andere Graph enthält die Residuen der logarithmierten mortality rate. Dafür wird eine lineare Regression durchgeführt, die als zu beschreibende Variable log(total), also den Logarithmus der gesamten Sterblichkeitsrate enthält und als beschreibende Variablen Dummy Variablen des Alters (40-75 Jahre). ```{r Figure 1, fig.dim=c(14,6)} #Figure 1 replizieren: hierfuer werden zwei Graphen benötigt (veteran rate/ residual log mortality rate), die schlussendlich in ein Figure geplottet werden #Werte fuer den Graph veteran rate ermitteln, dafuer werden die Daten, in denen nur die Maenner vorkommen verwendet. p1.1 = data_mcod_m %>% select(yob, vet) %>% #yob, vet Spalten werden ausgewaehlt group_by(yob) %>% #nach yob gruppiert summarise(vet = mean(vet)) #in vet Spalte wird Mittelwert von vet geschrieben #Dummy Variable fuer Variable "age" erstellen data_mcod_m <- dummy_cols(data_mcod_m, select_columns = 'age') #Spaltennamen getrennt mit "," in Zwischenspeicher kopieren, um diese anschliessend in "dummies" zu speichern writeClipboard(paste0('"',colnames(data_mcod_m),'"', collapse=", ")) dummies = c( "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48", "age_49", "age_50", "age_51", "age_52", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75") #dummy variablen mit "+" verbinden dummies.plus = paste0(dummies, collapse = '+') #Formel für den ols Schaetzer erstellen ols.formula = paste0("log(total) ~ ", dummies.plus) #lineare Regression laufen lassen und in Variable "ols" speichern ols = lm(ols.formula, data = data_mcod_m) #Mittelwert der Residuen je "yob" berechnen resid = data.frame(ols$residuals) %>% mutate(yob = data_mcod_m$yob) %>% group_by(yob) %>% summarise(mort_rate = mean(ols.residuals)) #Korrelation zwischen der mortality rate und der Veteranrate cor(resid$mort_rate, p1.1$vet) #Plot mit zwei y Achsen erstellen: plot.new() # Zusaetzlichen Platz fuer zweite y Achse erstellen par(mar = c(5, 4, 4, 4) + 0.3) #ersten Plot erstellen mit Punkten plot(resid$yob, resid$mort_rate, type = "p",pch = 22,ylim = c(-0.201, 0.152), xlim = c(1920, 1939.01), ylab = " ", xlab = "Year of Birth", lwd = 2, yaxs = "i", sub = "Age-adjusted log male mortality rate (Residuals from regression on age dummies) and male veteran rate, by year of birth") #gleichen Plot mit Linie erstellen lines(resid$yob, resid$mort_rate, lwd = 2) #Linien für das Gitternetz zeichnen abline(h = seq(-0.2, 0.15, 0.05), col = "gray", lty = 1, lwd = 1) #Achsenbeschriftung definieren axis(side = 1, seq(1920,1939,1)) axis(side = 2, seq(-0.2,0.15,0.05)) #zweiten Plot hinzufuegen ohne Achsen par(new = TRUE) plot(p1.1$yob, p1.1$vet, type = "p",pch = 4, axes = FALSE, xlab = "", ylab = "", lwd = 2, ylim = c(0, 0.901), yaxs = "i") lines(p1.1$yob, p1.1$vet, lwd = 2) # Zweite Achse hinzufuegen axis(side = 4, seq(0.0,0.9, 0.1)) mtext(" ", side = 4, line = 3) #Legende hizufuegen legend(x = "bottomleft", legend = c("Veteran Rate (right scale)", "Residual Log Mortality Rate (left scale)"), pch = c(4,22), lwd = c(2,2), bty = "n") ``` Es ist zu erkennen, dass es eindeutig einen positiven Zusammenhang zwischen der mortality rate und der Veteranenrate gibt. Es kann eine Korrelation von über 0.95 festgestellt werden. Besonders hervorzuheben sind dabei in der Grafik die Kohorten 1927–1929 und 1932–1935. Diese beiden Stellen markieren jeweils die Kohorten, die maßgeblich an dem Ende eines Krieges, und zwar des 2. Weltkrieges und des Koreakrieges, beteiligt waren. Allerdings wurden bei dieser Grafik Kohorteneffekte komplett ausgelassen. Einen Versuch einen Eindruck davon zu erhalten, bietet die Betrachtung der Frauen der verschiedenen Kohorten, da diese so gut wie nie einen Militärdienst absolvierten. # Abbildung 2 Der eben beschriebene Ansatz zur Erfassung von Kohorteneffekten durch die Betrachtung der Frauen wird nun in Figure 2 genutzt. Da davon ausgegangen wird, dass die Kohorteneffekte bei Frauen und Männern gleich sind, können diese Kohorteneffeke, die durch den Frauendatensatz ermittelt wurden von dem Datensatz der Männer abgezogen werden, um einen Eindruck zu erhalten, wie sich tatsächlich der Veteranenstatus auf die Sterblichkeit auswirkt. ```{r Figure 2, fig.dim=c(14,6)} #Figure 2 replizieren #age adjusted log mortality rate fuer Frauen berechnen: #Dummy Variablen fuer "age" von Frauen ertsellen data_mcod_f <- dummy_cols(data_mcod_f, select_columns = 'age') #Spaltennamen getrennt mit "," in Zwischenspeicher kopieren, um diese anschliessend in dummies zu speichern writeClipboard(paste0('"',colnames(data_mcod_f),'"', collapse=", ")) dummies = c( "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48", "age_49", "age_50", "age_51", "age_52", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75") #dummie variablen mit "+" verbinden dummies.plus = paste0(dummies, collapse = '+') #Formel für den ols schaetzer erstellen ols.formula = paste0("log(total) ~ ", dummies.plus) #lineare Regression laufen lassen und in Variable "ols" speichern ols_f = lm(ols.formula, data = data_mcod_f) #mittelwert der Residuen je "yob" berechnen resid_f = data.frame(ols_f$residuals) %>% mutate(yob = data_mcod_f$yob) %>% group_by(yob) %>% summarise(mort_rate = mean(ols_f.residuals)) #relative mort rate zwischen male und female berechnen (Differenzen in Differenzen Modell) resid_rel = data.frame(mort_rate = resid$mort_rate-resid_f$mort_rate) %>% mutate(yob =resid_f$yob) #Korrelation zwischen mortality rate und veteran rate berechnen cor(resid_rel$mort_rate, p1.1$vet) #Plot mit zwei y Achsen erstellen analog zu Figure 1 plot.new() # Zusaetzlichen Platz fuer zweite y Achse erstellen par(mar = c(5, 4, 4, 4) + 0.3) #ersten Plot mit Punkten erzeugen und Linie erzeugen plot(resid_rel$yob, resid_rel$mort_rate, type = "p",pch = 22,ylim = c(-0.062, 0.062), xlim = c(1920, 1939.01), ylab = " ", xlab = "Year of Birth", lwd = 2, yaxs = "i", bty = "o", sub = "Age-adjusted log male/female mortality rate (Residuals from regression on age dummies) and male veteran rate, by year of birth") lines(resid_rel$yob, resid_rel$mort_rate, lwd = 2) #Linien für das Gitternetz zeichnen abline(h = seq(-0.06, 0.06, 0.02), col = "gray", lty = 1, lwd = 1) #Achsenbeschriftung definieren axis(side = 1, seq(1920,1939,1)) axis(side = 2, seq(-0.06, 0.06,)) #zweiten Plot hinzufuegen ohne Achsen par(new = TRUE) plot(p1.1$yob, p1.1$vet, type = "p",pch = 4, axes = FALSE, xlab = "", ylab = "", lwd = 2, ylim = c(0, 0.901), yaxs = "i", bty = "o") lines(p1.1$yob, p1.1$vet, lwd = 2) # Zweite Achse hinzufuegen axis(side = 4, seq(0.0,0.9, 0.1)) mtext(" ", side = 4, line = 3) #Legende erstellen legend(x = "bottomleft", legend = c("Veteran Rate (right scale)", "Residual log male/female mortality rate (left scale)"), pch = c(4,22), lwd = c(2,2), bty = "n") ``` Auch die Figure 2 lässt, wie Figure 1, einen Zusammenhang zwischen der mortality rate und der Veteranrate vermuten. So gibt es wieder parallelen zwischen dem Ende des zweiten Weltkrieges und dem Koreakrieges und der damit verbundenen Reduktion der Veteran Rate und der Sterblichkeitsrate der betroffenen Kohorten. Allerdings hat hier sich die Korrelation der Veteranrate und der mortality rate auf ca. 0.85 reduziert. # Tabelle 2 Die eigentliche Analyse der Daten beginnt damit, den Veteraneffekt auf die Sterblichkeitsrate zu ermitteln. Dazu wird die Regression $$\log \bar M_{ct} = \alpha +\beta \bar V_c+ \delta (c)+ \lambda_{t-c}+ \varphi_t+u_{ct}$$ verwendet, wobei $\bar M_{ct}$ die Sterblichkeitsrate einer Kohorte c zu einer Zeit t ist, $\bar V_c$ die Veteran Rate einer Kohorte c, $\delta(c)$ eine lineare Funktion der Geburtsjahre und $\varphi_t$ unbeobachtete Kalenderjahreffekte sind. $u_{ct}$ ist hier der Störterm. In Tabelle 2 soll $\beta$ möglichst gut bestimmt werden, wobei hier nicht nur die gesamte Sterblichkeitsrate betrachtet wird, sondern auch nach den verschiedenen Todesursachen unterschieden wird. Die Tabelle 2 ist dabei in 5 Abschnitte unterteilt. Die ersten beiden Abschnitte stellen dabei eine leicht angepasste Vaiante von $\hat {\beta}$ dar, wobei $\hat {\beta}$ hier zur besseren Interpretierbarkeit noch durch dem Term $1000* \hat\beta \exp(\hat \mu_{ct})$ verändert wurden, wobei $\hat \mu_{ct}$ die vorhergesagten Werte der Regression sind. So lassen sich die Werte in Spalte 1 und 2 als Veteraneffekt/1000 interpretieren. Die Regressionen in Abschnitt (1) enthalten dabei nicht die Jahresdummys, in Abschnitt (2) dagegen schon. Der Abschnitt (3) enthält nochmal die Werte der Tabelle 1 und die Abschnitte (4) und (5) berechnen auf Basis des Abschnitts (2) und (3) die implizierten Sterblichkeitsraten für die Veteranen bzw. die Nicht-Veteranen. ```{r Table 2} #Tabelle 2 replizieren: #Variablen, die fuer die Regression genutzt werden reg_var = paste0("vet+yob+age_40+age_41+age_42+age_43+age_44+age_45+age_46+age_47+age_48+age_49+age_50+age_51+age_52+age_53+age_54+age_55+age_56+age_57+age_58+age_59+age_60+age_61+age_62+age_63+age_64+age_65+age_66+age_67+age_68+age_69+age_70+age_71+age_72+age_73+age_74+age_75", "") #Spalte 1 erstellen: Hier wird eine for-Schleife genutzt, die ueber die verschiedenen Todesursachen iteriert # Vektor mit den Todesursachen ueber den iteriert wird mort_cause = c("total","ischem","lungcan", "colon", "cerebro", "chronic+pneum", "accid+suicide") # Speicherort fuer den ermittelten Veteraneffekt und den dazugehoerigen Standardfehler vet_effect = c() sd_1 = c() # Spalten zu Daten hinzufuegen, in denen die Variablen "chronic" + "pneum" und "accid" + "suicide" addiert werden, da diese Todesursachen zusammengefasst betrachtet werden data_mcod_m = data_mcod_m %>% mutate("chronic+pneum" = chronic+pneum) %>% mutate("accid+suicide" = accid+suicide) # for Schleife iteriert über die verschiedenen Todesursachen for( i in mort_cause){ # Formel fuer die Regression wird erstellt ols_formula = paste0("log(", i ,") ~ ", reg_var) #Spalte mit aktueller Todesursache auswaehlen, diese wird fuer die Gewichtungen beonoetigt tmp = data_mcod_m %>%select(i) #Tabelle mit Gewichtungen erstellen weights_formula = data.frame(data_mcod_m$popmale*tmp/(1-tmp)) #erste Spalte des Vektors auswaehlen, da hier die Gewichtungen sind, so erhaelt man einen Vektor mit den entsprechenden Gewichtungen, der fuer die lineare Regression benoetigt wird vec = as.vector(weights_formula[,1]) #lineare Regression mit Gewichtungen wird berechnet ols <- lm(ols_formula, weights = vec, data = data_mcod_m) #relevante paramter extrahieren ols_summary = summary(ols) vet_coef = ols_summary$coefficients[,1][2] #Hier wird die die Variable "vet_coef", wie oben beschrieben, etwas angepasst, um eine bessere Interpretierbarkeit des Wertes zu erhalten exp_xb <- exp(predict(ols)) a <- vet_coef * exp_xb * 1000 vet <- mean(a) #Veteraneffekt fuer bestimmte Todesursache Vektor "vet_effect" hinzufuegen vet_effect = append(vet_effect, vet) #Standardabweichung berechnen: Hier werden die Fehler des cohort-level clustering verwendet, dafuer wird die Funktion "coeftest()" verwendet und die Standardabweichung von "vet" extrahiert sd_coef = coeftest(ols, vcov = vcovCL, cluster = ~yob)[2,2] #Standardabweichung, wird nun auch entsprechend dem Veteraneffekt angepasst b = sd_coef*sd_coef*exp_xb*exp_xb*1000*1000 c = mean(b) sd = sqrt(c) #Standardabweichung wird Vektor "sd_1" hinzugefuegt sd_1 = append(sd_1, sd) } # Die Fehler werden nun auf 2 Nachkommastellen gerundet sd_1 = round(sd_1,2) # Nun sollen die Fehler in einem Array gespeichert werden, wobei jeder Fehler mit 2 Klammern umschlossen wird also (), dies ist für die spätere Ausgabe wichtig sd_1_klammer = c() for (i in 0:(length(sd_1))){ sd_1_klammer[i] = paste0("(",sd_1[i],")") } # Spalte2 mit for-Schleife erstellen: Das Vorehen fuer die zweite Spalte ist das selbe, wie fuer die erste Spalte. Der einzige Unterschied ist lediglich, dass in der Regression nun auch Dummy Variablen fuer das Jahr mit einfließen #year dummies erstellen data_mcod_m <- dummy_cols(data_mcod_m, select_columns = 'year') #Spaltennamen mit "," Trennung in Zwischenablage kopieren writeClipboard(paste0('"',colnames(data_mcod_m),'"', collapse=", ")) #Jahresdummies wurden aus der Zwischenablage kopiert. Die restlichen Variablen wurden geloescht dummies_year = c("year_1968", "year_1969", "year_1970", "year_1971", "year_1972", "year_1973", "year_1974", "year_1975", "year_1976", "year_1977", "year_1978", "year_1979", "year_1980", "year_1981", "year_1982", "year_1983", "year_1984", "year_1985", "year_1986", "year_1987", "year_1988", "year_1989", "year_1990", "year_1991", "year_1992", "year_1993", "year_1994", "year_1995", "year_1996", "year_1997", "year_1998", "year_1999", "year_2000") #Variablen werden mit "+" zusammengefuegt dummies_year.plus = paste0(dummies_year, collapse = '+') #Jahresdummievariablen werden mit der Variablen, die auch in der Spalte 1 verwendet wurden, verknüft reg_var = paste0("vet+yob+age_40+age_41+age_42+age_43+age_44+age_45+age_46+age_47+age_48+age_49+age_50+age_51+age_52+age_53+age_54+age_55+age_56+age_57+age_58+age_59+age_60+age_61+age_62+age_63+age_64+age_65+age_66+age_67+age_68+age_69+age_70+age_71+age_72+age_73+age_74+age_75+", dummies_year.plus) # leerer Vektor als Speicherort fuer den Veteraneffekt wird erstellt vet_effect_year = c() # leerer Vektor als Speicherort fuer die zugehoerigen Standardabweichungen wird erstellt sd_2 = c() #ab hier ist es exact der selbe Code wie fuer die erste Spalte der Tabelle 2, daher wird hier nur Abschnittsweise kommentiert for( i in mort_cause){ #lineare Regression mit Gewichtungen ols_formula = paste0("log(", i ,") ~ ", reg_var) tmp = data_mcod_m %>%select(i) weights_formula = data.frame(data_mcod_m$popmale*tmp/(1-tmp)) vec = as.vector(weights_formula[,1]) ols <- lm(ols_formula, weights = vec, data = data_mcod_m) #relevante paramter extrahieren ols_summary = summary(ols) vet_coef = ols_summary$coefficients[,1][2] # Berechnung der mittleren Wirkung von vet und der Standardabweichung exp_xb <- exp(predict(ols)) a <- vet_coef * exp_xb * 1000 vet <- mean(a) vet_effect_year = append(vet_effect_year, vet) #Standardfehler berechnen sd_coef = coeftest(ols, vcov = vcovCL, cluster = ~yob)[2,2] b = sd_coef*sd_coef*exp_xb*exp_xb*1000*1000 c = mean(b) sd = sqrt(c) sd_2 = append(sd_2, sd) } # Die Fehler werden nun auf 2 Nachkommastellen gerundet sd_2 = round(sd_2,2) # Nun sollen die Fehler in einem Array gespeichert werden, wobei jeder Fehler mit 2 Klammern umschlossen wird also (), dies ist für die spätere Ausagbe wichtig sd_2_klammer = c() for (i in 0:(length(sd_2))){ sd_2_klammer[i] = paste0("(",sd_2[i],")") } #Spalte 3 erstellen # Hier wird wieder die Sterblichkeit durch verschiedene Ursachen berechnet. Dies sind die Zahlen, die auch schon in Tabelle 1 verwendet wurden total_m = mean(1000*data_mcod_m$total) ischem_m = mean(1000*data_mcod_m$ischem) lungcan_m = mean(1000*data_mcod_m$lungcan) colon_m = mean(1000*data_mcod_m$colon) cerebro_m = mean(1000*data_mcod_m$cerebro) respir_m = mean(1000*data_mcod_m$chronic)+ mean(1000*data_mcod_m$pneum) acc_suic_m = mean(1000*data_mcod_m$accid)+mean(1000*data_mcod_m$suicide) mean_rate = c(total_m, ischem_m, lungcan_m, colon_m, cerebro_m, respir_m, acc_suic_m) # Der Vektor mort_cause_table enthät die Namen der Todesursachen, die schlussendlich auch in der Tabelle erscheinen sollen mort_cause_table = c("All mortality causes", "Ischemic hert disease", "Lung cancer", "Colon cancer", "Cerebrovascular diseases", "Respiratory diseases", "Accidents and suicides") #Tabelle 2 wird erstellt tab2 = data.frame(mort_cause_table, veteran_effect = round(vet_effect,3), sd_1_klammer ,veteran_effect_year = round(vet_effect_year,3), sd_2_klammer, mean_rate = round(mean_rate,3)) #Implied mortality rate berechnen (x = implied mortality rate vets, mean_vet = mittlere vertanen rate, vet2 = veteran_effect_years): #Um dies zu berechenen gehe ich von der Formel mean_rate = mean_vet*x+(1-mean_vet)*(x-vet2) aus. Diese Formel wird nun einfach nach x umgestellt. Die jeweiligen zugehörogen Werte werden direkt aus tab2 entnommen #Die minimale Abweichung bei der implied mortality rate nonvet ist auf Rundungseffekte zurueckzufuehren # leeren Vektor als Speicherort fuer die implied mortality rate von Veteranen und Nicht-Veteranen bestimmen mort_rate_vet = c() mort_rate_nonvet = c() # Implied mortality mit dem oben beschriebenen Vorgehen fuer die verschiedenen Todesursachen berechnen mean_vet = mean(data_mcod_m$vet, na.rm = TRUE) for (i in 1:7){ vet = tab2[i,4] + tab2[i,6] - mean_vet * tab2[i,4] nonvet = vet - tab2[i,4] mort_rate_vet = append(mort_rate_vet, vet) mort_rate_nonvet = append(mort_rate_nonvet, nonvet) } # zu Tabelle 2 die implied mortality rates hinzufuegen tab2 = tab2 %>% mutate(impl_mort_rate_nonvet = round(mort_rate_nonvet,1), impl_mort_rate_vet = round(mort_rate_vet,1)) # Nun werden noch ein paar Zeilen hinzugefuegt zur Vervollstaendigung der Tabelle obs = c("Observations", 579, " ", 579, " " ," " ," ", " ", " " ) unre_age_dum = c("Unrestricted age dummies", "Yes", " ", "Yes", " ", " ", " ", " ") lin_cohort_trend = c("Linear cohort trend", "Yes", " ", "Yes", " ", " ", " ", " ") unre_year_dum = c("Unrestricted year dummies", "No", " ", "Yes", " ", " ", " ", " ") tab2 = rbind(tab2, obs, unre_age_dum, lin_cohort_trend, unre_year_dum) #Zur Markierung der Fussnote in der Tabelle footnote_header = paste0("Implied mortality rates" ,footnote_marker_alphabet(1, double_escape = T)) #Tabelle ordentlich darstellen mit Teilüberschriften und Fussnoten kable(tab2, booktabs = TRUE, escape = FALSE, col.names = c(" ","( 1 )"," ","( 2 )"," ","( 3 )","( 4 )","( 5 )"), align = c("lrlclccc"), caption = "Table 2 - Impact of Veteran Status on Mortality, Male-Only Sample") %>% add_header_above(c(" " = 1, "Veteran effect" = 4, "Mean rate" = 1, "Nonvets" = 1, "Vets" = 1)) %>% add_header_above(c(" " = 6, setNames(2, footnote_header)), escape = F) %>% kable_styling() %>% footnote(general = "The estimated standard errors in parentheses allow for cohort-level clustering. Each regression is weighted by the inverse of the sampling variance of the dependent variable.", alphabet = "Implied mortality rates based on the estimates in column 2", general_title = "Notes: ", footnote_as_chunk = T) ``` Durch die Werte in der ersten und zweiten Spalte der Tabelle 2 ist zu erkennen, dass es durchaus einen signifikanten Einfluss des Veteranstatus auf die Sterblichkeitsrate gibt, so ist diese bezogen auf alle Todesursachen bei Veteranen um 3,526 Tote/1000 höher als bei Nichtveteran. Betrachtet man zusätzlich noch die Jahreseffekte durch die Jahresdummies, so kommt man auf einen Wert von 3,424 Tote/1000, die durch den Veteranenstatus dazu kommen. Betrachtet man nun noch die Werte in Spalte 3, also die durchschnittliche Sterblichkeitsrate/1000, lassen sich die impliziten Sterblichkeitsraten für Veteranen und Nicht-Veteranen bestimmen, diese sind in Spalte 4 bzw. 5 zu finden. Diese beziehen sich jeweils auf die ermittelten Werte in Spalte 2. Also die Differenz der Werte in Spalte 5 und 6 ist genau der Wert in Spalte 2. Besonders hervorzuheben sind dabei die großen Unterschiede zwischen Veteranen und Nicht-Veteranen bei den Todesursachen "Ischemic heart disease" und "Lung cancer". # Tabelle 3 Um Kohorteneffekte der Sterblichkeitsraten besser abschätzen zu können, wird ein Differenzen-in Differenzen Modell herangezogen, welches durch die Formel $$\log \bar M_{sct} = \alpha +\beta \bar V_{sc}+ \omega_s + \delta_c+ \lambda_{st-c}+ \varphi_t+u_{sct}$$ beschrieben werden kann, s bezeichnet hierbei sex, also das Geschlecht der Person und $\omega_s$ ist eine Dummy Variable für das Geschlecht, die bei Männern 1 und bei Frauen 0 ist. Da Frauen so gut wie keinen Militärdienst geleistet haben, ist bei diesen $\bar V_{sc} = 0$. Die weiteren Variablen sind identisch zu dem Modell, welches für Tabelle 2 verwendet wurde. ```{r Table 3} #Tabelle 3 erstellen: Da das Vorgehen hier sehr aehnlich zu dem Vorgen in Tabelle 2 ist, werde ich hier besonders die Unterschiede kommentieren, wie das Vorbereiten der Daten. #year-dummy variablen fuer female Daten erstellen data_mcod_f <- dummy_cols(data_mcod_f, select_columns = 'year') # da hier die male und female Daten zusammengefuehrt werden, muessen diese vorerst auf eine Laenge gebracht werden. Ausserdem sollten beide Tabellen die selben Spaltennamen besitzen, dazu wird die Spalte "popmale" bzw. "popfemale" in "pop" umbenannt (I) und eine Dummyvariable "male" erstellt, die beinhalten, ob eine Person "male" ist (male = 1) oder female (male = 0) (II). Die Spalten "chronic+pneum" und "accid+suicide" aus dem male datensatz werden geloescht (III) und bei dem female Datensatz wird die "vet" Spalte hinzugefuegt, die es bei dem male Datensatz schon gibt, wo alle Werte auf null gesetzt werden (IV), da nur eine statistisch nicht relavente Gruppe an Frauen ueberhaupt einen Miliaerdienst geleistet hat data_m_3 = data_mcod_m %>% rename(pop = popmale) %>% #(I) mutate(male = 1) %>% #(II) select(-"chronic+pneum") %>% #(III) select(-"accid+suicide") #(III) data_f_3= data_mcod_f %>% rename(pop = popfemale) %>% #(I) mutate(male = 0) %>% #(II) mutate(vet = 0) #(IV) # Hier werden die bearbeiteten male und female Daten zusammengefuegt data_3 = rbind(data_m_3, data_f_3) %>% dummy_cols(select_columns = 'yob') %>% #yob dummies erstellen mutate("chronic+pneum" = chronic+pneum) %>% #Spalte "chonic+pneum" wird wieder hinzugefuegt, wo die ensprechenden Werte addiert werden mutate("accid+suicide" = accid+suicide) #Spalte "accid+suicide" wird wieder hinzugefuegt, wo die ensprechenden Werte addiert werden #Spaltennamen werden in Zwischenablage kopiert writeClipboard(paste0('"',colnames(data_3),'"', collapse=", ")) #yob dummies werden rauskopiert dummies_yob = c("yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939") #dummies_yob werden mit + getrennt dummies_yob.plus = paste0(dummies_yob, collapse = '+') #Spalte 1 #Regressionvariablen werden erstellt, wobei bei den Dummyvariabeln jeweils beruecksichtig wird, ob es sich um eine male oder female Person handelt reg_var = paste0("vet+yob+age_40*male+age_41*male+age_42*male+age_43*male+age_44*male+age_45*male+age_46*male+age_47*male+age_48*male+age_49*male+age_50*male+age_51*male+age_52*male+age_53*male+age_54*male+age_55*male+age_56*male+age_57*male+age_58*male+age_59*male+age_60*male+age_61*male+age_62*male+age_63*male+age_64*male+age_65*male+age_66*male+age_67*male+age_68*male+age_69*male+age_70*male+age_71*male+age_72*male+age_73*male+age_74*male+age_75*male+", dummies_yob.plus) # ab hier selbes Vorgehen wie bei Tabelle 2 vet3_effect = c() sd_1 = c() mort_cause = c("total","ischem","lungcan", "colon", "cerebro", "chronic+pneum", "accid+suicide") for( i in mort_cause){ ols_formula = paste0("log(", i ,") ~ ", reg_var) tmp = data_3 %>%select(i) weights_formula = data.frame(data_3$pop*tmp/(1-tmp)) vec = as.vector(weights_formula[,1]) #Ols Regression ausfuehren ols <- lm(ols_formula, weights = vec, data = data_3) #relevante paramter extrahieren und bearbeiten ols_summary = summary(ols) vet_coef = ols_summary$coefficients[,1][2] exp_xb <- exp(predict(ols)) a <- vet_coef * exp_xb * 1000 vet <- mean(a) vet3_effect = append(vet3_effect, vet) #Standardfehler berchnen sd_coef = coeftest(ols, vcov = vcovCL, cluster = ~yob)[2,2] b = sd_coef*sd_coef*exp_xb*exp_xb*1000*1000 c = mean(b) sd = sqrt(c) sd_1 = append(sd_1, sd) } # Die Fehler werden nun auf 2 Nachkommastellen gerundet sd_1 = round(sd_1,2) # Nun sollen die Fehler in einem Array gespeichert werden, wobei jeder Fehler mit 2 Klammern umschlossen wird also (), dies ist für die spätere Ausagbe wichtig sd_1_klammer = c() for (i in 0:(length(sd_2))){ sd_1_klammer[i] = paste0("(",sd_1[i],")") } #Spalte 2: Wird analog zu Spalte 1 erstellt, nur werden hier noch die year dummy Variablen hinzugefuegt #Regressionsvariablen mit year-dummys reg_var = paste0("vet+yob+age_40*male+age_41*male+age_42*male+age_43*male+age_44*male+age_45*male+age_46*male+age_47*male+age_48*male+age_49*male+age_50*male+age_51*male+age_52*male+age_53*male+age_54*male+age_55*male+age_56*male+age_57*male+age_58*male+age_59*male+age_60*male+age_61*male+age_62*male+age_63*male+age_64*male+age_65*male+age_66*male+age_67*male+age_68*male+age_69*male+age_70*male+age_71*male+age_72*male+age_73*male+age_74*male+age_75*male+", dummies_yob.plus, "+", dummies_year.plus) #ab hier wieder selbes Vorgehen wie vorher vet3_effect_year = c() sd_2 = c() mort_cause = c("total","ischem","lungcan", "colon", "cerebro", "chronic+pneum", "accid+suicide") for( i in mort_cause){ ols_formula = paste0("log(", i ,") ~ ", reg_var) ols_formula tmp = data_3 %>%select(i) weights_formula = data.frame(data_3$pop*tmp/(1-tmp)) vec = as.vector(weights_formula[,1]) #Ols Regression ausfuehren ols <- lm(ols_formula, weights = vec, data = data_3) #relevante paramter extrahieren und bearbeiten ols_summary = summary(ols) vet_coef = ols_summary$coefficients[,1][2] exp_xb <- exp(predict(ols)) a <- vet_coef * exp_xb * 1000 vet <- mean(a) vet3_effect_year = append(vet3_effect_year, vet) #Standardfehler berechnen sd_coef = coeftest(ols, vcov = vcovCL, cluster = ~yob)[2,2] b = sd_coef*sd_coef*exp_xb*exp_xb*1000*1000 c = mean(b) sd = sqrt(c) sd_2 = append(sd_2, sd) } # Die Fehler werden nun auf 2 Nachkommastellen gerundet sd_2 = round(sd_2,2) # Nun sollen die Fehler in einem Array gespeichert werden, wobei jeder Fehler mit 2 Klammern umschlossen wird also (), dies ist für die spätere Ausagbe wichtig sd_2_klammer = c() for (i in 0:(length(sd_2))){ sd_2_klammer[i] = paste0("(",sd_2[i],")") } #Tabelle 3 erstellen mit Spalten 1,2,3 tab3 = data.frame(mort_cause_table, veteran_effect = round(vet3_effect,3),sd_1_klammer, veteran_effect_year = round(vet3_effect_year,3), sd_2_klammer, mean_rate = round(mean_rate,3)) # Fuer die Spalten 4,5 wird wieder die implied mortality rate berechnet. Dies wird genau so gemacht, wie bei Tabelle 2 mean_vet = mean(data_mcod_m$vet, na.rm = TRUE) mort_rate_vet = c() mort_rate_nonvet = c() for (i in 1:7){ vet = tab3[i,4] + tab3[i,6] - mean_vet * tab3[i,4] nonvet = vet - tab3[i,4] mort_rate_vet = append(mort_rate_vet, vet) mort_rate_nonvet = append(mort_rate_nonvet, nonvet) } #implied mortality wird "tab3" hinzugefuegt tab3 = tab3 %>% mutate(impl_mort_rate_nonvet = round(mort_rate_nonvet,1), impl_mort_rate_vet = round(mort_rate_vet,1)) # Nun werden noch ein paar Zeilen hinzugefuegt zur Vervollstaendigung der Tabelle obs = c("Observations", 1158, " ", 1158, " " ," " ," ", " ", " " ) unre_s_age_dum = c("Unrestricted sex*age dummies", "Yes", " ", "Yes", " ", " ", " ", " ") unre_coho_dum = c("Unrestricted cohort dummies", "Yes", " ", "Yes", " ", " ", " ", " ") unre_year_dum = c("Unrestricted year dummies", "No", " ", "Yes", " ", " ", " ", " ") tab3 = rbind(tab3, obs, unre_s_age_dum, unre_coho_dum, unre_year_dum) #Markierung der Fussnote in der Tabelle footnote_header = paste0("Implied mortality rates" ,footnote_marker_alphabet(1, double_escape = T)) # Tabelle 3 ordentlich darstellen mit Teilueberschriften und Fussnoten kable(tab3, booktabs = TRUE,escape = FALSE, col.names = c(" ","( 1 )"," ","( 2 )"," ","( 3 )","( 4 )","( 5 )"), align = c("lrlclccc"), caption = "Table 3 - Impact of Veteran Status on Mortality, male and Female Sample") %>% add_header_above(c(" " = 1, "Veteran effect" = 4, "Mean rate" = 1, "Nonvets" = 1, "Vets" = 1)) %>% add_header_above(c(" " = 6, setNames(2, footnote_header)), escape = F) %>% kable_styling() %>% footnote(general = "The estimated standard errors in parentheses allow for cohort-level clustering. Each regression is weighted by the inverse of the sampling variance of the dependent variable.", alphabet = "Implied mortality rates based on the estimates in column 2", general_title = "Notes: ", footnote_as_chunk = T) ``` Die Ergebnisse, die durch das Differenzen in Differenzen Modell erhalten wurden, sind ähnlich zu den Ergebnissen des Modells aus Tabelle 2, wobei zu erwähnen ist das besonders bei den Todesursachen "Ischemic heart disease", "Lung cancer" und "Respiratory diseases" dem Veteranenstatus ein noch größerer Einfluss auf die Sterblichkeitsrate zu geschrieben wird. Der Wechsel des Vorzeichens bei den Todesursachen "Colon cancer", "Cerebrovascular diseases" und "Accidents and suicides" ist auf Grund des sowieso eher geringen Einflusses des Veteranstatus zu vernachlässigen. # Tabelle 5 Im nächsten Schritt der Analyse soll herausgefunden werden, welchen Einfluss der Veteranenstatus auf das Rauchen hat. Dazu soll folgende Regression ausgeführt werden: $$Smoke_{ict} = \alpha +\beta V_{ic}+\textbf{X}_{ict}\delta+\lambda_{t-c}+u_{ict},$$ dabei ist $Smoke_{ict}$ eine Dummy Variable, die 1 ist, wenn Person i in Kohorte c zur Zeit t ein Raucher war, $\textbf{X}_{ict}$ ist ein Vektor, der persönliche Eigenschaften, wie z.B. Bildung, Familienstand enthält, $\lambda_{t-c}$ gibt die unbeobachteten Alterseffekte an und $u_{ict}$ ist der Störterm. Die Regressionen werden für zwei Datensätze gemacht, einmal für ein Datensatz aus den Jahren 1967/68 und einmal aus dem Jahr 1990. Da es bei einem einfachen OLS-Schätzer, auf Grund der positiven Selektion beim Militärdienst, aber mit hoher Wahrscheinlichkeit zu einem verzehrten Ergebnis kommt, kann anstelle des OLS Schätzers auch ein Instrumentalvariablenschätzer verwendet werden, der das Korrelationsproblem $cor(V_{ic}, u_{ict}) \neq 0$ beheben soll. In Tabelle 5 werden zwei unterschiedliche Varianten des Instrumentalvariablenschätzers angewendet. Bei der ersten Variante (TSLS: men and women), werden die Daten der Männer und der Frauen verwendet und als Instrument für den Veteranenstatus werden Geburtsjahr dummies*male dummies verwendet. Bei der zweiten Variante des Instrumentalvariablenschätzers (TSLS: men only) werden nur die Daten für die Männer genommen. Bei der Regression werden außerdem die Geburtsjahrvariablen weggelassen und diese stattdessen als Instrument für den Veteranenstatus verwendet. ```{r Table 5} #Tabelle 5 erstellen: #OLS: men only #1967/68 #Daten vorbereiten: Dummy Variablen von "division" und "age" erstellen und anschliessend maennliche Personen herausfiltern und in "data_67_men" speichern data_67 = data_67 %>% dummy_cols(select_columns = 'division') %>% dummy_cols(select_columns = 'age') data_67_men = data_67 %>% filter(sex == 1) # Spaltennamen herausfinden und in dummies5 speichern, anschließend die benoetigten Variablen herausfiltern writeClipboard(paste0('"',colnames(data_67),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_28", "age_29", "age_30", "age_31", "age_32", "age_33", "age_34", "age_35", "age_36", "age_37", "age_38", "age_39", "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48") #Variablen mit + verbinden dummies5.plus = paste0(dummies5, collapse = '+') # Formel für den OLS Schätzer erstellen ols_formula = paste0("smoke100 ~ ", dummies5.plus) #ols ausführen und relevanten Parameter speichern ols = lm(ols_formula, data = data_67_men) ols_5_1 = coef(ols)[2] #Fehler berechnen sd_5_1 = coeftest(ols, vcov = vcovHC(ols, type = "HC3"))[2,2] #1990 #selbes vorgehen, wie vorher bei den Daten von 1967/90 #Daten vorbereiten data_90 = data_90 %>% dummy_cols(select_columns = 'division') %>% dummy_cols(select_columns = 'age') data_90_men = data_90 %>% filter(sex == 1) #Regressionvariablen erstellen writeClipboard(paste0('"',colnames(data_90),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75") dummies5.plus = paste0(dummies5, collapse = '+') #ols ausfueren und Fehler berechnen ols_formula = paste0("smoke100 ~ ", dummies5.plus) ols = lm(ols_formula, data = data_90_men) ols_5_2 = coef(ols)[2] sd_5_2 = coeftest(ols, vcov = vcovHC(ols, type = "HC3"))[2,2] #TSLS: men and women #Instrumentalvariablenschaetzer #1967/1968 #Daten vorbereiten -> yob, sex dummy erstellen data_67 = data_67 %>% dummy_cols(select_columns = 'yob') %>% dummy_cols(select_columns = "sex") #Binaere Variabeln fuer Instrumentalvariablenschaetzung ertstellen, die beinhalten ob Person in einem bestimmen Jahr geboren ist und maennlich ist myb20=(data_67$yob==1920)*(data_67$sex==1); myb21=(data_67$yob==1921)*(data_67$sex==1); myb22=(data_67$yob==1922)*(data_67$sex==1); myb23=(data_67$yob==1923)*(data_67$sex==1); myb24=(data_67$yob==1924)*(data_67$sex==1); myb25=(data_67$yob==1925)*(data_67$sex==1); myb26=(data_67$yob==1926)*(data_67$sex==1); myb27=(data_67$yob==1927)*(data_67$sex==1); myb28=(data_67$yob==1928)*(data_67$sex==1); myb29=(data_67$yob==1929)*(data_67$sex==1); myb30=(data_67$yob==1930)*(data_67$sex==1); myb31=(data_67$yob==1931)*(data_67$sex==1); myb32=(data_67$yob==1932)*(data_67$sex==1); myb33=(data_67$yob==1933)*(data_67$sex==1); myb34=(data_67$yob==1934)*(data_67$sex==1); myb35=(data_67$yob==1935)*(data_67$sex==1); myb36=(data_67$yob==1936)*(data_67$sex==1); myb37=(data_67$yob==1937)*(data_67$sex==1); myb38=(data_67$yob==1938)*(data_67$sex==1); myb39=(data_67$yob==1939)*(data_67$sex==1); #Binaervariabalen den Daten hinzufuegen data_67 = data_67 %>% mutate(myb20, myb21,myb22,myb23,myb24,myb25,myb26,myb27,myb28,myb29,myb30,myb31,myb32,myb33,myb34,myb35,myb36,myb37,myb38,myb39) #alle notwenigen Variablen fuer Regression erhalten writeClipboard(paste0('"',colnames(data_67),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_28", "age_29", "age_30", "age_31", "age_32", "age_33", "age_34", "age_35", "age_36", "age_37", "age_38", "age_39", "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48", "yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "sex_1", "sex_2") dummies5.plus = paste0(dummies5, collapse = '+') #Instrumente fuer die Regression erstellen und mit + verbinden iv = c("myb20", "myb21", "myb22", "myb23", "myb24", "myb25", "myb26", "myb27", "myb28", "myb29", "myb30", "myb31", "myb32", "myb33", "myb34", "myb35", "myb36", "myb37", "myb38", "myb39", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_28", "age_29", "age_30", "age_31", "age_32", "age_33", "age_34", "age_35", "age_36", "age_37", "age_38", "age_39", "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48", "yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "sex_1", "sex_2") iv_plus = paste0(iv, collapse = "+") # Instrumentalvariablenschaetzung durchfuehren, wobei erstmal die zugehörige ols Formel erstellt wird und diese anschliessend mit den Instrumenten zu der Instrumentalavariablenschaetzung zusammengefuegt wird ols_formula = paste0("smoke100 ~ ", dummies5.plus) iv_reg_formula = paste0(ols_formula, "|", iv_plus) #Instrumentalvaribalenschaetzung ausfuehren und benoetigten Parameter abspeichern iv_reg = ivreg(iv_reg_formula, data = data_67) vet3 = coef(iv_reg)[2] #Fehler berechnen sd_5_3 = coeftest(iv_reg, vcov = vcovHC(iv_reg, type = "HC1"))[2,2] #1990: Hier wird genau die selbe Vorgensweise verwendet, wie bei den Daten für 1967/68 data_90 = data_90 %>% dummy_cols(select_columns = 'yob') %>% dummy_cols(select_columns = "sex") #Binaere Variabeln fuer Instrumentalvariablenschaetzung erstellen myb20=(data_90$yob==1920)*(data_90$sex==1); myb21=(data_90$yob==1921)*(data_90$sex==1); myb22=(data_90$yob==1922)*(data_90$sex==1); myb23=(data_90$yob==1923)*(data_90$sex==1); myb24=(data_90$yob==1924)*(data_90$sex==1); myb25=(data_90$yob==1925)*(data_90$sex==1); myb26=(data_90$yob==1926)*(data_90$sex==1); myb27=(data_90$yob==1927)*(data_90$sex==1); myb28=(data_90$yob==1928)*(data_90$sex==1); myb29=(data_90$yob==1929)*(data_90$sex==1); myb30=(data_90$yob==1930)*(data_90$sex==1); myb31=(data_90$yob==1931)*(data_90$sex==1); myb32=(data_90$yob==1932)*(data_90$sex==1); myb33=(data_90$yob==1933)*(data_90$sex==1); myb34=(data_90$yob==1934)*(data_90$sex==1); myb35=(data_90$yob==1935)*(data_90$sex==1); myb36=(data_90$yob==1936)*(data_90$sex==1); myb37=(data_90$yob==1937)*(data_90$sex==1); myb38=(data_90$yob==1938)*(data_90$sex==1); myb39=(data_90$yob==1939)*(data_90$sex==1); data_90 = data_90 %>% mutate(myb20, myb21,myb22,myb23,myb24,myb25,myb26,myb27,myb28,myb29,myb30,myb31,myb32,myb33,myb34,myb35,myb36,myb37,myb38,myb39) #Regressionsvariablen und Instrumente definieren writeClipboard(paste0('"',colnames(data_90),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75", "yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "sex_1", "sex_2") dummies5.plus = paste0(dummies5, collapse = '+') iv = c("myb20", "myb21", "myb22", "myb23", "myb24", "myb25", "myb26", "myb27", "myb28", "myb29", "myb30", "myb31", "myb32", "myb33", "myb34", "myb35", "myb36", "myb37", "myb38", "myb39", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75", "yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "sex_1", "sex_2") iv_plus = paste0(iv, collapse = "+") ols_formula = paste0("smoke100 ~ ", dummies5.plus) iv_reg_formula = paste0(ols_formula, "|", iv_plus) #Regression und Fehlerrechnung durchfuehren iv_reg = ivreg(iv_reg_formula, data = data_90) vet4 = coef(iv_reg)[2] sd_5_4 = coeftest(iv_reg, vcov = vcovHC(iv_reg, type = "HC1"))[2,2] #TLSL: men only: hier wird exact das selbe Vorgehen wie bei dem TSLS Schätzer für men and women verwendet, nur werden hier nur die Daten mit maennlichen Personen verwendet #1967/1968 data_67_men = data_67_men %>% dummy_cols(select_columns = 'yob') #Regressionsvariablen und Instrumente definieren writeClipboard(paste0('"',colnames(data_67_men),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_28", "age_29", "age_30", "age_31", "age_32", "age_33", "age_34", "age_35", "age_36", "age_37", "age_38", "age_39", "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48") dummies5.plus = paste0(dummies5, collapse = '+') iv = c( "yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_28", "age_29", "age_30", "age_31", "age_32", "age_33", "age_34", "age_35", "age_36", "age_37", "age_38", "age_39", "age_40", "age_41", "age_42", "age_43", "age_44", "age_45", "age_46", "age_47", "age_48") iv_plus = paste0(iv, collapse = "+") ols_formula = paste0("smoke100 ~ ", dummies5.plus) iv_reg_formula = paste0(ols_formula, "|",iv_plus) #Iv Schaetzung laufen lassen und Fehler berchnen iv_reg = ivreg(iv_reg_formula, data = data_67_men) vet5 = coef(iv_reg)[2] sd_5_5 = coeftest(iv_reg, vcov = vcovHC(iv_reg, type = "HC1"))[2,2] #1990: Das Vorgehen wir nochmals fuer die Daten von 1990 wiederholt data_90_men = data_90_men %>% dummy_cols(select_columns = 'yob') #Regressionsvariablen und Instrumente definieren writeClipboard(paste0('"',colnames(data_90_men),'"', collapse=", ")) dummies5 = c("vet", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75") iv = c("yob_1920", "yob_1921", "yob_1922", "yob_1923", "yob_1924", "yob_1925", "yob_1926", "yob_1927", "yob_1928", "yob_1929", "yob_1930", "yob_1931", "yob_1932", "yob_1933", "yob_1934", "yob_1935", "yob_1936", "yob_1937", "yob_1938", "yob_1939", "black", "married", "edcat2", "edcat3", "edcat4", "division_1", "division_2", "division_3", "division_4", "division_5", "division_6", "division_7", "division_8", "division_9", "age_53", "age_54", "age_55", "age_56", "age_57", "age_58", "age_59", "age_60", "age_61", "age_62", "age_63", "age_64", "age_65", "age_66", "age_67", "age_68", "age_69", "age_70", "age_71", "age_72", "age_73", "age_74", "age_75") iv_plus = paste0(iv, collapse = "+") dummies5.plus = paste0(dummies5, collapse = '+') ols_formula = paste0("smoke100 ~ ", dummies5.plus) iv_reg_formula = paste0(ols_formula, "|",iv_plus) #IV Schaetzung durchfuehren und Fehler ermitteln iv_reg = ivreg(iv_reg_formula, data = data_90_men) vet6 = coef(iv_reg)[2] sd_5_6 = coeftest(iv_reg, vcov = vcovHC(iv_reg, type = "HC1"))[2,2] #Tabelle 5 erstellen: Nun werden alle erhaltenden Ergbenisse aufbereitet und in Vektoren zusammengefasst, um diese anschliessend in einer Tabelle zu speichern. Die Aufbereitung der Daten wird hier beispielhaft an der ersten Regression erlaeutert dataset = c("1967/68", "1990") #Regressionsergebnisse und Fehler werden auf drei Nachkommastellen gerundet OLS_men_only = c(round(ols_5_1,3), round(ols_5_2,3)) OLS_men_only_sd = c(round(sd_5_1,3),round(sd_5_2,3)) # Nun sollen die Fehler in einem Array gespeichert werden, wobei jeder Fehler mit 2 Klammern umschlossen wird also (), dies ist für die spätere Ausagbe wichtig OLS_men_only_sd_klammer = c() for (i in 0:(length(OLS_men_only_sd ))){ OLS_men_only_sd_klammer[i] = paste0("(",OLS_men_only_sd [i],")") } #TSLS: men and women TSLS_men_women = c(round(vet3,3),round(vet4,3)) TSLS_men_women_sd = c(round(sd_5_3,2),round(sd_5_4,2)) TSLS_men_women_sd_klammer = c() for (i in 0:(length(OLS_men_only_sd ))){ TSLS_men_women_sd_klammer[i] = paste0("(",TSLS_men_women_sd [i],")") } #TSLS: men TSLS_men = c(round(vet5,3),round(vet6,3)) TSLS_men_sd = c(round(sd_5_5,2),round(sd_5_6,2)) TSLS_men_sd_klammer = c() for (i in 0:(length(OLS_men_only_sd ))){ TSLS_men_sd_klammer[i] = paste0("(",TSLS_men_sd [i],")") } #Tabelle 5 mit Schaetzungsergebnissen wird erstellt tab5 = data.frame(dataset, OLS_men_only, OLS_men_only_sd_klammer, TSLS_men_women, TSLS_men_women_sd_klammer, TSLS_men, TSLS_men_sd_klammer) # Nun werden die Ergebnisse mit Hilfe des kable packages schoen aufbereitet angezeigt kable(tab5, booktabs = TRUE,escape = FALSE, col.names = c(" ","( 1 )"," ","( 2 )"," ","( 3 )"," "), align = c("lrlrlrl"), caption = "Table 5-Impact of Veteran Status on Smoking Behavior") %>% add_header_above(c(" " = 1, "OLS: men only" = 2, "TSLS:men and women" = 2, "TSLS:men only" = 2)) %>% kable_styling() %>% footnote(general = "Eicker-White standard errors are reported in parentheses.", general_title = "Notes: ", footnote_as_chunk = T) ``` In Tabelle 5 sind die Ergebnisse der vorher beschriebenen Schätzungen zu sehen. Besonders auffällig ist dabei der große Unterschied der Schätzungen des OLS Modells und der Instrumentalvariablenschätzungen (TSLS:= Two stage least square). Dieser Unterschied ist auf das vorher schon angesprochene Endogenitätsproblem zurückzuführen, da vermutlich eher weniger Personen, die Rauchen, das Auswahlverfahren des Militärs bestehen, und dadurch der eigentliche Effekt des Veteranenstatus eher systematisch unterschätzt wird. Diese Vermutung lässt sich durch die Betrachtung der Instrumentalvariablenschätzungen bestätigen. Die Unterschiede, welche bei den Varianten des Instrumentalvariablenschätzers auftreten, können vernachlässigt werden. Insgesamt lässt sich nach dem Instrumentalvariablenschätzer auf eine Erhöhung der Wahrscheinlichkeit zu Rauchen durch den Militärservice von 23.7-34.6 Prozentpunkte schließen. # Tabelle 6 Die Tabelle 6 dient dazu die Ergebnisse der vorherigen Analysen zusammenzufassen, um ein endgültiges Fazit zu ziehen. Da die Frage beantwortet werden soll, welchen Einfluss das Rauchen im Militär auf die Sterblichkeit der Veteranen hat, werden hier nur die Todesursachen "Ischemic heart disease" und "Lung cancer" betrachtet, da diese besonders oft durch Rauchen hervorgerufen werden. Die erste Spalte enthält dabei die Sterblichkeitsraten insgesamt für diese beiden Todesursachen. Diese wurden aus der ersten Tabelle bzw. Darstellung entnommen. Die zweite Spalte enthält, welchen Anteil der Todesfälle durch die Todesursache auf Rauchen zurückgeführt werden kann. Diese Werte wurden aus dem Paper entnommen. In Spalte 3 werden nun die Todesraten der Todesursachen pro 1000 Raucher durch das Rauchen berechnet. Dazu wird einfach der Wert in Spalte 1 mit dem Wert in Spalte 2 multipliziert und dies durch 750 geteilt, da es 750 Raucher/1000 Männer gibt. Der Wert in Spalte 3 wird nun mit dem Wert des TSLS-Schätzers (men and women, 1967/68) der Tabelle 5 multipliziert, wodurch die zusätzlichen Todesfälle/1000 durch das Veteranenrauchen ermittelt werden. In Spalte 5 und 6 werden nun noch diese Werte durch die Schätzer in Tabelle 2 und 3 (Spalte 2) geteilt, um den Anteil der Veterantoten durch das Rauchen im Militär herauszufinden. ```{r} #Tabelle 6 erstellen: #Betrachtete Todesursachen definieren causes = c("Ischemic heart disease","Lung cancer") #Vektor mit Sterblichkeitsrate der definierten Todesursachen erstellen (aus Tabelle 1) death_rate_1000 = c(ischem_m, lungcan_m) #Vektor mit Prozentangaben, wie viel der Todesfaelle auf Rauchen zurueckzufuehren sind percent_death_smoking = c(40, 87) # Die Todesraten pro 1000 Raucher durchs Rauchen fuer die beiden Krankheiten death_rate_1000_smoking = c(ischem_m*0.4/750*1000, lungcan_m*0.87/750*1000) #Zusaetzliche Tote pro 1000 wegen dem Rauchen der Veteranen add_death_1000 = c(death_rate_1000_smoking[1]*tab5[1,4], death_rate_1000_smoking[2]*tab5[1,4]) #prozentuale Anteil der Veterantoten durch Rauchen. Einmal auf die Ergebnisse der Tabelle 2 bezogen und einmal auf die Ergebnisse der Tabelle 3 percent_vet_effect_smoking_2 = c(add_death_1000[1]/as.numeric(tab2[2,4])*100, add_death_1000[2]/as.numeric(tab2[3,4])*100) percent_vet_effect_smoking_3 = c(add_death_1000[1]*100/as.numeric(tab3[2,4]), add_death_1000[2]*100/as.numeric(tab3[3,4])) # Nun werden die Ergebisse in einer Tabelle zusammengefassst tab6 = data.frame(causes, round(death_rate_1000,2), percent_death_smoking, round(death_rate_1000_smoking,2), round(add_death_1000,3), round(percent_vet_effect_smoking_2,0), round(percent_vet_effect_smoking_3,0)) # "tab6" wird nun schoen mit dem kable-package angezeigt kable(tab6, booktabs = TRUE, escape = FALSE, col.names = c(" ","( 1 )","( 2 )","( 3 )","( 4 )","( 5 )","( 6 )"), align = c("ccccccc"), caption = "Table 6 - Percentage of Veteran Mortality Effect Explained by Military Induced Smoking") %>% add_header_above(c(" " =1,"Deaths rate per 1000 (Table 1)" = 1,"% Deaths caused by smoking " = 1,"Death rate per 1000 smokers due to smoking" = 1,"Additional deaths per 1000 due to veteran smoking" = 1,"% Veteran effect due to smoking (Table 2)" = 1,"% Veteran effect due to smoking (Table 3)"= 1)) %>% kable_styling() ``` Besonders die Ergebnisse in Spalte 5 und 6 tragen zu einer Einschätzung bei, welche Folgen das Rauchen im Militär hat. So lässt sich durch Anwendung der hier verwendeten Modelle vermuten, dass 64-79% der zusätzliche Todesfälle durch "Ischemic heat disease" bei Veteranen (verglichen mit nicht Veteranen) auf das Rauchen im Militär zurückzuführen sind. Bei der Todesursache "Lung cancer" sind es 35-58%.