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
#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)
#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.

# 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).

#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)
## [1] 0.9638621
#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.

#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)
## [1] 0.8491059
#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.

#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 (<Fehler>), 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 (<Fehler>), 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(" ","&#40 1 &#41"," ","&#40 2 &#41"," ","&#40 3 &#41","&#40 4 &#41","&#40 5 &#41"), 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)
Table 2 - Impact of Veteran Status on Mortality, Male-Only Sample

Implied mortality ratesa
Veteran effect
Mean rate
Nonvets
Vets
( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 )
All mortality causes 3.562 (0.97) 3.424 (0.97) 15.406 13.2 16.6
Ischemic hert disease 0.965 (0.27) 0.74 (0.21) 3.979 3.5 4.2
Lung cancer 0.986 (0.18) 1.027 (0.2) 1.851 1.2 2.2
Colon cancer 0.024 (0.04) 0.016 (0.04) 0.388 0.4 0.4
Cerebrovascular diseases 0.043 (0.05) -0.015 (0.05) 0.642 0.7 0.6
Respiratory diseases 0.255 (0.09) 0.309 (0.11) 0.91 0.7 1
Accidents and suicides 0.233 (0.08) 0.097 (0.05) 0.645 0.6 0.7
Observations 579 579
Unrestricted age dummies Yes Yes
Linear cohort trend Yes Yes
Unrestricted year dummies No Yes
Notes: 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.
a Implied mortality rates based on the estimates in column 2

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.

#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 (<Fehler>), 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 (<Fehler>), 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(" ","&#40 1 &#41"," ","&#40 2 &#41"," ","&#40 3 &#41","&#40 4 &#41","&#40 5 &#41"), 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)
Table 3 - Impact of Veteran Status on Mortality, male and Female Sample

Implied mortality ratesa
Veteran effect
Mean rate
Nonvets
Vets
( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 )
All mortality causes 3.847 (0.44) 3.754 (0.43) 15.406 12.9 16.7
Ischemic hert disease 1.147 (0.1) 0.919 (0.09) 3.979 3.4 4.3
Lung cancer 1.643 (0.22) 1.688 (0.23) 1.851 0.7 2.4
Colon cancer -0.094 (0.02) -0.086 (0.02) 0.388 0.4 0.4
Cerebrovascular diseases -0.054 (0.03) -0.041 (0.03) 0.642 0.7 0.6
Respiratory diseases 0.809 (0.12) 0.836 (0.13) 0.91 0.4 1.2
Accidents and suicides -0.092 (0.02) -0.114 (0.02) 0.645 0.7 0.6
Observations 1158 1158
Unrestricted sex*age dummies Yes Yes
Unrestricted cohort dummies Yes Yes
Unrestricted year dummies No Yes
Notes: 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.
a Implied mortality rates based on the estimates in column 2

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.

#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 (<Fehler>), 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(" ","&#40 1 &#41"," ","&#40 2 &#41"," ","&#40 3 &#41"," "), 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)
Table 5-Impact of Veteran Status on Smoking Behavior
OLS: men only
TSLS:men and women
TSLS:men only
( 1 ) ( 2 ) ( 3 )
1967/68 0.078 (0.005) 0.276 (0.03) 0.237 (0.14)
1990 0.126 (0.004) 0.346 (0.02) 0.301 (0.03)
Notes: Eicker-White standard errors are reported in parentheses.

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.

#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(" ","&#40 1 &#41","&#40 2 &#41","&#40 3 &#41","&#40 4 &#41","&#40 5 &#41","&#40 6 &#41"), 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()
Table 6 - Percentage of Veteran Mortality Effect Explained by Military Induced Smoking
Deaths rate per 1000 (Table 1)
% Deaths caused by smoking
Death rate per 1000 smokers due to smoking
Additional deaths per 1000 due to veteran smoking
% Veteran effect due to smoking (Table 2)
% Veteran effect due to smoking (Table 3)
( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 )
Ischemic heart disease 3.98 40 2.12 0.586 79 64
Lung cancer 1.85 87 2.15 0.593 58 35

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%.