library(AER)
data("Fatalities")

Fatalities |> head(5) |> kable()

#help(Fatalities)

glimpse(Fatalities)

Fatalities |> 
  as_tibble() |> 
  select(-state) |> #for readability
  modelsummary::datasummary_skim()

fatalities <- as_tibble(Fatalities) |>
  remove_missing() |> 
  mutate(
    state = str_to_upper(state),
    year = as.integer(paste(year)),
    jail = (jail == "yes"),
    jail_name = ifelse(jail, "Jail penalty", "No jail penalty")
  ) 

fatal_jail <- fatalities |> 
  ggplot(aes(x = jail_name, y = fatal)) + 
  # geom_hline(aes(yintercept = mean(fatal, na.rm = TRUE))) +
  geom_jitter(width = 0.25) +
  labs(
    title = "Jail penaly for drunk driving and fatal accidents",
    y = "Number of fatal accidents per state",
    x = "State law on drunk driving"
  )

fatal_jail_log <- fatal_jail + 
  scale_y_log10() + 
  labs(y = "Number of fatal accidents\nper state (log scale)")

# Evolution of law
law_evol <- fatalities |> 
  group_by(year) |> 
  # summarise(prop_jail = mean(jail, na.rm = TRUE))
  summarise(prop_jail = mean(jail, na.rm = TRUE)) |> 
  ggplot(aes(year, prop_jail)) +
  geom_line() +
  labs(
    title = "Adoption of jail sentence for drunk driving",
    y = "Proportion of states with a jail law", 
    x = NULL
  ) 

# Making the map
states_sf <- tigris::states(
    cb = TRUE, resolution = "20m", year = 2024, progress_bar = FALSE) |>
  tigris::shift_geometry() |> 
  rename(state = STUSPS)

fatalities_sf <- fatalities |> 
  filter(jail) |> 
  group_by(state) |> 
  mutate(first_year = min(year, na.rm = TRUE)) |> 
  ungroup() |> 
  filter(year == first_year) |> 
  dplyr::full_join(states_sf, by = join_by(state)) |> 
  sf::st_as_sf()

law_map <- fatalities_sf |> 
  ggplot() +
  geom_sf(aes(fill = first_year), color = "white", linewidth = 0.1) +
  scale_mediocre_c(pal = "portal", gradient = "left") +
  theme_mediocre_map() +
  labs(
    title = "First year of adoption of the 'jail' law",
    subtitle = "In the data set",
    fill = NULL
  ) +
  theme(axis.text.y = element_blank(), 
        axis.text.x = element_blank())

graph_balance <- function(balance_var) {
  fatalities |> 
    ggplot(aes(x = {{ balance_var }}, fill = jail_name, color = jail_name)) + 
    geom_density() +
    labs(
      title = paste("Balance plot for", substitute(balance_var)),
      fill = NULL,
      color = NULL,
      y = "Density"
    )
}

balance_table <- fatalities |> 
  select(-state) |> 
  modelsummary::datasummary_balance(formula = ~jail)

first_year <- fatalities |>
  filter(jail) |>
  group_by(state) |>
  mutate(treat_year = min(year, na.rm = TRUE)) |>
  ungroup() |>
  filter(year == treat_year) |>
  select(treat_year, state)

dat <- fatalities |>
  left_join(first_year, by = join_by(state)) |>
  mutate(
    first_year = replace_na(treat_year, 0)
  )

event_study_reg <- feols(
  log(fatal) ~ sunab(treat_year, year) | state + year,
  data = dat
)

event_study_graph <- event_study_reg |> 
  tidy(conf.int = TRUE) |>
  mutate(term = as.integer(str_remove(term, "year::"))) |> 
  rbind(c(-1, rep(0, 6))) |> 
  ggplot(aes(x = term, y = estimate)) + 
  geom_point() + 
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) + 
  geom_vline(xintercept = -1) + 
  labs(
    title = "Event study graph: law and impact on fatalities",
    x = "Number of years relative to the introduction of the law",
    y = "Estimate"
  )

reg_ctrl <- feols(
  data = fatalities, 
  log(fatal) ~ jail + log(unemp) + log(income) | state + year
)

coef_plot <- modelsummary::modelplot(list(`No controls` = reg, `Controls` = reg_ctrl)) +
  geom_vline(xintercept = 0) +
  labs(title = "Coefficient plot", caption = "R package: modelsummary")

distrib_coef_plot <- reg_ctrl |> 
  broom::tidy() |> 
  ggplot(aes(y = term)) +
  ggdist::stat_halfeye(
    aes(xdist = distributional::dist_student_t(
      df = df.residual(reg), mu = estimate, sigma = std.error)),
    fill = colors_mediocre$complementary,
    color = colors_mediocre$base,
    # size = 10, 
    alpha = 0.6
  ) + 
  geom_vline(xintercept = 0) +
  labs(
    title = "Distribution of estimates",
    x = "Point estimate",
    y = NULL,
    caption = "R package: ggdist"
  )

ex_duration <- readRDS("~/Documents/Teaching/data_viz_summer/content/slides/data/ex_duration.RDS")
ex_duration <- readRDS("data/ex_duration.RDS")

raw <- ex_duration |> 
  ggplot(aes(x = date, y = duration)) + 
  geom_point() + 
  scale_y_log10() +
  labs(
    title = "Evolution of the duration of items in time",
    x = NULL, 
    y = "Duration (in s)"
  )

opacity <- ex_duration |> 
  ggplot(aes(x = date, y = duration)) + 
  geom_point(alpha = 0.01) + 
  scale_y_log10() +
  labs(
    title = "Evolution of the duration of items in time, by channel",
    x = NULL, 
    y = "Duration (in s)"
  ) + 
  facet_wrap(~ channel)

heat_map <- ex_duration |> 
  ggplot(aes(x = date, y = duration)) + 
  geom_bin2d(bins = 70) + 
  scale_y_log10() +
  labs(
    title = "Evolution of the duration of items in time, by channel",
    x = NULL, 
    y = "Duration (in s)"
  ) + 
  facet_wrap(~ channel) + 
  scale_mediocre_c(pal = "coty", gradient = "left") + 
  labs(fill = "Number of observations per tile")

binscatter <- ex_duration |> 
  ggplot(aes(x = date, y = duration)) + 
  geom_point(alpha = 0.01) + 
  scale_y_log10() +
  stat_summary_bin(
    fun.y = 'mean', 
    bins = 20, 
    geom = "point", 
    color = colors_mediocre[["complementary"]],
    size = 3
  ) +
  labs(
    title = "Scatter and binscatter plot", 
    x = NULL, 
    y = "Duration (in s)"
  )