CSSR-Presentation.Rmd

---
title: "Drink and Thrive? <br> German Economics of Binge-Drinking"
subtitle: "CSSR Presentation"
author: "Torben and Alex"
date: "02 December 2016"
output: ioslides_presentation
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)

# Clear Global environment
rm(list=ls())

# Setting Working directory
try(setwd("/home/torben/GIT/Pair_Assignment_2"), silent = TRUE)
try(setwd("D:/Eigene Datein/Dokumente/Uni/Hertie/Materials/Collaborative Social Science Data Analysis/CSSR_Project"), silent = TRUE)

source("main.R")

# Collect packages/libraries we need:
packages <- c("fontcm", "psych", "plotly", "Zelig")
# can we add a short list so we have an overview for what we need the packages?
# psych allows for summary statistics
# plotly makes html graphs

# install packages if not installed before
for (p in packages) {
  if (p %in% installed.packages()[,1]) {
    require(p, character.only=T)
  }
  else {
    install.packages(p, repos="http://cran.rstudio.com", dependencies = TRUE)
    require(p, character.only=T)
  }
}
rm(p, packages)

# Initiate Label Names
statelist_code <- c("DE-BW", "DE-BY", "DE-BE", "DE-BB", "DE-HB", "DE-HH", 
                    "DE-HE", "DE-MV", "DE-NI", "DE-NW", "DE-RP", "DE-SL",
                    "DE-SN", "DE-ST", "DE-SH", "DE-TH")

statelist_name <- c("Baden-Württemberg", "Bayern", "Berlin", "Brandenburg", 
                    "Bremen", "Hamburg", "Hessen", "Mecklenburg-Vorpommern", 
                    "Niedersachsen", "Nordrhein-Westfalen", "Rheinland-Pfalz", 
                    "Saarland", "Sachsen", "Sachsen-Anhalt", 
                    "Schleswig-Holstein", "Thüringen")
statelist.code.short <- substr(statelist_code, 4, 5)
statelist_cities <- c("DE-BE", "DE-HB", "DE-HH")
statelist_other <- c("DE-BW", "DE-BY", "DE-BB", "DE-HE", "DE-MV", "DE-NI", "DE-NW", "DE-RP", "DE-SL",
                    "DE-SN", "DE-ST", "DE-SH", "DE-TH")
statelist_other_wo_treatment <- c("DE-BY", "DE-BB", "DE-HE", "DE-MV", "DE-NI", "DE-NW", "DE-RP", "DE-SL",
                    "DE-SN", "DE-ST", "DE-SH", "DE-TH")
```


## Structure

- Introduction
- Data 
- Descriptive Statistics
- Empirical Strategy
- Results


## Trend Sport of the 21st Century?

<img src="presentation_pictures/drunk.jpeg" width="90%" height="90%"/>


## Background

**Binge Drinking -- A Problem?**

German crime statistics: violent crime increasingly committed under the influence of alcohol

Adolescents especially affected

**Research Questions**

Similar to *Marcus and Siedler, 2015*:

What explains regional differences? Correlation with economic situation?

Night sales ban on alcohol: any effect on binge drinking in Baden-Wuerttemberg?

## Data

------------------------------------------------------------------------
Variables                       Reasoning
------------------------------  ----------------------------------------
F10.0 cases                     Number of people hospitalized with 
per 1000                        acute alcohol intoxication

K70 cases                       Number of people hospitalized diagnosed
per 1000                        with alcohol-related liver disease

GDP per capita                  Value and change between years as 
                                indicator of economic development

(Youth) unemployment rate       Indicator of economic situation for 
                                age groups

Population density              Indicator for urbanization

Treatment dummies               Capture the policy measure effect
(local & temporal)

------------------------------------------------------------------------

## Data Visuals - F10.0 per Age Group

```{r echo = FALSE, fig.width=8, fig.height=5} 
# Bar diagram: F10.0 cases per 1000 per age group in DE-DE
PL1 <- TOTAL %>% filter(GENDER == "all", AGE!="all", STATE=="DE-DE")
ggplot(data=PL1, aes(x=AGE, y=F100_p1000)) +
  geom_bar(stat="identity", fill = "#1E90FF") +
  xlab("Age Group") +
  ylab("F10.0 cases per 1000 people") +
  ggtitle("") +
  #theme_bw() +           # do we want it to be bw or grey?
  guides(fill=FALSE) +
  theme(text  = element_text(colour = "#797979"),
        axis.text.x  = element_text(angle=90, vjust=0.5),
        axis.title.y=element_text(margin=margin(10,10,0,10)))
```

## Data Visuals - K70 per Age Group

```{r echo = FALSE, fig.width=8, fig.height=5} 
# Bar diagram: K70 cases per 1000 per age group
ggplot(data=PL1, aes(x = AGE, y = K70_p1000, fill = AGE)) +
  geom_bar(stat="identity", fill = "#1E90FF") +
  xlab("Age Group") +
  ylab("K70 cases per 1000 people") +
  ggtitle("") +
  guides(fill=FALSE) +
  theme(text  = element_text(colour = "#797979"),
        axis.text.x  = element_text(angle=90, vjust=0.5),
        axis.title.y=element_text(margin=margin(10,10,0,10)))
```


## Data Visuals - F10.0 Cases in Germany

```{r echo = FALSE, fig.width=4, fig.height=5, out.extra='style="float:left"'} 
spplot(map.Germany.2000.F100_p1000, 
       zcol = "value", 
       col.regions = colors,
       main=list(label="2000: F10.0 diagnoses\nper 1000 pp all age groups"),
       at = seq(0, 2.2, by = 0.1))
```
```{r echo = FALSE, fig.width=4, fig.height=5, out.extra='style="float:right"'} 
spplot(map.Germany.2014.F100_p1000, 
       zcol = "value", 
       col.regions = colors,
       main=list(label="2014: F10.0 diagnoses\nper 1000 pp all age groups"),
       at = seq(0, 2.2, by = 0.1),
       colorkey=FALSE)
```


## Data Visuals - F10.0 Increase Over Time

```{r echo = FALSE} 
#  data frame for plot
DS0 <- TOTAL %>% filter(GENDER == "all", AGE=="all") %>% 
  select(STATE, YEAR, F100_p1000, F102_p1000, K70_p1000)
DS0$GROUP <- "Non-City States"
DS0$GROUP[DS0$STATE == "DE-BW"] <- "Treatment"
DS0$GROUP[DS0$STATE == "DE-BE" | DS0$STATE == "DE-HB" | DS0$STATE == "DE-HH"] <- "City"

# Start with Treatment
F100Plot <- plot_ly(data = filter(DS0, STATE == "DE-BW"),
                    x = ~YEAR,
                    y = ~F100_p1000,
                    type = "scatter",
                    mode = "lines",
                    name = "Treatment State BW",
                    #color = ~GROUP[1],
                    legendgroup = ~GROUP[1],
                    line = list(color = "red", width = 4)) %>%
  layout(title = "", xaxis = list(title = "Year", titlefont = list(size = 18)),
         yaxis = list(title = "F10.0 cases per 1000 people", titlefont = list(size = 18)),
         legend = list(font = list(size = 14)))

# Non-City States
for (i in 1:12){
  F100Plot <- add_trace(F100Plot, 
                        data = filter(DS0, STATE == statelist_other_wo_treatment[i]), 
                        text = ~STATE[1], 
                        legendgroup = ~GROUP[1], 
                        #color = ~GROUP[1], 
                        name = paste("State", substr(statelist_other_wo_treatment[i], 4, 5)),
                        line = list(color = colorRampPalette(brewer.pal(9, "Greys"))(16)[i+4], width = 2))
}
# City States
for (i in 1:3){
  F100Plot <- add_trace(F100Plot, 
                        data = filter(DS0, STATE == statelist_cities[i]), 
                        #text = ~STATE[1], 
                        #legendgroup = ~GROUP[1], 
                        #color = colorRampPalette(brewer.pal(5, "Reds"))(5)[i+2], #~STATE[1], 
                        line = list(color = colorRampPalette(brewer.pal(5, "Blues"))(5)[i+2], width = 2),
                        name = paste("City State", substr(statelist_cities[i], 4, 5)))
}
F100Plot
```

## Data Visuals - F10.0 Increase - 15-19 y/o

```{r echo = FALSE} 
#  data frame for plot
DS0 <- TOTAL %>% filter(GENDER == "all", AGE=="15-19y") %>% 
  select(STATE, YEAR, F100_p1000, F102_p1000, K70_p1000)
DS0$GROUP <- "Non-City States"
DS0$GROUP[DS0$STATE == "DE-BW"] <- "Treatment"
DS0$GROUP[DS0$STATE == "DE-BE" | DS0$STATE == "DE-HB" | DS0$STATE == "DE-HH"] <- "City"

# Start with Treatment
F100Plot <- plot_ly(data = filter(DS0, STATE == "DE-BW"),
                    x = ~YEAR,
                    y = ~F100_p1000,
                    type = "scatter",
                    mode = "lines",
                    name = "Treatment State BW",
                    #color = ~GROUP[1],
                    legendgroup = ~GROUP[1],
                    line = list(color = "red", width = 4)) %>%
  layout(title = "", xaxis = list(title = "Year", titlefont = list(size = 18)),
         yaxis = list(title = "F10.0 cases per 1000 people", titlefont = list(size = 18)),
         legend = list(font = list(size = 14)))

# Non-City States
for (i in 1:12){
  F100Plot <- add_trace(F100Plot, 
                        data = filter(DS0, STATE == statelist_other_wo_treatment[i]), 
                        text = ~STATE[1], 
                        legendgroup = ~GROUP[1], 
                        #color = ~GROUP[1], 
                        name = paste("State", substr(statelist_other_wo_treatment[i], 4, 5)),
                        line = list(color = colorRampPalette(brewer.pal(9, "Greys"))(16)[i+4], width = 2))
}
# City States
for (i in 1:3){
  F100Plot <- add_trace(F100Plot, 
                        data = filter(DS0, STATE == statelist_cities[i]), 
                        #text = ~STATE[1], 
                        #legendgroup = ~GROUP[1], 
                        #color = colorRampPalette(brewer.pal(5, "Reds"))(5)[i+2], #~STATE[1], 
                        line = list(color = colorRampPalette(brewer.pal(5, "Blues"))(5)[i+2], width = 2),
                        name = paste("City State", substr(statelist_cities[i], 4, 5)))
}
F100Plot
```


## Empiricial Strategy

- Dependent variable: Alcohol-related hospitalization (per 1000)

**Question 1: Differences between states**

Method: Multiple regression

- *Variation*: Levels and differences
- *Robustness check:* Different years and long vs. short-term alcohol diagnosis

**Question 2: Effect of night-sale-ban**

Method: Difference-in-Difference

- *Variation*: Control variables and number of periods
- *Robustness check:* different control group compositions

## Results Question 1
```{r Result1, echo = FALSE}
# Simulation of model 1 
Z1.2014 <- zelig(F100_p1000 ~ GDP_P_C + UR.LF + PD, cite = FALSE,
            data = M1.2014, model = "normal")

# Simulation of relationship F10.0 diagnoses and GDP per capita
setZ.GDP.PP <- setx(Z1.2014, GDP_P_C = 15:46)
simZ.GDP.PP <- sim(Z1.2014, x = setZ.GDP.PP)

# Simulation of relationship F10.0 diagnoses and unemployment rate
setZ.UR.LF <- setx(Z1.2014, UR.LF = 5:21)
simZ.UR.LF <- sim(Z1.2014, x = setZ.UR.LF)

# Simulation of model 2 
Z2.F100 <- zelig(F100.D ~ GDP.D + UR.LF.D -1, cite = FALSE,
                 data = ungroup(M2), model = "normal")
  
# Simulation of relationhip F10.0 diagnoses and GDP per capita
setZ2.GDP.D <- setx(Z2.F100, GDP.D = -4:4)
simZ2.GDP.D <- sim(Z2.F100, x = setZ2.GDP.D)


# Simulation of relationhip F10.0 diagnoses and GDP per capita
setZ2.UR.LF.D <- setx(Z2.F100, UR.LF.D = -3:4)
simZ2.UR.LF.D <- sim(Z2.F100, x = setZ2.UR.LF.D)
```


```{r R1GraphTop, echo = FALSE, fig.width=4, fig.height=2.8, out.extra='style="float:left"'}
ci.plot(simZ.GDP.PP, ylab ="Expected conditional F10.0", xlab = "GDP level",
        main = "Model 1: F10.0 Diagnoses 2014")
# title(xlab="GDP level", line=-1, cex.lab=1)
ci.plot(simZ.UR.LF, ylab ="Expected conditional F10.0", xlab = "Unemployment level",
        main = "Model 1: F10.0 Diagnoses 2014")

```

```{r R1GraphBelow, echo = FALSE, fig.width=4, fig.height=2.8, out.extra='style="float:left"'}
ci.plot(simZ2.GDP.D, ylab ="Expected conditional F10.0", xlab = "GDP change",
        main = "Model 2: Changes F10.0 Diagnoses")
ci.plot(simZ2.UR.LF.D, ylab ="Expected conditional F10.0", xlab = "Unemployment change",
        main = "Model 2: Changes F10.0 Diagnoses")
```

## Result Night Sale Ban
F10.0 Diagnoses for 15-19 year olds (control: 13 non-city states)
```{r Result Night Ban, results = "asis", echo = FALSE} 
stargazer::stargazer(mod5.age15.13,
                     covariate.labels = c("BW night sale ban",
                                          "Youth unemployment rate", "GDP per capita"),
                     model.numbers  = FALSE,
                     dep.var.labels = c("F10.0 cases per 1000 PP"),                           
                     title = 'Model 5: multi-period panel differences-in-differences',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)
```



## Conclusion

Marcus and Sielder (2014) find about 7% decrease, we find about 20% decrease

```{r echo = FALSE} 
plot_ly(bancomparison, x = ~YEAR, y = ~WBAN, type = 'bar', name = 'With Ban') %>%
  add_trace(y = ~WOBAN, name = 'Without Ban') %>%
  layout(xaxis = list(title = 'Year', titlefont = list(size = 18)), 
         yaxis = list(title = 'F10.0 Cases for 15 to 19 year-olds in DE-BW',
                      titlefont = list(size = 18)),
         legend = list(x = 0.8, y = 1, font = list(size = 18)),
         barmode = 'group')
```

## Conclusion 2

Problems of our research design:

- Availability of data 
    - without paying only annual and state-level data available
    - 14 years x 16 states = 224 obs. 
    - access to monthly and county level data would allow
    - 168 months x 402 counties = 67536 obs.
- limited amount of control variables
    - few meaningful state level variables available
    - but: less of a problem in Diff-in-Diff method


## Appendix 1

```{r Model1, results = "asis", echo = FALSE}
mod1.labels <- c("GDP per capita", "Unemployment rate", "Beer tax", "Population density", "(Intercept)")
stargazer::stargazer(mod1.2000, mod1.2007, mod1.2014,
                     covariate.labels = mod1.labels,
                     column.labels = c("2000", "2007", "2014"),
                     model.numbers  = FALSE,
                     dep.var.labels = "F10.0 Diagnoses per 1000 capita",
                     title = 'Regression results for Model 1 for 2000, 2007 and 2014',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)
```

## Appendix 2

```{r Model2, results = "asis", echo = FALSE} 
stargazer::stargazer(mod2.F100, mod2.F102, mod2.K70,
#                     column.labels = c("F10.0", "F10.2", "K70"),
                     covariate.labels = c("GDP Change", "Unemployment Change", "Beer Tax Change", "(Intercept)"),
                     model.numbers  = FALSE,
                     dep.var.labels = c("Change in F10.0","Change in F10.2","Change in K70"),
                     title = 'Regression results for Model 2 with first differenced data',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)
```

## Appendix 3

```{r Model3, results = "asis", echo = FALSE} 
stargazer::stargazer(mod3.16, mod3.13, 
                     column.labels = c("All States", "All but City States"),
                     covariate.labels = c("BW night sale ban", "Post-treatment dummy", "Interaction", "(Intercept)"),
                     model.numbers  = FALSE,
                     dep.var.labels = c("F10.0 cases per 1000 PP"),                    
                     title = 'Model 3: simple differences-in-differences',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)

```

## Appendix 4

```{r Model4, results = "asis", echo = FALSE} 
stargazer::stargazer(mod4, 
#                     column.labels = c("All States"),
                     covariate.labels = c("BW night sale ban", "Post-sales ban dummy", 
                                          "GDP per capita", "Youth unemployment rate", 
                                          "Interaction dummy", "(Intercept)"),
                     model.numbers  = FALSE,
                     dep.var.labels = c("F10.0 cases per 1000"),                           
                     title = 'Model 4: simple diff-in-diff with controls',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)
```

## Appendix 5

```{r Model5, results = "asis", echo = FALSE} 

stargazer::stargazer(mod5, mod5.age15.16, mod5.age15.13,
                     column.labels = c("All States", "15-19y", "15-19y + control: non-city"),
                     covariate.labels = c("BW night sale ban", "Youth unemployment rate", "GDP per capita"),
                     model.numbers  = FALSE,
                     dep.var.labels = c("F10.0 cases per 1000 PP"),                           
                     title = 'Model 5: multi-period panel differences-in-differences',
                     omit.stat = c("rsq", "f"),
                     digits = 2, type = 'html', header = FALSE)
```

## Appendix 6

```{r Maps2, message=FALSE, echo = FALSE} 

#library(googleVis)

#require(datasets)

#states <- data.frame(state.name, state.x77)

#GeoStates <- gvisGeoChart(states, "state.name", "Illiteracy",
#                            options=list(region="DE", 
#                            displayMode="regions", 
#                            resolution="provinces",
#                            width=600, height=400))

#plot(GeoStates)
```