Let's check how many units are there in every district and which districts have treated units

We can use the base R table() function, which performs categorical tabulation of data with the variable and its frequency, to check how many people live on each district.

table(soda_tax_df$district) %>%
  knitr::kable(col.names = c("District", "Inhabitants")) %>% # create kable table
  kableExtra::kable_styling() # view kable table

Even more, we can create a two-way cross-table to check how many treated units live in each district.

table(soda_tax_df$treatment, soda_tax_df$district) %>% 
  knitr::kable() %>% # create kable table
  kableExtra::kable_styling() # view kable table

Let's check the distribution of sugar-added drinks consumption at different points in time

We can also utilize summary() to check the distribution of our variables. In this case, we can look at the distribution of our variable of interest at the two points in time:

summary(soda_tax_df$pre_tax)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   400.6   506.6   509.7   610.8  1687.6

ggplot(soda_tax_df, aes(x = pre_tax)) +
  geom_vline(xintercept = 509.7, linetype = "dashed") +
  geom_density(fill = "#f1a340", alpha = 0.3) +
  theme_minimal() +
  labs(x = "Pre-tax soda consumtion (oz)",
       y = "Density")+
  ggtitle("Soda consumption before the tax, in all districts")

summary(soda_tax_df$post_tax)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0   327.0   445.9   449.8   566.1  1705.6

ggplot(soda_tax_df, aes(x = post_tax)) +
  geom_density(fill = "#998ec3", alpha = 0.3) +
  geom_vline(xintercept = 449.8, linetype = "dashed") +
  theme_minimal() +
  labs(x = "Post-tax soda consumtion (oz)",
       y = "Density") +
  ggtitle("Soda consumption after the tax, in all districts")

a. Calculating without time

after_model <- lm(post_tax ~ treatment, data = soda_tax_df)

stargazer(after_model, type = "text")

We could read this result substantively as: those who lived in districts were the tax was implemented consumed on average 146.9 ounces less of sugar-added drinks per week compared to those who lived in districts were the tax was not put in place. This calculation would give us a comparison of the treatment and control groups after treatment.

To believe the results of our after_model, we would need to believe that the mean ignorability of treatment assignment assumption is fulfilled. We would have to think carefully about possible factors that could differentiate our treatment and control groups. We use a treatment indicator based on the districts where the measure was able to be implemented. Treatment was not fully randomly assigned, so there may be lots of potential confounders that create baseline differences in the scores for those living in Old Eagleton compared to those in Snake Lounge, which also affect the after-treatment comparisons.

If we think about the mechanics behind this naive calculation, we are just comparing the average observed outcomes for those treated and not treated after the tax was imposed:

soda_tax_df %>% 
  dplyr::group_by(treatment) %>% 
  dplyr::summarize(mean = mean(post_tax)) %>% 
  knitr::kable(col.names = c("Treatment", "Average after tax")) %>%
  kableExtra::kable_styling()

b. Including the time dimension (First differences on treatment indicator)

We can introduce the time component to our calculation by incorporating the pre-treatment levels of sugar-added drink consumption, which gives us the diff-in-diff estimand. We could calculate this in a fairly straightforward manner by creating a variable assessing the change:

• change: The difference in sugar-added drink consuption between post- and pre-tax

soda_tax_df <- soda_tax_df %>%
  dplyr::mutate(change = post_tax - pre_tax) #

did_model <- lm(change ~ treatment, data = soda_tax_df)

stargazer(did_model, after_model, type = "text")

We could read this result substantively as: those who lived in districts were the tax was implemented consumed on average 149.7 ounces less of sugar-added drinks per week compared to those who lived in districts were the tax was not put in place. This calculation would give us the change, or difference, in sugar-added drink consumption for treatment and control groups.

To believe the results of our did_model, we would need to believe that there are parallel trends between the two groups.

Note: when simulating the data the post_tax was defined as: \(post_tax = 15 + pre\_tax - 150(treatment) + error\)

c. Including the time dimension (Regression formulation of the DiD model)

\[Y_{it} = β_0 + β_1D^*_i + β_2P_t + β_{DD}D^∗_i × P_t + q_{it} \]

We need to convert our data to a long format to render the time and treatment dummy variables. This is how our data look like:

head(soda_tax_df, 6)

## # A tibble: 6 x 6
##      id district     treatment pre_tax post_tax change
##   <dbl> <chr>            <dbl>   <dbl>    <dbl>  <dbl>
## 1     1 Snake Lounge         0   1688.    1706.  17.9 
## 2     2 Snake Lounge         0    427.     438.  11.0 
## 3     3 Snake Lounge         0    566.     560.  -6.80
## 4     4 Snake Lounge         0    607.     624.  17.0 
## 5     5 Snake Lounge         0    573.     607.  34.2 
## 6     6 Snake Lounge         0    496.     502.   5.72

We will utilize the pivot_longer() function from tidyr to format our data frame.

soda_tax_df_long <- soda_tax_df %>%
  dplyr::select(-change) %>% # select the columns we are interested in
  tidyr::pivot_longer(cols = c(pre_tax, post_tax), names_to = "period", values_to = "soda_drank") %>% # grab columns, put names to new variable period and values to new variable soda_drank
  dplyr::mutate(after_tax = ifelse(period == "post_tax", 1, 0)) # create dummy for period

  head(soda_tax_df_long, 6)

## # A tibble: 6 x 6
##      id district     treatment period   soda_drank after_tax
##   <dbl> <chr>            <dbl> <chr>         <dbl>     <dbl>
## 1     1 Snake Lounge         0 pre_tax       1688.         0
## 2     1 Snake Lounge         0 post_tax      1706.         1
## 3     2 Snake Lounge         0 pre_tax        427.         0
## 4     2 Snake Lounge         0 post_tax       438.         1
## 5     3 Snake Lounge         0 pre_tax        566.         0
## 6     3 Snake Lounge         0 post_tax       560.         1

We can see that under our long format, we have two entries for every individual. We have our variable after_tax, which represents \(P_t\), where 0 and 1 are pre- and post-tax periods respectively. We can now render our regression based on the formula:

did_long <- lm(soda_drank ~ treatment + after_tax + treatment * after_tax, data = soda_tax_df_long) #running our model

did_long_clustered_se <- coeftest(did_long, vcov=vcovHC(did_long,type="HC0",cluster="district")) #clustering out standard errors at the district level

stargazer(did_long_clustered_se, type = "text")

(We can think about the results of this regression in terms of the lab's slide 10: we are mostly interested in the coefficient for the interaction term in our regression equation).

soda_tax_df_long %>% 
  dplyr::group_by(period, treatment) %>% #group by period and treatment
  dplyr::summarize(soda_consumed = mean(soda_drank)) %>% #render averages
  tidyr::pivot_wider(names_from = period, values_from = soda_consumed) %>%
  dplyr::select(treatment, pre_tax, post_tax) %>%
  dplyr::arrange(desc(treatment)) %>%
  knitr::kable(col.names = c("Treatment", "Pre-tax", "Post-tax"),
               digits = 2) %>%
  kableExtra::kable_styling(full_width = F) %>%
  kableExtra::add_header_above(c("", "Period" = 2)) #add header for period

soda_tax_df %>% 
  select(-id, -district) %>% 
  group_by(treatment) %>%
  summarize_all(funs(mean)) %>% 
  pivot_longer(cols = c(pre_tax, post_tax), 
               names_to = "period", 
               values_to = "soda_drank") %>% 
  mutate(time = ifelse(period == "post_tax", 2, 1)) %>% 
  group_by(treatment) %>% 
  ggplot(aes(x = factor(time), y = soda_drank, color = factor(treatment))) +
  geom_point() +
  geom_line(aes(x = time, y = soda_drank)) +
  labs(x = "Time periods", y = "Ounces of soda drank per week", color = "Treatment group")+
  theme_minimal() 

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.