This document looks at a bunch of DDOT data that’s super cool but a little buried on their website. We’ll do 3 things:

To view the code in this document, click the “show” buttons on the right. Some of this analysis also uses supporting functions I put in a separate set of files; you can find all of the code here: https://github.com/pete-rodrigue/ddot-project-evals

knitr::opts_chunk$set(echo = TRUE, fig.width = 8, 
                      message=F, warning=F,
                      root.dir = "C:/Users/edwar/Documents/GitHub/ddot-project-evals")

setwd("C:/Users/edwar/Documents/GitHub/ddot-project-evals")
source("load_data.R")
source("make_plots.R")

# Libraries we need:
library(readr)
library(ggplot2)
library(dplyr)
library(stringr)
library(purrr)
library(broom)
library(sf)
library(leaflet)
library(tidyr)
library(tidyverse)
library(plotly)
library(zoo)
# google sheet w/ the data:
url <- "https://docs.google.com/spreadsheets/d/e/2PACX-1vQtDsTt9VVyn-wptsMT6vxtCM2J1dPYWji90OOWj1v4c4pPz6Kwa6mQtOk096EU3npvJ9g4fcjnwLpp/pub?gid=1445637369&single=true&output=csv"

# load the data from the google sheet:
df <- 
  read_csv(url, show_col_types = F) %>%
  filter(`Project name` != "15 St NW") %>%
  mutate(
    `Project short name` = 
      case_when(
          `Project name` == "16th Street NW Bus Priority Project" ~ "16 ST NW",
          `Project name` == "9th Street NW Protected Bike Lane" ~ "9 ST NW",
          `Project name` == "G Street NW Protected Bike Lane" ~ "G ST NW",
          `Project name` == "Minnesota Avenue SE Bus Priority Project" ~ "MN AVE SE",
          `Project name` == "Pennsylvania Avenue SE Bus Priority Project" ~ "PA AVE NW",
          `Project name` == "Virginia Avenue NW Protected Bike Lane" ~ "VA AVE NW",
          .default = NA
        )
  )

DDOT’s bus and bike project evaluations

First let’s look at some of the results from DDOT’s detailed before-after analyses of their projects. These all find huge reductions in crashes after project installation. The dashed lines below show the average drop in crashes by type. These projects led to either decreased driver travel times or extremely similar driving times.


# to_plot <-
  df %>%
  filter(Outcome %in% c("Injury Crashes", "All Crashes", "Pedestrian Injury")) %>%
  mutate(Change = Change*-1) %>%
  ggplot() +
  geom_bar(aes(x=`Project short name`, y=Change, fill=Outcome), stat="identity", position="dodge") +
  labs(y="% decrease in crashes", x="") +
  ggtitle("All projects saw decreases in crashes") +
  geom_hline(aes(yintercept=mean(Change[Outcome == "All Crashes"], na.rm=T)), 
             color="#F8766D", linetype="dashed", linewidth=.8, alpha=.5) +
  geom_hline(aes(yintercept=mean(Change[Outcome == "Injury Crashes"], na.rm=T)), 
             color="#00BA38", linetype="dashed", linewidth=.8, alpha=.5) +
  geom_hline(aes(yintercept=mean(Change[Outcome == "Pedestrian Injury"], na.rm=T)), 
             color="#619CFF", linetype="dashed", linewidth=.8, alpha=.5) +
  theme_minimal() +
  theme(legend.title = element_blank())

to_plot <-
  df %>%
  filter(str_detect(Outcome, "Private motor vehicle")) %>%
  group_by(`Project short name`, Outcome) %>%
  summarize(change = mean(Change, na.rm=T)) %>%
  ungroup() %>%
  mutate(Outcome = ifelse(Outcome == "Private motor vehicle AM peak travel time (seconds)", 
                          "AM traffic travel time",
                          "PM traffic travel time"))

ggplot(to_plot) +
  geom_point(
    aes(x=change, color=Outcome, y=`Project short name`), size=3, alpha=.7
  ) +
  geom_vline(aes(xintercept=0), color="black") +
  geom_vline(aes(xintercept=mean(change[Outcome == "AM traffic travel time"], na.rm=T)), color="#F8766D", linetype="dashed") +
  geom_vline(aes(xintercept=mean(change[Outcome == "PM traffic travel time"], na.rm=T)), color="#00BFC4", linetype="dashed") +
  labs(x="Change in travel time in seconds", y="") +
  theme_minimal() +
  xlim(c(min(to_plot$change)-10, max(to_plot$change)+25)) +
  theme(legend.title = element_blank()) +
  ggtitle("Drivers' travel times faster on average after project completion")

rm(to_plot); rm(url)
outcome_vars <-
  c("All Crashes", "Vehicle Crashes", "Bike Crashes", "Pedestrian Crashes", 
    "Injury Crashes", "Serious Injury Crashes", "Driver Injury", "Bicyclist Injury", "Pedestrian Injury",
    "Driver Serious Injury", "Bicyclist Serious Injury", "Pedestrian Serious Injury")


# load crash data:
cd <- load_crash_data()
# crash_details <- read_csv("data/Crash_Details_Table.csv", show_col_types = F)
# load project shapefiles:
ddot_projs <- 
  st_read("data/ddot-evaluated-projects.geojson", quiet=T) |>
  st_transform(4326) |>
  filter(project_name != "Minnesota Ave Bike Lane Project")

get_crashes_in_corridors <- function(line_sf, crashes_sf, buffer_in_meters) {
  mybuffer <-
    line_sf |>
    st_transform(6487) |>               # EPSG:6487 — NAD83 / Maryland (metres)
    st_buffer(dist = buffer_in_meters,  # X m buffer
              endCapStyle = "FLAT"
              ) |>  
    st_transform(st_crs(ddot_projs))    # reproject back to match cd's CRS
  
  crashes_near_projects <- st_join(crashes_sf, mybuffer, join = st_within, left = FALSE)
  
  crashes_near_projects
}

group_corridor_crashes_by_month <- function(df, group_by_vars = NA, outcome_vars) {
    df                                      |>
    st_drop_geometry()                      |>
    group_by(across(all_of(group_by_vars))) |>
    summarize(
      `All Crashes` = n(),
      across(all_of(outcome_vars), \(x) sum(x, na.rm = TRUE)),
      .groups = "drop_last"
      )                                     |>
    ungroup()
}

crashes_near_projects <- get_crashes_in_corridors(line_sf = ddot_projs, crashes_sf = cd, buffer_in_meters = 4)

cd_in_projs_summed <- 
  group_corridor_crashes_by_month(df = crashes_near_projects, 
                                  group_by_vars = c("project_name", "month", "year", "date"),
                                  outcome_vars = outcome_vars)
# plot the projects:
source("make_plots.R")
# create a list that has the project names and pre & post dates:
proj_dates <-
  list(
     "Minnesota Avenue SE Bus Priority Project" = 
       list(
        period_1         = list(label = "Pre",  start = "2017-02-01", end = "2020-01-31"),
        period_2         = list(label = "Post", start = "2023-11-01", end = "2024-10-31")
       ),
     "Pennsylvania Avenue SE Bus Priority Project" = 
       list(
        period_1         = list(label = "Pre",  start = "2018-06-01", end = "2022-05-31"),
        period_2         = list(label = "Post", start = "2024-05-01", end = "2024-08-31") 
       ),
     "9th Street NW Protected Bike Lane" = 
       list(
        period_1         = list(label = "Pre",  start = "2016-02-01", end = "2020-01-31"),
        period_2         = list(label = "Post", start = "2023-10-01", end = "2024-09-30") 
       ),
     "Virginia Avenue NW Protected Bike Lane" = 
       list(
        period_1         = list(label = "Pre",  start = "2018-04-01", end = "2022-03-31"),
        period_2         = list(label = "Post", start = "2023-04-01", end = "2024-03-31") 
       ),
     "21st Street NW Protected Bike Lane" = 
       list(
        period_1         = list(label = "Pre",  start = "2017-02-01", end = "2020-01-31"),
        period_2         = list(label = "Post", start = "2023-11-01", end = "2024-10-31") 
       ),
     "16th Street NW Bus Priority Project" = 
       list(
        period_1         = list(label = "Pre",  start = "2017-02-01", end = "2020-01-31"),
        period_2         = list(label = "Post", start = "2023-01-01", end = "2024-12-31") 
       ),
     "G Street NW Protected Bike Lane" = 
       list(
        period_1         = list(label = "Pre",  start = "2018-03-01", end = "2020-02-29"),
        period_2         = list(label = "Post", start = "2021-03-01", end = "2024-02-29") 
       )
  )

Below is a quick look at all of the crashes in these project corridors. Note that DDOT seems to map all of these crashes to the roadway center line. I just drew the corridor center lines by hand in Felt GIS, and then made a little buffer aroud each center line in R. I assigned all of the crashes that fell within that buffer to the corridor.

plot_leaflet(
  lines_sf = ddot_projs, 
  # polygons_sf = ddot_projs_buf, 
  line_labels = c("felt.feature", "project_name"),
  points_sf = crashes_near_projects
  )



Replicating DDOT project evaluations

Below, we try to replicate those results. The data actually agrees pretty well, given I only had DDOT’s PDF that briefly describes their methodology, not their actual code or a detailed how-to. This suggests the results are relatively robust to little choices about how to do the analysis.

Crashes go down, sometimes by quite a lot, after these projects get installed.

eval_summary <- data.frame("Project" = NA, "Source" = NA, "Outcome" = NA, "% change" = NA)
eval_objs    <- list()
for (proj in unique(cd_in_projs_summed$project_name)) {  
  
  if(proj %in% c("15th St NW/SW Protected Bike Lane", "21st Street NW Protected Bike Lane")) {next}
  
  
  periods <- proj_dates[[proj]]
  p1      <- periods$period_1
  p2      <- periods$period_2
  
  rvs <- plot_crashes_over_time(
    data             = filter(cd_in_projs_summed, year >= 2016 &
                              project_name == proj),      
    smoothing_factor = 4,
    period_1         = p1,
    period_2         = p2,
    title            = proj
  )
  
  plot_table <-
    rvs[[3]] %>% 
      dplyr::left_join(select(df, `Project name`, Outcome, Change) %>% filter(`Project name` == proj), by=c("series"="Outcome")) %>%
      rename(pct_change_ddot = Change) %>%
      select(-`Project name`) 
  
  eval_objs <- c(eval_objs, list(list(proj, rvs[[1]], plot_table)))
  
  eval_summary <- bind_rows(
    eval_summary, 
    data.frame("Project" = rep(proj, 4), 
               "Source"  = c("This analysis", "This analysis", "DDOT", "DDOT"),
               "Outcome" = rep(c("All Crashes", "Injury Crashes"), 2), 
               "% change" = c(plot_table$pct_change[plot_table$series == "All Crashes"],
                              plot_table$pct_change[plot_table$series == "Injury Crashes"],
                              plot_table$pct_change_ddot[plot_table$series == "All Crashes"],
                              plot_table$pct_change_ddot[plot_table$series == "Injury Crashes"])
               )
  )
  
}
eval_summary <- eval_summary[!is.na(eval_summary$Outcome),]
eval_summary |>
  ggplot(aes(x = X..change, y = Project, color=Source)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) +
  geom_point(size = 4.5, alpha=.6) +
  facet_wrap(~ Outcome, ncol = 1, scales = "free_y") +
  scale_x_continuous(labels = scales::label_percent(scale = 1)) +
  labs(x = "% Change", y = NULL, color = "Outcome") +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor   = element_blank(),
    strip.text         = element_text(face = "bold"),
    legend.position    = "bottom",
    legend.title       = element_blank()
  ) +
  ggtitle("We're more or less able to replicate DDOT's results")

You can view detailed tables and time series by clicking the little arrows below. The line charts are smoothed, showing a 4 month moving average.

for (myl in eval_objs) {
  cat("\n<details>\n<summary>", myl[[1]], "</summary>\n\n")
  print(myl[[2]])
  cat("\n")
  print(myl[[3]] %>% 
          rename(Outcome = series, `% change` = pct_change, `DDOT % change` = pct_change_ddot, 
                 `Pre (per month)` = Pre_per_month, `Post (per month)` = Post_per_month) %>% knitr::kable())
  cat("\n</details>\n")
}
16th Street NW Bus Priority Project

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 28.58 22.12 -22.6 -35
Injury Crashes 9.25 5.33 -42.3 NA
Serious Injury Crashes 0.36 0.42 15.4 NA
Vehicle Crashes 28.50 22.12 -22.4 NA
Bike Crashes 1.31 0.54 -58.5 -56
Pedestrian Crashes 1.61 1.38 -14.7 -20
Driver Injury 6.75 3.83 -43.2 NA
Bicyclist Injury 1.00 0.29 -70.8 -71
Pedestrian Injury 1.50 1.17 -22.2 -13
9th Street NW Protected Bike Lane

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 13.17 7.92 -39.9 -43
Injury Crashes 4.94 1.25 -74.7 NA
Serious Injury Crashes 0.21 0.08 -60.0 NA
Vehicle Crashes 13.15 7.92 -39.8 NA
Bike Crashes 0.46 0.75 63.6 -63
Pedestrian Crashes 0.88 0.50 -42.9 -50
Driver Injury 3.62 0.42 -88.5 NA
Bicyclist Injury 0.33 0.42 25.0 -77
Pedestrian Injury 0.96 0.42 -56.5 -56
G Street NW Protected Bike Lane

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 0.88 0.28 -68.3 -72
Injury Crashes 0.21 0.08 -60.0 NA
Serious Injury Crashes 0.04 0.00 -100.0 NA
Vehicle Crashes 0.88 0.28 -68.3 NA
Bike Crashes 0.08 0.03 -66.7 NA
Pedestrian Crashes 0.17 0.00 -100.0 -92
Driver Injury 0.00 0.06 Inf NA
Bicyclist Injury 0.04 0.03 -33.3 NA
Pedestrian Injury 0.17 0.00 -100.0 -86
Minnesota Avenue SE Bus Priority Project

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 14.28 8.67 -39.3 -47
Injury Crashes 6.44 2.42 -62.5 NA
Serious Injury Crashes 0.22 0.08 -62.5 NA
Vehicle Crashes 14.28 8.58 -39.9 NA
Bike Crashes 0.11 0.08 -25.0 -94
Pedestrian Crashes 0.28 0.17 -40.0 -40
Driver Injury 6.11 2.25 -63.2 NA
Bicyclist Injury 0.11 0.00 -100.0 -100
Pedestrian Injury 0.22 0.17 -25.0 -18
Pennsylvania Avenue SE Bus Priority Project

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 5.15 3.0 -41.7 -32
Injury Crashes 1.83 0.5 -72.7 NA
Serious Injury Crashes 0.02 0.0 -100.0 NA
Vehicle Crashes 5.12 3.0 -41.5 NA
Bike Crashes 0.17 0.5 200.0 NA
Pedestrian Crashes 0.33 0.0 -100.0 -100
Driver Injury 1.44 0.0 -100.0 NA
Bicyclist Injury 0.10 0.5 380.0 NA
Pedestrian Injury 0.29 0.0 -100.0 -100
Virginia Avenue NW Protected Bike Lane

Outcome Pre (per month) Post (per month) % change DDOT % change
All Crashes 2.98 2.42 -18.9 -50
Injury Crashes 1.10 0.83 -24.5 NA
Serious Injury Crashes 0.08 0.08 0.0 NA
Vehicle Crashes 2.98 2.42 -18.9 NA
Bike Crashes 0.12 0.08 -33.3 NA
Pedestrian Crashes 0.19 0.08 -55.6 -57
Driver Injury 0.81 0.75 -7.7 NA
Bicyclist Injury 0.12 0.08 -33.3 NA
Pedestrian Injury 0.15 0.00 -100.0 -100



Using additional data to expand on DDOT project evaluations

DDOT has looked at a few projects in depth. But we can also use DC’s crash data to look at every (or almost every) protected and buffered bike lane in DDOT’s “bike lane” shapefile that was built in the late 2010s and early 2020s.

Click here to read a little more about my methodology

Quick methodology:

  • Took the “Bicycle Lanes” shape file from Open Data DC and subset it to just “protected”, “dual protected”, “buffered”, and “dual buffered” bike lanes. We don’t really care about painted bike lanes, since those don’t have super robust safety benefits.
  • Unioned together all of the subblock segments in that data set by route and build year, since different pieces of some routes were built in different years. I hand-labeled when each piece of the bike lanes were installed by using satellite imagery and google street view (this was a pain). You can find my hand-labeled data set here.
  • Overlaid the crash data, and assigned any crash that’s within 4 meters of the roadway segment’s center line to the corridor (DDOT codes the crashes to the center line, rather than wherever the MPD officer poked at the map on their iPad).
  • Finally, we just plot the number of crashes per month on each route, before and after the protected or buffered bike lane was installed. I’m only showing routes that have at least 8 subblock segments and were installed between 2017 and 2024. I’m also not showing the lanes we just looked at above!
  • I’m using 2017 to 2019 as the “pre” years, unless the bike lane was installed in 2017, 18, or 19. In that case, I use 2016 up to the installation year as the “pre” years.
  • I’m using 2024 and 2025 as the “post” years, unless the bike lane was installed in 2024. In that case, I use 2025 as the “post” year.

The blue shading in the time series below shows the “before” period. The red shading shows the “after” period. I ignore peak COVID years (2020-2023) when calculating the before/after change, since those years were atypical.

Here’s a quick map of the routes. Hover over each route to see the route name, the install year, and whether it’s coded as protected (1) or buffered (0) in the Open Data DC data.

# which years to include in the analysis
ANALYSIS_YEARS = c(2016, 2017, 2018, 2019, 2024, 2025)
POST_YEARS     = c(2024, 2025)
# whether to analyze the data at the subblock level (SUBBLOCKKEY) or 
# block level (BLOCKKEY):
ANALYSIS_LEVEL = "SUBBLOCKKEY"
# which treatment groups to include, options include
# switcher     : units that got PBLs between 2020 and 2023 inclusive 
# never_taker  : units that never got PBLs
# late_taker   : units that got PBLs in 2024 or 2025
# always_taker : units that got PBLs before 2019
TREATMENT_GROUPS_TO_INCLUDE = c("switcher", "never_taker")
# whether to include buffered bike lanes as well or just fully protected bike lanes
PBLS_ONLY = FALSE
# a bbox for map plotting so we don't plot tons of data at once:
PLOT_BBOX = c(xmin = -77.020, ymin = 38.94, xmax = -77.010, ymax = 38.944)

# load data:
cd                <- load_crash_data()
pbl               <- load_pbl_data(PBLS_ONLY) 
cabi_sf           <- load_cabi_data()[[1]]
annual_cabi_rides <- load_cabi_data()[[2]]
rs                <- load_subblocks()
mm_trips <- st_read("data/DC_RideReport_MM_Trips.geojson", quiet=T) 
segment_lookup <- 
  readr::read_csv(file = "data/trips_subblocks_matched.csv", show_col_types = F)

mm_trips_tbl <-
  mm_trips                               %>%
  st_drop_geometry()                     %>%
  mutate(
    mean_2019_trips  = (trips_2019_Q1 + trips_2019_Q2 + trips_2019_Q3 + trips_2019_Q4) / 4,
    mean_2024_trips  = (trips_2024_Q1 + trips_2024_Q2 + trips_2024_Q3 + trips_2024_Q4) / 4,
    mean_2025_trips  = (trips_2025_Q1 + trips_2025_Q2 + trips_2025_Q3 + trips_2025_Q4) / 4,
    mean_24_25_trips = (mean_2024_trips + mean_2025_trips) / 2
  )                                      %>%
  select(OBJECTID, starts_with("mean_")) %>%
  replace(is.na(.), 0)


pbl <- pbl |>
  mutate(
    ROUTENAME = if_else(ROUTENAME == "Ramp-36001642", "IRVING ST NW", ROUTENAME),
    ROUTENAME = if_else(ROUTENAME == "KENYON ST NW" , "IRVING ST NW", ROUTENAME),
    ROUTENAME = if_else(ROUTENAME == "IRVING ST NE" , "IRVING ST NW", ROUTENAME),
    ROUTENAME = if_else(ROUTENAME == "PARK PL NW"   , "5TH ST NW"   , ROUTENAME),
    ) 

pbl_unioned <- pbl                                     |>
  group_by(ROUTENAME, build_year)                      |>
  summarise(
    min_install = min(build_year, na.rm = TRUE),
    max_install = max(build_year, na.rm = TRUE),
    `Has some protected lane segments` = if_else(
      max(!is.na(BIKELANE_PROTECTED) + !is.na(BIKELANE_DUAL_PROTECTED)) == 1,
          1, 0),
    n_segments  = n(),
    .groups = "drop" 
  )                                                    |>
  filter(n_segments >= 8)

bl_subblock_trips <- 
  pbl                                                  |>
  group_by(ROUTENAME, build_year)                      |>
  mutate(
    min_install = min(build_year, na.rm = TRUE),
    max_install = max(build_year, na.rm = TRUE),
    `Has some protected lane segments` = if_else(
      max(!is.na(BIKELANE_PROTECTED) + !is.na(BIKELANE_DUAL_PROTECTED)) == 1,
          1, 0),
    n_segments  = n(),
    .groups = "drop" 
  )                                                    |>
  ungroup()                                            |>
  filter(n_segments >= 8)                              |>
  select(SUBBLOCKKEY, build_year,
         max_install, ROUTENAME)                       |>
  mutate(route_unique_id =
           paste0(ROUTENAME, " (", max_install, ")"))  |>
  st_drop_geometry()                                   |>
  dplyr::left_join(segment_lookup, by="SUBBLOCKKEY")   |>
  dplyr::left_join(mm_trips_tbl, by = "OBJECTID")      
plot_leaflet(lines_sf = pbl_unioned, 
             line_labels = c("ROUTENAME", "build_year", "Has some protected lane segments"))

Note that the categorization of the bike lane types is a little unintuitive. The 17th St NW bike lane, for example, is classified as “buffered” not “protected,” which sounds less safe, but the actual bike lane is great:

The 17th St NW bike lane.
The 17th St NW bike lane.

(Photo by Martha Wilson for PopVille)

This bike lane was associated with some of the biggest safety improvements (see below).

High-level findings:

crashes_near_projects <- get_crashes_in_corridors(line_sf = pbl_unioned, crashes_sf = cd, buffer_in_meters = 4)
crashes_near_projects$route_unique_id <- paste0(crashes_near_projects$ROUTENAME, " (", crashes_near_projects$max_install, ")")


cd_in_projs_summed <- 
  group_corridor_crashes_by_month(df = crashes_near_projects, 
                                  group_by_vars = c("route_unique_id", "min_install", "max_install", "month", "year", "date", "Has some protected lane segments"),
                                  outcome_vars = outcome_vars)

already_shown <- c("21ST ST NW (2022)", "9TH ST NW (2023)", "G ST NW (2020)",  "MINNESOTA AVE SE (2023)", "PENNSYLVANIA AVE SE (2023)",  "VIRGINIA AVE NW (2022)")
bl_summary <- data.frame("Project" = NA, "Outcome" = NA, "% change" = NA, "MM change" = NA)
bl_objs    <- list()
for (proj in unique(cd_in_projs_summed$route_unique_id)) {  
  
  temp_df <- filter(cd_in_projs_summed, year >= 2016 & route_unique_id == proj)
  
  temp_df$`Serious Injury Crashes` = temp_df$`Driver Serious Injury` + temp_df$`Pedestrian Serious Injury` + temp_df$`Bicyclist Serious Injury`
  
  proj_year <- unique(temp_df$min_install)
  if(proj_year < 2017) {next}
  if (temp_df$route_unique_id[1] %in% already_shown) {next}
  
  if (proj_year < 2020) {
    p1sy = 2016; p1ey = proj_year - 1; p2sy = 2024; p2ey = 2025
  } else if (proj_year == 2024) {
    p1sy = 2017; p1ey = 2019; p2sy = 2025; p2ey = 2025
  } else {
    p1sy = 2017; p1ey = 2019; p2sy = 2024; p2ey = 2025
  }
    
  p1      <- list(label = "Pre" ,  start = paste0(p1sy, "-01-01"), end = paste0(p1ey, "-12-31"))
  p2      <- list(label = "Post",  start = paste0(p2sy, "-01-01"), end = paste0(p2ey, "-12-31"))
  
  rvs <- plot_crashes_over_time(
    data             = temp_df,      
    smoothing_factor = 4,
    period_1         = p1,
    period_2         = p2,
    title            = proj
  )
  
  plot_table <- 
    rvs[[3]] %>% 
      rename(Outcome = series,
             `Before, per month` = Pre_per_month,
             `After, per month`  = Post_per_month,
             `% change`          = pct_change)
  
  mm_data <- 
    filter(bl_subblock_trips, route_unique_id==proj)                     |> 
    summarize(mean_2019_trips  = round(mean(mean_2019_trips , na.rm=T)), 
              mean_24_25_trips = round(mean(mean_24_25_trips, na.rm=T))) |>
    mutate(mm_increase     = mean_24_25_trips - mean_2019_trips,
           mm_pct_increase = round(100*(mean_24_25_trips - mean_2019_trips) / mean_2019_trips))
  
  bl_objs <- c(bl_objs, 
               list(list("project_name"=proj, 
                         "plotly_chart"=rvs[[1]], 
                         "table"=plot_table,
                         "mm_data" = mm_data
                         )
                    )
               )
  
  bl_summary <- bind_rows(
    bl_summary, 
    data.frame("Project" = rep(proj, 2), 
               "Outcome" = c("All Crashes", "Injury Crashes"), 
               "% change" = c(plot_table$`% change`[plot_table$Outcome == "All Crashes"],
                              plot_table$`% change`[plot_table$Outcome == "Injury Crashes"])
               )
  )
}
bl_summary <- bl_summary[!is.na(bl_summary$Outcome),]
bl_summary |>
  mutate(Project = reorder(Project, -X..change)) |>
  ggplot(aes(x = X..change, y = Project, color = Outcome)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) +
  geom_point(size = 3) +
  geom_vline(
    data = bl_summary    |>
      group_by(Outcome)  |>
      summarise(mean_change = mean(X..change, na.rm = TRUE), .groups = "drop"),
    aes(xintercept = mean_change, color = Outcome),
    linetype = "solid",
    linewidth = 0.8, 
    alpha=.5
  ) +
  scale_x_continuous(labels = scales::label_percent(scale = 1)) +
  labs(x = "% Change", y = NULL, color = "Outcome") +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor   = element_blank(),
    strip.text         = element_text(face = "bold"),
    legend.position    = "bottom"
  ) +
  ggtitle("These projects lead to massive declines in crashes")

bl_summary |>
  mutate(Project = reorder(Project, -X..change)) |>
  filter(!(Project %in% c("K ST NW (2021)", "C ST NE (2023)"))) |>
  ggplot(aes(x = X..change, y = Project, color = Outcome)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) +
  geom_point(size = 3) +
  geom_vline(
    data = bl_summary    |>
      filter(!(Project %in% c("K ST NW (2021)", "C ST NE (2023)"))) |>
      group_by(Outcome)  |>
      summarise(mean_change = mean(X..change, na.rm = TRUE), .groups = "drop"),
    aes(xintercept = mean_change, color = Outcome),
    linetype = "solid",
    linewidth = 0.8, 
    alpha=.5
  ) +
  scale_x_continuous(labels = scales::label_percent(scale = 1)) +
  labs(x = "% Change", y = NULL, color = "Outcome") +
  theme_minimal(base_size = 11) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor   = element_blank(),
    strip.text         = element_text(face = "bold"),
    legend.position    = "bottom"
  ) +
  ggtitle("These projects lead to massive declines in crashes")


You can view detailed tables and time series by clicking the little arrows below. The line charts are smoothed, showing a 4 month moving average. Below the table, you can also see the change in “micromobility trips” (ebikes and scooters) from RideReport.com.


for (myl in bl_objs) {
  cat("<details>\n<summary>", myl$project_name, "</summary>")
  print(myl$plotly_chart)
  print(myl$table %>% knitr::kable())
  cat(paste0("<strong>Change in micromobility usage</strong><br>", 
               "2019 micromobility trips: ", format(myl$mm_data$mean_2019_trips, big.mark = ",") , "<br>",
               "Avg. of '24 & '25 trips:  ", format(myl$mm_data$mean_24_25_trips, big.mark = ","), "<br>",
               "Increase in trips:        ", format(myl$mm_data$mm_increase, big.mark = ",")     ,  "<br>",
               "% increase:               ", format(myl$mm_data$mm_pct_increase, big.mark = ",") ,  "<br>"
               )
        )
  cat("\n\n</details>\n")
}
14TH ST NW (2020)

Outcome Before, per month After, per month % change
All Crashes 4.06 3.42 -15.8
Injury Crashes 1.58 0.75 -52.6
Serious Injury Crashes 0.06 0.08 50.0
Vehicle Crashes 4.06 3.42 -15.8
Bike Crashes 0.53 0.62 18.4
Pedestrian Crashes 0.42 0.08 -80.0
Driver Injury 0.81 0.25 -69.0
Bicyclist Injury 0.42 0.42 0.0
Pedestrian Injury 0.36 0.08 -76.9

Change in micromobility usage
2019 micromobility trips: 586
Avg. of ’24 & ’25 trips: 9,032
Increase in trips: 8,446
% increase: 1,441

17TH ST NW (2021)

Outcome Before, per month After, per month % change
All Crashes 6.47 2.92 -54.9
Injury Crashes 2.53 0.79 -68.7
Serious Injury Crashes 0.06 0.00 -100.0
Vehicle Crashes 6.42 2.92 -54.5
Bike Crashes 0.61 0.33 -45.5
Pedestrian Crashes 0.39 0.08 -78.6
Driver Injury 1.72 0.42 -75.8
Bicyclist Injury 0.50 0.29 -41.7
Pedestrian Injury 0.31 0.08 -72.7

Change in micromobility usage
2019 micromobility trips: 366
Avg. of ’24 & ’25 trips: 3,757
Increase in trips: 3,391
% increase: 927

19TH ST SE (2024)

Outcome Before, per month After, per month % change
All Crashes 1.33 0.08 -93.8
Injury Crashes 0.47 0.08 -82.4
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 1.33 0.08 -93.8
Bike Crashes 0.03 0.00 -100.0
Pedestrian Crashes 0.06 0.00 -100.0
Driver Injury 0.42 0.08 -80.0
Bicyclist Injury 0.03 0.00 -100.0
Pedestrian Injury 0.03 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 12
Avg. of ’24 & ’25 trips: 812
Increase in trips: 800
% increase: 6,667

4TH ST NW (2020)

Outcome Before, per month After, per month % change
All Crashes 1.36 0.79 -41.8
Injury Crashes 0.53 0.12 -76.3
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 1.36 0.79 -41.8
Bike Crashes 0.08 0.08 0.0
Pedestrian Crashes 0.03 0.04 50.0
Driver Injury 0.42 0.04 -90.0
Bicyclist Injury 0.08 0.08 0.0
Pedestrian Injury 0.00 0.00 NaN

Change in micromobility usage
2019 micromobility trips: 118
Avg. of ’24 & ’25 trips: 2,997
Increase in trips: 2,879
% increase: 2,440

5TH ST NW (2022)

Outcome Before, per month After, per month % change
All Crashes 1.64 1.00 -39.0
Injury Crashes 1.08 0.17 -84.6
Serious Injury Crashes 0.06 0.00 -100.0
Vehicle Crashes 1.64 1.00 -39.0
Bike Crashes 0.08 0.04 -50.0
Pedestrian Crashes 0.08 0.04 -50.0
Driver Injury 0.94 0.17 -82.4
Bicyclist Injury 0.06 0.00 -100.0
Pedestrian Injury 0.08 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 2
Avg. of ’24 & ’25 trips: 710
Increase in trips: 708
% increase: 35,400

6TH ST NE (2017)

Outcome Before, per month After, per month % change
All Crashes 1.42 0.88 -38.2
Injury Crashes 0.50 0.25 -50.0
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 1.42 0.88 -38.2
Bike Crashes 0.08 0.17 100.0
Pedestrian Crashes 0.08 0.00 -100.0
Driver Injury 0.33 0.12 -62.5
Bicyclist Injury 0.08 0.12 50.0
Pedestrian Injury 0.08 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 42
Avg. of ’24 & ’25 trips: 3,143
Increase in trips: 3,101
% increase: 7,383

C ST NE (2023)

Outcome Before, per month After, per month % change
All Crashes 1.39 1.58 14
Injury Crashes 0.42 0.46 10
Serious Injury Crashes 0.03 0.00 -100
Vehicle Crashes 1.39 1.58 14
Bike Crashes 0.00 0.04 Inf
Pedestrian Crashes 0.00 0.00 NaN
Driver Injury 0.42 0.42 0
Bicyclist Injury 0.00 0.04 Inf
Pedestrian Injury 0.00 0.00 NaN

Change in micromobility usage
2019 micromobility trips: 0
Avg. of ’24 & ’25 trips: 1,139
Increase in trips: 1,139
% increase: Inf

COLUMBIA RD NW (2024)

Outcome Before, per month After, per month % change
All Crashes 4.50 3.50 -22.2
Injury Crashes 1.78 1.08 -39.1
Serious Injury Crashes 0.08 0.08 0.0
Vehicle Crashes 4.47 3.42 -23.6
Bike Crashes 0.67 0.33 -50.0
Pedestrian Crashes 0.58 0.50 -14.3
Driver Injury 0.69 0.50 -28.0
Bicyclist Injury 0.58 0.17 -71.4
Pedestrian Injury 0.50 0.42 -16.7

Change in micromobility usage
2019 micromobility trips: 247
Avg. of ’24 & ’25 trips: 4,301
Increase in trips: 4,054
% increase: 1,641

FLORIDA AVE NE (2019)

Outcome Before, per month After, per month % change
All Crashes 7.33 4.75 -35.2
Injury Crashes 3.50 0.96 -72.6
Serious Injury Crashes 0.19 0.00 -100.0
Vehicle Crashes 7.33 4.75 -35.2
Bike Crashes 0.25 0.42 66.7
Pedestrian Crashes 0.31 0.04 -86.4
Driver Injury 2.97 0.75 -74.8
Bicyclist Injury 0.22 0.21 -6.2
Pedestrian Injury 0.28 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 234
Avg. of ’24 & ’25 trips: 8,906
Increase in trips: 8,672
% increase: 3,706

GRANT CIR NW (2017)

Outcome Before, per month After, per month % change
All Crashes 0.33 0.17 -50
Injury Crashes 0.33 0.00 -100
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 0.33 0.17 -50
Bike Crashes 0.00 0.00 NaN
Pedestrian Crashes 0.08 0.00 -100
Driver Injury 0.17 0.00 -100
Bicyclist Injury 0.00 0.00 NaN
Pedestrian Injury 0.17 0.00 -100

Change in micromobility usage
2019 micromobility trips: 0
Avg. of ’24 & ’25 trips: 759
Increase in trips: 759
% increase: Inf

HAYES ST NE (2021)

Outcome Before, per month After, per month % change
All Crashes 1.00 0.17 -83.3
Injury Crashes 0.36 0.04 -88.5
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 1.00 0.17 -83.3
Bike Crashes 0.00 0.00 NaN
Pedestrian Crashes 0.00 0.00 NaN
Driver Injury 0.36 0.04 -88.5
Bicyclist Injury 0.00 0.00 NaN
Pedestrian Injury 0.00 0.00 NaN

Change in micromobility usage
2019 micromobility trips: 0
Avg. of ’24 & ’25 trips: 23
Increase in trips: 23
% increase: Inf

I ST SE (2024)

Outcome Before, per month After, per month % change
All Crashes 1.58 1.08 -31.6
Injury Crashes 0.97 0.50 -48.6
Serious Injury Crashes 0.06 0.00 -100.0
Vehicle Crashes 1.58 1.08 -31.6
Bike Crashes 0.14 0.00 -100.0
Pedestrian Crashes 0.08 0.00 -100.0
Driver Injury 0.81 0.50 -37.9
Bicyclist Injury 0.08 0.00 -100.0
Pedestrian Injury 0.08 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 114
Avg. of ’24 & ’25 trips: 2,451
Increase in trips: 2,337
% increase: 2,050

IRVING ST NW (2020)

Outcome Before, per month After, per month % change
All Crashes 2.94 0.83 -71.7
Injury Crashes 1.25 0.17 -86.7
Serious Injury Crashes 0.14 0.00 -100.0
Vehicle Crashes 2.94 0.83 -71.7
Bike Crashes 0.14 0.08 -40.0
Pedestrian Crashes 0.08 0.00 -100.0
Driver Injury 1.03 0.12 -87.8
Bicyclist Injury 0.14 0.04 -70.0
Pedestrian Injury 0.06 0.00 -100.0

Change in micromobility usage
2019 micromobility trips: 65
Avg. of ’24 & ’25 trips: 2,602
Increase in trips: 2,537
% increase: 3,903

K ST NW (2021)

Outcome Before, per month After, per month % change
All Crashes 2.56 2.62 2.7
Injury Crashes 0.94 1.08 14.7
Serious Injury Crashes 0.08 0.08 0.0
Vehicle Crashes 2.53 2.62 3.8
Bike Crashes 0.28 0.54 95.0
Pedestrian Crashes 0.08 0.17 100.0
Driver Injury 0.69 0.46 -34.0
Bicyclist Injury 0.17 0.50 200.0
Pedestrian Injury 0.08 0.12 50.0

Change in micromobility usage
2019 micromobility trips: 369
Avg. of ’24 & ’25 trips: 4,293
Increase in trips: 3,924
% increase: 1,063

MOUNT OLIVET RD NE (2024)

Outcome Before, per month After, per month % change
All Crashes 4.97 4.00 -19.6
Injury Crashes 1.94 1.00 -48.6
Serious Injury Crashes 0.08 0.25 200.0
Vehicle Crashes 4.94 4.00 -19.1
Bike Crashes 0.11 0.00 -100.0
Pedestrian Crashes 0.22 0.25 12.5
Driver Injury 1.69 0.75 -55.7
Bicyclist Injury 0.06 0.00 -100.0
Pedestrian Injury 0.17 0.25 50.0

Change in micromobility usage
2019 micromobility trips: 10
Avg. of ’24 & ’25 trips: 2,631
Increase in trips: 2,621
% increase: 26,210

NEW JERSEY AVE NW (2020)

Outcome Before, per month After, per month % change
All Crashes 3.56 1.96 -44.9
Injury Crashes 1.78 0.50 -71.9
Serious Injury Crashes 0.03 0.00 -100.0
Vehicle Crashes 3.53 1.96 -44.5
Bike Crashes 0.08 0.12 50.0
Pedestrian Crashes 0.11 0.08 -25.0
Driver Injury 1.61 0.33 -79.3
Bicyclist Injury 0.06 0.08 50.0
Pedestrian Injury 0.11 0.08 -25.0

Change in micromobility usage
2019 micromobility trips: 80
Avg. of ’24 & ’25 trips: 2,408
Increase in trips: 2,328
% increase: 2,910

SHERMAN CIR NW (2018)

Outcome Before, per month After, per month % change
All Crashes 0.17 0.04 -75.0
Injury Crashes 0.12 0.04 -66.7
Serious Injury Crashes 0.00 0.00 NaN
Vehicle Crashes 0.17 0.04 -75.0
Bike Crashes 0.12 0.00 -100.0
Pedestrian Crashes 0.00 0.00 NaN
Driver Injury 0.00 0.04 Inf
Bicyclist Injury 0.12 0.00 -100.0
Pedestrian Injury 0.00 0.00 NaN

Change in micromobility usage
2019 micromobility trips: 0
Avg. of ’24 & ’25 trips: 472
Increase in trips: 472
% increase: Inf

WARDER ST NW (2022)

Outcome Before, per month After, per month % change
All Crashes 1.06 0.79 -25.0
Injury Crashes 0.39 0.12 -67.9
Serious Injury Crashes 0.03 0.04 50.0
Vehicle Crashes 1.06 0.79 -25.0
Bike Crashes 0.03 0.08 200.0
Pedestrian Crashes 0.00 0.04 Inf
Driver Injury 0.36 0.04 -88.5
Bicyclist Injury 0.03 0.04 50.0
Pedestrian Injury 0.00 0.04 Inf

Change in micromobility usage
2019 micromobility trips: 0
Avg. of ’24 & ’25 trips: 1,284
Increase in trips: 1,284
% increase: Inf

WEST VIRGINIA AVE NE (2021)

Outcome Before, per month After, per month % change
All Crashes 4.03 2.79 -30.7
Injury Crashes 1.50 0.67 -55.6
Serious Injury Crashes 0.06 0.12 125.0
Vehicle Crashes 4.03 2.79 -30.7
Bike Crashes 0.25 0.17 -33.3
Pedestrian Crashes 0.08 0.08 0.0
Driver Injury 1.25 0.38 -70.0
Bicyclist Injury 0.14 0.21 50.0
Pedestrian Injury 0.11 0.08 -25.0

Change in micromobility usage
2019 micromobility trips: 128
Avg. of ’24 & ’25 trips: 5,480
Increase in trips: 5,352
% increase: 4,181



Comparing streets with these new projects against other DC streets

Maybe the decline in crashes was not particular to these corridors, and in fact happened across DC during this time? Maybe it’s all the new speed cameras, or the fact DC is towing chronic speeders? The declines we observe on bike lane streets could just be part of a bigger, District-wide trend.

To try and answer that question, let’s compare these streets with new bike lane projects to other streets in DC, and see if they have similar or different crash trends.

end_of_treatment = 2023
min_year         = 2017

rs$length_m <- as.numeric(st_length(rs))

crashes_near_subblocks <- 
  get_crashes_in_corridors(line_sf = rs, crashes_sf = cd, buffer_in_meters = 4) |>
  st_drop_geometry() |>
    # For crashes in multiple buffers, randomly keep one. this will add noise but if anything that will bias our effect sizes downwards.
    group_by(CRIMEID) |>
    slice_sample(n = 1) |>         # randomly pick one buffer per crash
    ungroup() |>
  # group_by(SUBBLOCKKEY, year, month, date) |>
  group_by(SUBBLOCKKEY, year) |>
  summarise(
    across(ends_with(" Crashes"), \(x) sum(x, na.rm = TRUE)),
    across(ends_with(" Injury"), \(x) sum(x, na.rm = TRUE)),
      `Total Crashes`      = n(),
      TOTALTRAVELLANES     = first(TOTALTRAVELLANES),
      TOTALPARKINGLANES    = first(TOTALPARKINGLANES),
      TOTALTRAVELLANEWIDTH = first(TOTALTRAVELLANEWIDTH),
      WARD_ID              = first(WARD_ID),
      DCFUNCTIONALCLASS    = first(DCFUNCTIONALCLASS),
    .groups   = "drop"
  ) |>
  ungroup() |>
  dplyr::left_join(
  st_drop_geometry(pbl) %>% 
    mutate(
           pbl = ifelse((!is.na(BIKELANE_PROTECTED)) | (!is.na(BIKELANE_DUAL_PROTECTED)), 1, 0),
           bbl = ifelse(is.na(BIKELANE_PROTECTED)    &  is.na(BIKELANE_DUAL_PROTECTED), 1, 0)
           ) %>%
    rename(route_name = ROUTENAME) %>%
    mutate(unique_route_id = paste(route_name, build_year)) %>%
    select(SUBBLOCKKEY, build_year, pbl, bbl, route_name, unique_route_id), 
  by = "SUBBLOCKKEY") |>
  dplyr::left_join(st_drop_geometry(rs) %>% select(SUBBLOCKKEY, length_m), by="SUBBLOCKKEY")


# Identify variable groups programmatically
crash_vars         <- names(crashes_near_subblocks)[
                        grepl("Crashes$|Injury$", names(crashes_near_subblocks))
                      ]
time_invariant_vars <- c("TOTALTRAVELLANES", "TOTALPARKINGLANES", "TOTALTRAVELLANEWIDTH",
                          "WARD_ID", "DCFUNCTIONALCLASS", "build_year", "pbl", "bbl",
                          "route_name", "unique_route_id", "length_m")



crashes_near_subblocks_filled <- crashes_near_subblocks |>
  filter(year >= min_year) |>
  # Complete over all SUBBLOCKKEY × existing year-month combinations
  # using nesting() so we don't invent year-month combos that don't exist in the data
  complete(
    SUBBLOCKKEY, 
    # nesting(year, month, date),
    # nesting(year),
    year = min_year:2025,
    fill = setNames(as.list(rep(0, length(crash_vars))), crash_vars)
  ) |>
  # Fill time-invariant characteristics from whichever rows have data for that subblock
  group_by(SUBBLOCKKEY) |>
  fill(all_of(time_invariant_vars), .direction = "downup") |>
  mutate(
    mean_crashes_per_800_meter_pre = mean(
      `Total Crashes`[year %in% seq(min_year, 2017, 1)] / length_m[year %in% seq(min_year, 2017, 1)] * 800,
      na.rm = TRUE
    )
  ) |>
  ungroup() |>
  mutate(
    mean_crashes_per_meter_pre_bin = ntile(mean_crashes_per_800_meter_pre, 3)
  ) |>
  mutate(
    across(
      all_of(crash_vars),
      \(x) x / length_m * 800,
      .names = "{.col} per 800 meters"
    )
  ) |>
  filter(build_year %in% 2020:end_of_treatment | is.na(build_year)) |>
  filter(year %in% c(min_year:2019, (end_of_treatment+1):2025)) |>
  mutate(
    treated = if_else(is.na(build_year), 0, 1),
    post    = if_else(year > end_of_treatment, 1, 0)
  )

Here are the changes in means across different variables. The bike lane streets start from a higher level of crashes and, with the exception of bike crashes, see a bigger drop in crashes. (This is likely due to a large increase in bike traffic).

crash_vars_to_plot <- c("Total Crashes per 800 meters", "Injury Crashes per 800 meters", 
                         "Vehicle Crashes per 800 meters", "Driver Injury per 800 meters",
                         "Pedestrian Crashes per 800 meters", "Pedestrian Injury per 800 meters",
                        "Bike Crashes per 800 meters", "Bicyclist Injury per 800 meters")

panel_pairs <- list(
  "Overall"     = list(total = "Total Crashes per 800 meters",    injury = "Injury Crashes per 800 meters"),
  "People driving"      = list(total = "Vehicle Crashes per 800 meters",  injury = "Driver Injury per 800 meters"),
  "People walking"  = list(total = "Pedestrian Crashes per 800 meters", injury = "Pedestrian Injury per 800 meters"),
  "People bicycling"        = list(total = "Bike Crashes per 800 meters",     injury = "Bicyclist Injury per 800 meters")
)

plot_data <- map_dfr(names(panel_pairs), \(panel) {
  pair <- panel_pairs[[panel]]
  crashes_near_subblocks_filled |>
    select(treated, post, total = all_of(pair$total), injury = all_of(pair$injury)) |>
    mutate(panel = panel)
}) |>
  mutate(
    treated = if_else(treated == 1, "Bike lane installs", "Other streets"),
    post    = factor(if_else(post == 1, "Post", "Pre"), levels = c("Pre", "Post")),
    panel   = factor(panel, levels = names(panel_pairs))
  ) |>
  group_by(treated, post, panel) |>
  summarise(
    total  = mean(total,  na.rm = TRUE),
    injury = mean(injury, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(plot_data, aes(x = post, group = treated, color = treated)) +
  geom_line(aes(y = total), linewidth = 0.8, alpha = 0.4, linetype = "dashed") +
  geom_point(aes(y = total), size = 3, alpha = 0.4) +
  geom_line(aes(y = injury), linewidth = 0.8) +
  geom_point(aes(y = injury), size = 3) +
  facet_wrap(~ panel, 
             scales = "free_y",
             ncol = 4) +
  scale_color_manual(values = c("Other streets" = "grey50", "Bike lane installs" = "#4575b4")) +
  scale_y_continuous(expand = expansion(mult = 0.2)) +
  labs(
    x        = NULL,
    y        = "Mean annual crashes per subblock",
    color    = NULL,
    title    = "Pre vs. post crash rates by mode",
    subtitle = "Solid = injury crashes, faded/dashed = total crashes"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    legend.position  = "bottom",
    panel.grid.minor = element_blank(),
    strip.text       = element_text(face = "bold")
  )

We can conduct permutation tests to get a sense for the statistical significance of these differences in means. Did streets with bike lanes see bigger declines in crashes, other than bike crashes, than other streets? The answer seems to be yes. Note that the statistical significance of the “bike crashes increased” results are weaker (not significant at the 10% level).

# ── Observed DiD statistic ─────────────────────────────────────────────────────
compute_did <- function(df, outcome) {
  df |>
    group_by(treated, post) |>
    summarise(mean_val = mean(.data[[outcome]], na.rm = TRUE), .groups = "drop") |>
    pivot_wider(names_from = post, values_from = mean_val, names_prefix = "post_") |>
    mutate(change = post_1 - post_0) |>
    summarise(did = diff(change)) |>   # treated change minus control change
    pull(did)
}

# ── Pre-compute once ───────────────────────────────────────────────────────────
subblock_treatment <- crashes_near_subblocks_filled |> distinct(SUBBLOCKKEY, treated)
df_no_treated      <- crashes_near_subblocks_filled |> select(-treated)

# ── Permutation test function ──────────────────────────────────────────────────
run_permutation_test <- function(outcome_var, df, df_with_treated, lookup, 
                                 n_permutations = 500,
                                 alternative = "less") {  # "less" or "greater"
  
  observed_did <- compute_did(df_with_treated, outcome_var)
  
  permuted_dids <- map_dbl(1:n_permutations, \(i) {
    permuted <- lookup |> mutate(treated = sample(treated))
    df |>
      left_join(permuted, by = "SUBBLOCKKEY") |>
      compute_did(outcome_var)
  })
  
  # One-sided p-value
  p_value <- if (alternative == "less") {
    mean(permuted_dids <= observed_did)   # H1: DiD < 0 (crashes went down)
  } else {
    mean(permuted_dids >= observed_did)   # H1: DiD > 0 (bike crashes went up)
  }
  
  clean_label <- str_remove(outcome_var, " per 800 meters")
  
  tibble(did = permuted_dids) |>
    ggplot(aes(x = did)) +
    geom_histogram(bins = 50, fill = "grey70", color = "white") +
    geom_vline(xintercept = observed_did, color = "#542788", linewidth = 1) +
    labs(
      x        = "Permuted DiD",
      y        = "Count",
      title    = clean_label,
      subtitle = sprintf("DiD = %.4f, p = %.3f (%s)", observed_did, p_value, alternative)
    ) +
    theme_minimal(base_size = 9) +
    theme(
      plot.title    = element_text(face = "bold", size = 9),
      plot.subtitle = element_text(size = 8)
    )
}

# ── Run with appropriate one-sided test per outcome ────────────────────────────
bike_vars  <- c("Bike Crashes per 800 meters", "Bicyclist Injury per 800 meters")

n = 500
plots <- map(crash_vars_to_plot, \(v) {
  run_permutation_test(
    outcome_var     = v,
    df              = df_no_treated,
    df_with_treated = crashes_near_subblocks_filled,
    lookup          = subblock_treatment,
    n_permutations  = n,
    alternative     = if_else(v %in% bike_vars, "greater", "less")
  )
})

patchwork::wrap_plots(plots, ncol = 2) +
  patchwork::plot_annotation(
    title    = "Permutation tests: DiD estimates (one-sided)",
    subtitle = sprintf("%d permutations per outcome", n)
  )

What’s going on with people crashing on bicycles on the bike lane streets? Part of the answer seems to be that there are just way, way more people biking on these streets than there were previously. The increase is much larger on streets that got bike lanes than on other streets.

Here are the micromobility trip totals for all road segments in DC:

to_plot <-
  mm_trips %>%
  st_drop_geometry() %>%
  summarise(across(starts_with("trips_"), sum, na.rm = TRUE)) %>%
  tidyr::pivot_longer(col = everything(), names_to = "period", values_to = "trips") %>%
    mutate(
    year    = as.integer(stringr::str_extract(period, "\\d{4}")),
    quarter = as.integer(stringr::str_extract(period, "(?<=Q)\\d")),
    time    = year + (quarter - 1) / 4,
    label   = paste0(year, " Q", quarter)
  )

ggplot(to_plot, aes(x = time, y = trips)) +
  geom_line(linewidth = 0.8, color="#58933e") +
  geom_point(size = 2.5, color="#58933e") +
  scale_x_continuous(breaks = to_plot$time[seq(1, nrow(to_plot), 4)], labels = to_plot$label[seq(1, nrow(to_plot), 4)]) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x=NULL, y = "Trips") +
  ggtitle("Micromobility trips up")

But as the chart below shows, rides have increased way more on streets that had bike lanes installed:

crashes_near_subblocks_filled                                                  |>
  dplyr::left_join(segment_lookup, by="SUBBLOCKKEY")                           |>
  dplyr::left_join(mm_trips_tbl, by = "OBJECTID")                              |> 
  group_by(treated)                                                            |>
  summarize(
    mean_2019_trips = mean(mean_2019_trips, na.rm=T),
    mean_24_25_trips = mean(mean_24_25_trips, na.rm=T)
  )                                                                            |>
  tidyr::pivot_longer(cols = -treated, names_to = "post", values_to = "trips") |>
  mutate(post = if_else(post=="mean_2019_trips", 0, 1))                        |>
  
  mutate(
    treated = if_else(treated == 1, "Bike lane streets", "Other streets"),
    post    = factor(if_else(post == 1, "After", "Before"), levels = c("Before", "After"))
  ) |>
  ggplot(aes(x = post, y = trips, fill = treated)) +
  geom_col(position = "dodge", alpha = 0.8) +
  geom_text(aes(label = round(trips, 0)),
            position = position_dodge(width = 0.9),
            vjust = -0.5, size = 3.5) +
  scale_fill_manual(values = c("Other streets" = "grey50", "Bike lane streets" = "#4575b4")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    x    = NULL,
    y    = "Trips",
    fill = NULL,
    title = "Trips before and after bike lane installation"
  ) +
  theme_minimal(base_size = 11) +
  theme(
    legend.position  = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank()
  )



Click here to see some empirical CDFs

Here are the empirical cumulative distribution functions:

# crashes_near_subblocks_filled |>
#   mutate(treated = if_else(is.na(build_year), "Other streets", "Bike lane streets")) |>
#   ggplot(aes(x = log(`Total Crashes per 800 meters`+.0001), fill = treated)) +
#   geom_density(binwidth = 1, alpha = 0.4, position = "identity") +
#   scale_fill_manual(values = c("Other streets" = "#b35806", "Bike lane streets" = "#542788")) 

# crashes_near_subblocks_filled |>
#   mutate(treated = if_else(is.na(build_year), "Other streets", "Bike lane streets")) |>
#   mutate(post    = if_else(year >= 2024, "2. post", "1. pre")) |>
#   ggplot(aes(x = log(`Total Crashes per 800 meters`+.0001), fill = treated)) +
#   geom_density(binwidth = 1, alpha = 0.4, position = "identity") +
#   scale_fill_manual(values = c("Other streets" = "#b35806", "Bike lane streets" = "#542788")) +
#   facet_wrap(~post, ncol=1)


ecdf <-
  crashes_near_subblocks_filled |>
  mutate(treated = if_else(is.na(build_year), "Other streets", "Bike lane streets")) |>
  ggplot(aes(x = `Total Crashes per 800 meters`, color = treated)) +
  stat_ecdf(geom = "step") +
  scale_color_manual(values = c("Other streets" = "#b35806", "Bike lane streets" = "#542788")) +
  labs(color = NULL) +
  theme(legend.position = "top") +
  ggtitle("Empirical CDF by group, all years pooled")

ggplotly(ecdf) %>% layout(hovermode = "x")
ecdf <-
  crashes_near_subblocks_filled |>
  mutate(treated = if_else(is.na(build_year), "Other streets", "Bike lane streets")) |>
  mutate(post    = if_else(post == 1, "2. After", "1. Before")) |>
  ggplot(aes(x = `Total Crashes per 800 meters`, color = treated)) +
  stat_ecdf(geom = "step") +
  scale_color_manual(values = c("Other streets" = "#b35806", "Bike lane streets" = "#542788")) +
  labs(color = NULL) +
  facet_wrap(~post) +
  ggtitle("Empirical CDF by group, before vs. after")

ggplotly(ecdf) %>% layout(hovermode = "x")
# 
# ecdf <-
#   crashes_near_subblocks_filled |>
#   mutate(treated = if_else(is.na(build_year), "Other streets", "Bike lane streets")) |>
#   mutate(post    = if_else(post == 1, "2. After", "1. Before")) |>
#   ggplot(aes(x = `Injury Crashes per 800 meters`, color = treated)) +
#   stat_ecdf(geom = "step") +
#   scale_color_manual(values = c("Other streets" = "#b35806", "Bike lane streets" = "#542788")) +
#   labs(color = NULL) +
#   facet_wrap(~post) +
#   ggtitle("Empirical CDF by group, before vs. after")
# 
# ggplotly(ecdf) %>% layout(hovermode = "x")