How the model works

The model answers two questions:

  1. Which counted locations have more pedestrian crashes than their foot traffic would predict? Crashes are assigned to the nearest pedestrian count site. A negative binomial safety performance function (SPF) estimates the expected crashes for a site with that pedestrian volume, and an Empirical Bayes (EB) adjustment (Highway Safety Manual, Part B) blends the prediction with what was observed so a single random crash does not dominate the ranking. Sites are ranked by excess expected crashes (EB estimate minus SPF prediction).
  2. Where do pedestrian crashes cluster, counted or not? All pedestrian crashes are binned into a hexagon grid, weighted by severity when a severity field is available, and tested with the Getis-Ord Gi* statistic. Adjacent significant cells are merged into hot spot zones.

Count records that share a station number but differ only by a direction suffix (e.g. 50-00013N_EB, 50-00013N_NB, …) are treated as legs of one site and summed.

Everything you may need to change is in the configuration chunk below.

0. Configuration

CONFIG <- list(
  # --- Input files --------------------------------------------------------
  count_file    = "*.xlsxCombined_Files.csv",  # if not found, any file ending in "xlsxCombined_Files.csv" is used
  location_file = "MetrocLoc_Combined.csv",    # count-site coordinates: Location, Latitude, Longitude
  crash_file    = "Metrocrashes.csv",          # crash export (Michigan Traffic Crash Facts layout)
  factor_file   = "mdot_factor_table.csv",     # OPTIONAL: CountMonth, DayOfWeek, Seasonal_Factor, DOW_Factor

  # --- Pedestrian counts --------------------------------------------------
  ped_mode_pattern = "^ped",                       # Mode values in the file: Ped, Bike, Escooter, Unclassified, Total
  leg_suffix_regex = "_(N|S|E|W|NE|NW|SE|SW)B$",   # "82-00097N_WB" -> site "82-00097N"
  default_district = "University",                 # last-resort fill for a blank District

  # --- Crash columns: NULL = auto-detect, or type the exact column name ---
  # (read.csv turns "Crash: Pedestrian" into "Crash..Pedestrian", "Crash Year" into "Crash.Year", etc.)
  crash_cols = list(lat = NULL, lon = NULL, ped_flag = NULL, severity = NULL, year = NULL,
                    date = NULL, id = NULL, speed = NULL, intersection = NULL, city = NULL, region = NULL),
  crashes_are_ped_only = FALSE,     # TRUE if the crash file is already filtered to pedestrian crashes
  study_years = NULL,               # NULL = derived from the crash years

  # --- Spatial settings ---------------------------------------------------
  crs_meters    = 3078,   # Michigan GeoRef, units = meters
  match_dist_m  = 76,     # crash -> nearest count site if within 76 m (250 ft, HSM intersection influence area)
  hex_size_m    = 400,    # hot spot grid cell size (edge to edge)
  hex_context_m = 1200,   # grid cells farther than this from any crash are left out of the hot spot test
  limit_sites_to_crash_area = TRUE,  # only model count sites inside the regions/area the crash file covers

  # --- Severity and speed -------------------------------------------------
  # KABCO weights (relative). Used only if the crash file has a severity column such as "Worst Injury in Crash".
  severity_weights = c(K = 10, A = 10, B = 4, C = 2, O = 1, Unknown = 1),
  high_speed_mph   = 40,  # crashes at or above this speed limit are flagged as high-speed

  # --- Tiering --------------------------------------------------------------
  high_tier_share = 0.10,  # sites in the top 10% by excess crashes (and above expected) are "High"

  out_dir = "outputs"
)
dir.create(CONFIG$out_dir, showWarnings = FALSE)

1. Helper functions

`%||%` <- function(a, b) if (is.null(a)) b else a   # built into R >= 4.4; defined here for older R

# Use the configured file, or fall back to the first file matching a pattern
resolve_file <- function(path, fallback_pattern = NULL) {
  if (file.exists(path)) return(path)
  if (!is.null(fallback_pattern)) {
    hits <- list.files(pattern = fallback_pattern)
    if (length(hits) > 0) {
      message(sprintf("'%s' not found; using '%s'.", path, hits[1]))
      return(hits[1])
    }
  }
  stop(sprintf("File not found: '%s' (working directory: %s)", path, getwd()))
}

# Numeric conversion that tolerates "1,234" and blanks
to_num <- function(x) suppressWarnings(as.numeric(gsub(",", "", as.character(x))))

# Most common non-missing value
mode_of <- function(x) {
  x <- x[!is.na(x)]
  if (length(x) == 0) NA_character_ else names(sort(table(x), decreasing = TRUE))[1]
}

# Return a column name: the override if supplied, otherwise the first regex that matches
find_col <- function(df, override, patterns, key, required = TRUE) {
  if (!is.null(override)) {
    if (!override %in% names(df)) stop(sprintf("CONFIG$crash_cols$%s = '%s' is not a column in the crash file.", key, override))
    return(override)
  }
  for (p in patterns) {
    hit <- grep(p, names(df), ignore.case = TRUE, value = TRUE)
    if (length(hit) > 0) return(hit[1])
  }
  if (required) stop(sprintf(
    "Could not auto-detect the '%s' column. Set CONFIG$crash_cols$%s. Columns available: %s",
    key, key, paste(names(df), collapse = ", ")))
  NULL
}

# Keep points inside Michigan; fix longitudes recorded without the minus sign
clean_coords <- function(df, what) {
  n0 <- nrow(df)
  df <- df %>%
    mutate(Latitude  = to_num(Latitude),
           Longitude = to_num(Longitude),
           Longitude = if_else(Longitude > 0, -Longitude, Longitude)) %>%
    filter(between(Latitude, 41.5, 48.5), between(Longitude, -90.5, -82))
  message(sprintf("%s: kept %d of %d rows with valid Michigan coordinates.", what, nrow(df), n0))
  if (n0 > 0 && nrow(df) / n0 < 0.5)
    warning(sprintf("%s: over half the rows were dropped. Are the coordinates projected (e.g. State Plane) instead of lat/long?", what))
  df
}

# Interpret a yes/no style flag ("Pedestrian Involved" / "No Pedestrian Involved", Y/N, 1/0, TRUE/FALSE)
is_yes <- function(x) {
  if (is.logical(x)) return(coalesce(x, FALSE))
  num <- to_num(x)
  txt <- str_to_lower(str_trim(as.character(x)))
  out <- ifelse(!is.na(num), num > 0,
                txt %in% c("y", "yes", "true", "t", "x") |
                  (str_detect(txt, "pedestrian|involved") & !str_detect(txt, "^no|not |non|none")))
  coalesce(out, FALSE)
}

# Map injury text to KABCO
parse_kabco <- function(x) {
  x <- str_to_upper(str_trim(as.character(x)))
  case_when(
    str_detect(x, "\\(K\\)|^K$|FATAL")                        ~ "K",
    str_detect(x, "\\(B\\)|^B$|NON-?INCAP|MINOR")             ~ "B",  # checked before A: "non-incapacitating"
    str_detect(x, "\\(A\\)|^A$|SERIOUS|INCAPACITATING")       ~ "A",
    str_detect(x, "\\(C\\)|^C$|POSSIBLE")                     ~ "C",
    str_detect(x, "\\(O\\)|^O$|NO INJURY|PDO|PROPERTY")       ~ "O",
    TRUE                                                      ~ "Unknown"
  )
}

# Getis-Ord Gi* z-scores with binary weights; nb must include each cell itself
gi_star <- function(x, nb) {
  n    <- length(x)
  xbar <- mean(x)
  s    <- sqrt(sum(x^2) / n - xbar^2)
  if (!is.finite(s) || s == 0) return(rep(0, n))
  w    <- lengths(nb)
  lag  <- vapply(nb, function(j) sum(x[j]), numeric(1))
  (lag - xbar * w) / (s * sqrt((n * w - w^2) / (n - 1)))
}

2. Pedestrian counts

2.1 Load and de-duplicate

counts_raw <- read.csv(resolve_file(CONFIG$count_file, "xlsxCombined_Files\\.csv$"),
                       stringsAsFactors = FALSE, na.strings = c("", "NA", "-"))

# Make sure descriptive columns exist so later steps never fail on a missing field
for (col in c("Community", "CountyName", "District", "OnRoad", "Pathway", "LocationDescription")) {
  if (!col %in% names(counts_raw)) counts_raw[[col]] <- NA_character_
}

counts_raw <- counts_raw %>%
  mutate(across(where(is.character), str_trim),
         across(c(Community, CountyName, District, OnRoad, Pathway, LocationDescription), as.character),
         MasterLocalId = as.character(MasterLocalId),
         Daily_Total   = to_num(Daily_Total))

counts_valid <- filter(counts_raw, !is.na(MasterLocalId), !is.na(Mode))
counts       <- distinct(counts_valid)

tibble(raw_rows             = nrow(counts_raw),
       blank_id_or_mode     = nrow(counts_raw) - nrow(counts_valid),
       exact_duplicates     = nrow(counts_valid) - nrow(counts),
       unique_rows          = nrow(counts),
       fan_out_factor       = round(nrow(counts_valid) / nrow(counts), 1))
## # A tibble: 1 × 5
##   raw_rows blank_id_or_mode exact_duplicates unique_rows fan_out_factor
##      <int>            <int>            <int>       <int>          <dbl>
## 1   349151             3540           336917        8694           39.8

2.2 QA checks

count(counts, Mode, name = "rows")
##           Mode rows
## 1         Bike 1889
## 2     Escooter  408
## 3          Ped 2202
## 4        Total 3011
## 5 Unclassified 1184
# Same site, day, pathway and mode recorded more than once with different values (summed later)
counts %>%
  count(MasterLocalId, StartDateTime, Pathway, Mode, name = "n") %>%
  summarise(repeated_records = sum(n > 1))
##   repeated_records
## 1              364
# Do the "Total" rows equal the sum of the individual modes?
group_cols <- setdiff(names(counts), c("Mode", "Daily_Total"))

counts %>%
  filter(Mode != "Total") %>%
  group_by(across(all_of(group_cols))) %>%
  summarise(sum_of_modes = sum(Daily_Total, na.rm = TRUE), .groups = "drop") %>%
  inner_join(counts %>%
               filter(Mode == "Total") %>%
               select(all_of(group_cols), reported_total = Daily_Total),
             by = group_cols, relationship = "many-to-many") %>%
  mutate(diff = reported_total - sum_of_modes) %>%
  summarise(records_checked = n(),
            exact_match     = sum(diff == 0, na.rm = TRUE),
            mismatched      = sum(diff != 0, na.rm = TRUE),
            missing_values  = sum(is.na(diff)))
## # A tibble: 1 × 4
##   records_checked exact_match mismatched missing_values
##             <int>       <int>      <int>          <int>
## 1            3011        2693        318              0

2.3 Keep pedestrian counts and build site-level daily volumes

ped_counts <- counts %>%
  filter(str_detect(Mode, regex(CONFIG$ped_mode_pattern, ignore_case = TRUE)))

if (nrow(ped_counts) == 0) {
  warning("No Mode matched CONFIG$ped_mode_pattern; falling back to Mode == 'Total' (this includes bicycles).")
  ped_counts <- counts %>% filter(Mode == "Total")
}

# Fill blank District from the county, then from the county code at the start of the ID
known <- counts %>%
  filter(!is.na(District)) %>%
  mutate(prefix = str_extract(MasterLocalId, "^\\d+"))
county_lookup <- known %>%
  filter(!is.na(CountyName)) %>%
  count(CountyName, District) %>%
  group_by(CountyName) %>% slice_max(n, n = 1, with_ties = FALSE) %>% ungroup() %>%
  select(CountyName, district_by_county = District)
prefix_lookup <- known %>%
  filter(!is.na(prefix)) %>%
  count(prefix, District) %>%
  group_by(prefix) %>% slice_max(n, n = 1, with_ties = FALSE) %>% ungroup() %>%
  select(prefix, district_by_prefix = District)

n_blank_district <- sum(is.na(ped_counts$District))

ped_counts <- ped_counts %>%
  mutate(prefix = str_extract(MasterLocalId, "^\\d+")) %>%
  left_join(county_lookup, by = "CountyName") %>%
  left_join(prefix_lookup, by = "prefix") %>%
  mutate(
    District  = coalesce(District, district_by_county, district_by_prefix, CONFIG$default_district),
    site_id   = str_remove(MasterLocalId, CONFIG$leg_suffix_regex),
    CountDate = as_date(parse_date_time(StartDateTime,
                                        orders = c("Ymd HMS", "Ymd HM", "Ymd", "mdY HMS", "mdY IMSp", "mdY HM", "mdY"),
                                        quiet = TRUE)),
    Pathway_Type = case_when(
      str_detect(Pathway, regex("crosswalk", ignore_case = TRUE)) ~ "Crosswalk",
      str_detect(Pathway, regex("bike lane", ignore_case = TRUE)) ~ "Bike Lane",
      str_detect(Pathway, regex("roadway",   ignore_case = TRUE)) ~ "Roadway",
      str_detect(Pathway, regex("sidewalk",  ignore_case = TRUE)) ~ "Sidewalk",
      TRUE ~ "Other / unspecified")
  ) %>%
  select(-prefix, -district_by_county, -district_by_prefix) %>%
  filter(!is.na(CountDate), !is.na(Daily_Total), Daily_Total >= 0)

message(sprintf("Filled District for %d pedestrian rows that were blank.", n_blank_district))
## Filled District for 134 pedestrian rows that were blank.
# One row per site per day: sum across legs, directions and pathways counted at that site
site_days <- ped_counts %>%
  group_by(site_id, CountDate) %>%
  summarise(Community   = first(Community),
            CountyName  = first(CountyName),
            District    = first(District),
            OnRoad      = paste(unique(na.omit(OnRoad)), collapse = " / "),
            Description = first(LocationDescription),
            leg_ids     = paste(sort(unique(MasterLocalId)), collapse = "; "),
            Daily_Peds  = sum(Daily_Total),
            .groups = "drop") %>%
  mutate(CountMonth = as.character(month(CountDate, label = TRUE, abbr = TRUE)),
         DayOfWeek  = as.character(wday(CountDate,  label = TRUE, abbr = TRUE)))

# Optional seasonal / day-of-week adjustment
if (file.exists(CONFIG$factor_file)) {
  factors <- read.csv(CONFIG$factor_file, stringsAsFactors = FALSE) %>%
    mutate(CountMonth = as.character(CountMonth), DayOfWeek = as.character(DayOfWeek))
  site_days <- left_join(site_days, factors, by = c("CountMonth", "DayOfWeek"))
  message("Applied adjustment factors from ", CONFIG$factor_file)
} else {
  message("No factor table found (", CONFIG$factor_file, "); volumes are unadjusted average daily counts.")
  site_days <- mutate(site_days, Seasonal_Factor = 1, DOW_Factor = 1)
}
## No factor table found (mdot_factor_table.csv); volumes are unadjusted average daily counts.
site_days <- site_days %>%
  mutate(across(c(Seasonal_Factor, DOW_Factor), ~ coalesce(as.numeric(.x), 1)),
         Adjusted_Peds = Daily_Peds * Seasonal_Factor * DOW_Factor)

site_volume <- site_days %>%
  group_by(site_id) %>%
  summarise(Community    = first(Community),
            CountyName   = first(CountyName),
            District     = first(District),
            OnRoad       = mode_of(OnRoad),
            Description  = first(Description),
            leg_ids      = paste(sort(unique(unlist(str_split(leg_ids, "; ")))), collapse = "; "),
            ped_volume   = mean(Adjusted_Peds),   # average daily pedestrians
            days_counted = n(),
            first_count  = min(CountDate),
            last_count   = max(CountDate),
            .groups = "drop")

tibble(count_sites = nrow(site_volume),
       site_days   = nrow(site_days),
       sites_counted_one_day = sum(site_volume$days_counted == 1))
## # A tibble: 1 × 3
##   count_sites site_days sites_counted_one_day
##         <int>     <int>                 <int>
## 1         380      1904                   220
count(site_volume, District, name = "sites")
## # A tibble: 7 × 2
##   District   sites
##   <chr>      <int>
## 1 Bay           27
## 2 Grand         59
## 3 Metro         86
## 4 North         61
## 5 Southwest     23
## 6 Superior      20
## 7 University   104

2.4 Exploratory summaries

ped_counts %>%
  group_by(District, Pathway_Type) %>%
  summarise(avg_daily = mean(Daily_Total), .groups = "drop") %>%
  ggplot(aes(District, avg_daily)) +
  geom_col(fill = "steelblue") +
  facet_wrap(~ Pathway_Type, scales = "free_y") +
  labs(title = "Average daily pedestrian count by region and pathway type",
       subtitle = "Mean per count record",
       x = "Region", y = "Average daily pedestrians") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(site_days, aes(District, Daily_Peds, fill = District)) +
  geom_boxplot(outlier.color = "red", outlier.size = 2, alpha = 0.7) +
  scale_y_continuous(trans = "log1p", breaks = c(0, 10, 100, 1000, 10000), labels = scales::comma) +
  labs(title = "Daily pedestrian volume by region (QA/QC)",
       subtitle = "Red dots are potential anomalies; log scale",
       x = "Region", y = "Daily pedestrians per site") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), legend.position = "none")

site_volume %>%
  group_by(District) %>%
  slice_max(ped_volume, n = 10, with_ties = FALSE) %>%
  ungroup() %>%
  transmute(District, site_id, Community, OnRoad,
            avg_daily_peds = round(ped_volume), days_counted) %>%
  knitr::kable(caption = "Top 10 counted sites per region by average daily pedestrians")
Top 10 counted sites per region by average daily pedestrians
District site_id Community OnRoad avg_daily_peds days_counted
Bay 77-00006N St Clair M 29 86 1
Bay 77-00007N Marine City Brown St / M 29 61 1
Bay 77-00075N Marine City Ward St / M 29 53 1
Bay 26-00002N Beaverton Ross St 42 1
Bay 26-00003N Beaverton Ross St 37 1
Bay 32-00001N Port Austin M 25 Shoulders 11 3
Bay 74-00003N Sanilac M 25 Shoulders 9 3
Bay 79-00002N Millington M 15 Shoulders 7 3
Bay 79-00001N Columbia M 25 Shoulders 5 3
Bay 77-00067N Ira Twp - St. Clair Church Rd 4 1
Grand 61-00002N Muskegon Shared use path on Port City Blvd, Muskegon 272 3
Grand 41-00077N Walker Northridge Dr 98 1
Grand 34-00007N Ionia Dexter St 59 1
Grand 41-00080N Walker 3 Mile Rd 42 1
Grand 41-00075N Walker Fruit Ridge Ave 31 1
Grand 70-00052N Grand Haven N Beacon Blvd 26 2
Grand 03-00023N Allegan Monroe St 19 1
Grand 70-00054N Grand Haven Jackson St 19 2
Grand 70-00001N Zeeland 196 BR 17 1
Grand 08-00004N Hastings M-43 16 1
Metro 82-00022N Detroit Warren Ave 397 3
Metro 82-00131N Detroit St Antoine St 252 1
Metro 82-00128N Detroit Gratiot Ave 202 1
Metro 82-00007N Detroit US-12 (Michigan Ave) 159 3
Metro 82-00020N Detroit Livernois Ave 106 1
Metro 82-00063N Detroit US-12 (Michigan Ave) 97 3
Metro 82-00106N Detroit Gratiot Ave 94 1
Metro 63-00095N Hazel Park Dequindre St 94 2
Metro 63-00079N Hazel Park Dequindre Rd 92 2
Metro 82-00095N Detroit S Chrysler Dr 90 2
North 15-00008N Charlevoix W Crosswalk on Park Ave at US 31 3893 1
North 15-00006N Charlevoix N Crosswalk on US 31 at Bridge Park Dr 2671 2
North 15-00004N Charlevoix W Crosswalk on E Mason St at US 31 2352 2
North 15-00013N Charlevoix N Crosswalk on US 31 at E Mason St 2000 2
North 15-00005N Charlevoix Mid-Block Crossing on US 31 1901 2
North 15-00003N Charlevoix W Crosswalk on Antrim St at US 31 / N Crosswalk on US 31 at Antrim St 1638 2
North 15-00012N Charlevoix E Crosswalk on E Antrim St at US 31 1026 2
North 15-00011N Charlevoix Bridge St 1010 1
North 35-00014N East Tawas Newman St 660 2
North 45-00008N Suttons Bay Twp - Leelanau Crosswalk on Jefferson just west of M22 644 2
Southwest 11-00012N New Buffalo (WC 01) New Buffalo Welcome Center 652 91
Southwest 12-00001N Coldwater (WC 02) Coldwater Welcome Center 172 30
Southwest 11-00034N St Joseph Forres St 50 2
Southwest 11-00032N St Joseph Central Ave 36 2
Southwest 39-00038N Oshtemo Twp - Kalamazoo Maple Hill Dr 31 1
Southwest 11-00028N St Joseph M-63 16 2
Southwest 11-00002N Benton Twp - Berrien Cinema Dr. 15 3
Southwest 11-00003N Benton Twp - Berrien M 139 12 3
Southwest 11-00005N Benton Twp - Berrien E Napier Ave 11 3
Southwest 11-00009N St Joseph Twp - Berrien Lakeshore Dr 10 3
Superior 21-00010N Escanaba 12th Ave 845 3
Superior 17-00026N Sault Ste Marie Portage Ave 242 4
Superior 17-00025N Sault Ste Marie Portage Ave 162 4
Superior 17-00027N Sault Ste Marie Ashmun St 124 4
Superior 17-00019N Sault Ste Marie Ashmun St 111 7
Superior 17-00032N Sault Ste Marie Ashmun St 101 1
Superior 17-00028N Sault Ste Marie Ashmun St 89 4
Superior 17-00033N Sault Ste Marie Ashmun St 86 1
Superior 17-00017N Sault Ste Marie Ashmun St 79 7
Superior 17-00031N Sault Ste Marie Spruce St 70 1
University 23-00017N Eaton Rapids Main St 834 1
University 23-00014N Eaton Rapids Hamlin St 704 1
University 23-00016N Eaton Rapids Main St 598 1
University 23-00015N Eaton Rapids Hamlin St 483 1
University 58-00002N Petersburg (WC 03) Dundee Welcome Center 387 30
University 47-00002N Howell Grand River Ave 284 3
University 33-00021N East Lansing Crosswalk over Burcham Dr. W of Hagadorn 268 3
University 81-00003N Ann Arbor Sidewalk/Roadway of N side Bonisteel 237 44
University 33-00020N East Lansing Crosswalk over Burcham Dr. E of Hagadorn 165 3
University 81-00004N Ann Arbor Sidewalk/Roadway of S side Bonisteel 114 44

3. Count-site locations

# Location IDs may be leg-level ("82-00097N_WB") or site-level ("82-00097N"); both map to the site.
locs <- read.csv(resolve_file(CONFIG$location_file), stringsAsFactors = FALSE) %>%
  transmute(site_id = str_remove(str_trim(as.character(Location)), CONFIG$leg_suffix_regex),
            Latitude, Longitude) %>%
  clean_coords("Count-site locations") %>%
  group_by(site_id) %>%
  summarise(Latitude = mean(Latitude), Longitude = mean(Longitude), .groups = "drop")
## Count-site locations: kept 167 of 282 rows with valid Michigan coordinates.
sites <- inner_join(site_volume, locs, by = "site_id")

message(sprintf("%d of %d counted sites have coordinates.", nrow(sites), nrow(site_volume)))
## 12 of 380 counted sites have coordinates.
if (nrow(sites) == 0) {
  stop("No count sites matched the location file. Example count site IDs: ",
       paste(head(site_volume$site_id), collapse = ", "),
       " | example location IDs: ", paste(head(locs$site_id), collapse = ", "))
}

sites_sf <- sites %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE) %>%
  st_transform(CONFIG$crs_meters)

4. Pedestrian crashes

crash_raw <- read.csv(resolve_file(CONFIG$crash_file), stringsAsFactors = FALSE)
cc <- CONFIG$crash_cols

col_lat   <- find_col(crash_raw, cc$lat,  c("latitude", "(^|[._ ])lat($|[._ ])"), "lat")
col_lon   <- find_col(crash_raw, cc$lon,  c("longitude", "(^|[._ ])(lon|long|lng)($|[._ ])"), "lon")
col_ped   <- if (isTRUE(CONFIG$crashes_are_ped_only)) NULL else
  find_col(crash_raw, cc$ped_flag,
           c("pedestrian.*(involved|flag|crash)", "^crash.*pedestrian", "^pedestrian$", "pedestrian", "(^|[._ ])ped($|[._ ])"),
           "ped_flag")
col_sev   <- find_col(crash_raw, cc$severity, c("worst.*injury", "severity", "kabco", "injury"), "severity", required = FALSE)
col_year  <- find_col(crash_raw, cc$year, c("crash.*year", "^year$"), "year", required = FALSE)
col_date  <- find_col(crash_raw, cc$date, c("crash.*date", "^date$"), "date", required = FALSE)
col_id    <- find_col(crash_raw, cc$id,   c("crash.*(id|instance)", "^(crash_?)?id$", "instance"), "id", required = FALSE)
col_speed <- find_col(crash_raw, cc$speed, c("speed.*limit.*(site|crash)", "^speed.*limit"), "speed", required = FALSE)
col_int   <- find_col(crash_raw, cc$intersection, c("^crash.*intersection$", "^intersection"), "intersection", required = FALSE)
col_city  <- find_col(crash_raw, cc$city, c("city.*township", "^city", "municipal"), "city", required = FALSE)
col_reg   <- find_col(crash_raw, cc$region, c("mdot.*region", "region"), "region", required = FALSE)

# Show what was detected so it can be checked
tibble(field = c("latitude", "longitude", "pedestrian flag", "severity", "year", "date", "crash id",
                 "speed limit", "intersection", "city/township", "MDOT region"),
       column_used = c(col_lat, col_lon, col_ped %||% "(file already pedestrian-only)",
                       col_sev %||% "(none)", col_year %||% "(none)", col_date %||% "(none)",
                       col_id %||% "(row number)", col_speed %||% "(none)", col_int %||% "(none)",
                       col_city %||% "(none)", col_reg %||% "(none)")) %>%
  knitr::kable(caption = "Crash columns detected")
Crash columns detected
field column_used
latitude Crash.Latitude
longitude Crash.Longitude
pedestrian flag Crash..Pedestrian
severity (none)
year Crash.Year
date (none)
crash id Crash.Instance
speed limit Speed.Limit.at.Crash.Site
intersection Crash..Intersection
city/township City.or.Township
MDOT region MDOT.Region..2015.Present. 
col_or <- function(col, f, default) if (!is.null(col)) f(crash_raw[[col]]) else default

crashes <- tibble(
  crash_id   = col_or(col_id, as.character, as.character(seq_len(nrow(crash_raw)))),
  Latitude   = crash_raw[[col_lat]],
  Longitude  = crash_raw[[col_lon]],
  is_ped     = col_or(col_ped, is_yes, TRUE),
  severity   = col_or(col_sev, parse_kabco, "Unknown"),
  crash_date = col_or(col_date, function(x) as_date(parse_date_time(x, orders = c("Ymd", "mdY", "Ymd HMS", "mdY HMS", "mdY IMSp"), quiet = TRUE)), as.Date(NA)),
  crash_year = col_or(col_year, function(x) as.integer(to_num(x)), NA_integer_),
  speed_mph  = col_or(col_speed, to_num, NA_real_),
  is_intersection = col_or(col_int, function(x) {
    t <- str_to_lower(str_trim(as.character(x)))
    str_detect(t, "intersection|^y|^true|^1$") & !str_detect(t, "non|^no")
  }, NA),
  city       = col_or(col_city, function(x) str_trim(sub(".*:", "", x)), NA_character_),   # "Macomb County: Warren" -> "Warren"
  region     = col_or(col_reg, as.character, NA_character_)
) %>%
  mutate(crash_year = coalesce(as.integer(year(crash_date)), crash_year),
         speed_mph  = if_else(between(speed_mph, 5, 85), speed_mph, NA_real_))

message(sprintf("%d of %d crash records are pedestrian crashes.", sum(crashes$is_ped), nrow(crashes)))
## 4732 of 7918 crash records are pedestrian crashes.
ped_crashes <- crashes %>%
  filter(is_ped) %>%
  distinct(crash_id, .keep_all = TRUE) %>%
  clean_coords("Pedestrian crashes") %>%
  mutate(severity   = factor(severity, levels = names(CONFIG$severity_weights)),
         sev_weight = unname(CONFIG$severity_weights[as.character(severity)]),
         is_KA      = severity %in% c("K", "A"),
         high_speed = coalesce(speed_mph >= CONFIG$high_speed_mph, FALSE))
## Pedestrian crashes: kept 4677 of 4686 rows with valid Michigan coordinates.
if (nrow(ped_crashes) == 0) stop("No pedestrian crashes with valid coordinates. Check the pedestrian flag and coordinate columns.")

has_severity <- !is.null(col_sev) && any(ped_crashes$severity != "Unknown")
if (!has_severity) {
  message("No crash severity column found: every crash is weighted equally. ",
          "Add 'Worst Injury in Crash' (or another KABCO field) to the crash export to weight fatal/serious crashes.")
}
## No crash severity column found: every crash is weighted equally. Add 'Worst Injury in Crash' (or another KABCO field) to the crash export to weight fatal/serious crashes.
# Points to highlight on maps: fatal/serious crashes if severity exists, otherwise high-speed crashes
ped_crashes$highlight <- if (has_severity) ped_crashes$is_KA else ped_crashes$high_speed
highlight_label <- if (has_severity) "Fatal/serious ped crashes" else
  sprintf("Ped crashes on %d+ mph roads", CONFIG$high_speed_mph)

# Study period length (used for crash rates)
yrs <- ped_crashes$crash_year[!is.na(ped_crashes$crash_year)]
n_years <- if (!is.null(CONFIG$study_years)) CONFIG$study_years else
  if (length(yrs) > 0) diff(range(yrs)) + 1 else { warning("No crash dates found; assuming a 1-year study period."); 1 }
message(sprintf("Study period: %s year(s)%s", n_years,
                if (length(yrs)) sprintf(" (%d-%d)", min(yrs), max(yrs)) else ""))
## Study period: 5 year(s) (2020-2024)
ped_crashes %>%
  summarise(ped_crashes        = n(),
            intersection_pct   = round(100 * mean(is_intersection, na.rm = TRUE), 1),
            median_speed_limit = median(speed_mph, na.rm = TRUE),
            high_speed_pct     = round(100 * mean(high_speed), 1))
## # A tibble: 1 × 4
##   ped_crashes intersection_pct median_speed_limit high_speed_pct
##         <int>            <dbl>              <dbl>          <dbl>
## 1        4677               44                 35           30.6
if (has_severity) count(ped_crashes, severity, name = "ped_crashes")

crash_sf <- ped_crashes %>%
  st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326, remove = FALSE) %>%
  st_transform(CONFIG$crs_meters)
ped_crashes %>%
  filter(!is.na(speed_mph)) %>%
  mutate(location = case_when(is_intersection %in% TRUE  ~ "Intersection",
                              is_intersection %in% FALSE ~ "Non-intersection",
                              TRUE ~ "Unknown")) %>%
  ggplot(aes(factor(speed_mph), fill = location)) +
  geom_bar() +
  scale_fill_manual(values = c("Intersection" = "#2b8cbe", "Non-intersection" = "#fdae61", "Unknown" = "grey70")) +
  labs(title = "Pedestrian crashes by speed limit", x = "Speed limit at crash site (mph)",
       y = "Pedestrian crashes", fill = NULL) +
  theme_minimal()

5. Assign crashes to count sites

Count sites outside the area the crash file covers are removed first (otherwise they would look artificially safe). Each crash then goes to its nearest count site, and only if that site is within 76 m, so one crash is never counted at several nearby sites.

n_sites_all <- nrow(sites_sf)

if (CONFIG$limit_sites_to_crash_area) {
  crash_regions <- unique(na.omit(crash_sf$region))
  if (length(crash_regions) > 0) {
    sites_sf <- filter(sites_sf, str_to_lower(District) %in% str_to_lower(crash_regions))
  }
  coverage <- crash_sf %>% st_union() %>% st_convex_hull() %>% st_buffer(1000)
  sites_sf <- sites_sf[lengths(st_intersects(sites_sf, coverage)) > 0, ]
  message(sprintf("%d of %d located count sites are inside the crash data coverage%s.",
                  nrow(sites_sf), n_sites_all,
                  if (length(crash_regions)) paste0(" (region: ", paste(crash_regions, collapse = ", "), ")") else ""))
}
## 12 of 12 located count sites are inside the crash data coverage (region: Metro).
if (nrow(sites_sf) == 0) stop("No count sites fall inside the area covered by the crash file.")

nearest <- st_nearest_feature(crash_sf, sites_sf)
crash_sf <- crash_sf %>%
  mutate(nearest_site   = sites_sf$site_id[nearest],
         dist_to_site_m = as.numeric(st_distance(crash_sf, sites_sf[nearest, ], by_element = TRUE)),
         assigned_site  = if_else(dist_to_site_m <= CONFIG$match_dist_m, nearest_site, NA_character_))

site_crashes <- crash_sf %>%
  st_drop_geometry() %>%
  filter(!is.na(assigned_site)) %>%
  group_by(assigned_site) %>%
  summarise(ped_crashes          = n(),
            ka_crashes           = sum(is_KA),
            fatal_crashes        = sum(severity == "K"),
            severity_score       = sum(sev_weight),
            intersection_crashes = sum(is_intersection, na.rm = TRUE),
            high_speed_crashes   = sum(high_speed),
            median_speed_mph     = median(speed_mph, na.rm = TRUE),
            .groups = "drop")

sites_sf <- sites_sf %>%
  left_join(site_crashes, by = c("site_id" = "assigned_site")) %>%
  mutate(across(c(ped_crashes, ka_crashes, fatal_crashes, severity_score,
                  intersection_crashes, high_speed_crashes), ~ coalesce(as.numeric(.x), 0)))

tibble(ped_crashes_total     = nrow(crash_sf),
       matched_to_count_site = sum(!is.na(crash_sf$assigned_site)),
       pct_matched           = round(100 * mean(!is.na(crash_sf$assigned_site)), 1),
       sites_in_model        = nrow(sites_sf),
       sites_with_crashes    = sum(sites_sf$ped_crashes > 0))
## # A tibble: 1 × 5
##   ped_crashes_total matched_to_count_site pct_matched sites_in_model
##               <int>                 <int>       <dbl>          <int>
## 1              4677                     3         0.1             12
## # ℹ 1 more variable: sites_with_crashes <int>

6. Site screening model (SPF + Empirical Bayes)

model_df <- sites_sf %>%
  st_drop_geometry() %>%
  mutate(log_vol = log(ped_volume + 1))

if (sum(model_df$ped_crashes) == 0) stop("No crashes matched any count site; try a larger CONFIG$match_dist_m.")
if (sum(model_df$ped_crashes) < 20 || nrow(model_df) < 30)
  warning("Few sites or crashes: SPF estimates will be unstable. Treat the ranking as indicative.")

fit_spf <- function(df) {
  nb <- tryCatch(suppressWarnings(MASS::glm.nb(ped_crashes ~ log_vol, data = df)),
                 error = function(e) NULL)
  if (!is.null(nb) && is.finite(nb$theta) && nb$theta < 1e3) {
    return(list(fit = nb, theta = nb$theta, type = "Negative binomial"))
  }
  # Fallback: Poisson mean model + method-of-moments overdispersion
  pois <- glm(ped_crashes ~ log_vol, family = poisson, data = df)
  mu   <- fitted(pois)
  k    <- max(sum((df$ped_crashes - mu)^2 - mu) / sum(mu^2), 0.01)
  list(fit = pois, theta = 1 / k, type = "Poisson + moment-based overdispersion")
}

spf  <- fit_spf(model_df)
beta <- unname(coef(spf$fit)["log_vol"])

model_df <- model_df %>%
  mutate(
    predicted        = as.numeric(predict(spf$fit, newdata = model_df, type = "response")),
    eb_weight        = 1 / (1 + predicted / spf$theta),          # HSM: w = 1 / (1 + k * mu), k = 1/theta
    eb_expected      = eb_weight * predicted + (1 - eb_weight) * ped_crashes,
    excess_crashes   = eb_expected - predicted,                  # potential for safety improvement
    crash_rate_per_M = if_else(ped_volume > 0,
                               ped_crashes / (ped_volume * 365 * n_years) * 1e6, NA_real_),
    pct_high_speed   = if_else(ped_crashes > 0, 100 * high_speed_crashes / ped_crashes, NA_real_),
    excess_pct       = percent_rank(excess_crashes),
    concern_tier     = case_when(
      excess_crashes > 0 & excess_pct >= 1 - CONFIG$high_tier_share ~ "High",
      excess_crashes > 0                                            ~ "Elevated",
      TRUE                                                          ~ "Expected or below"),
    concern_tier     = factor(concern_tier, levels = c("High", "Elevated", "Expected or below"))
  ) %>%
  arrange(desc(excess_crashes), desc(ka_crashes)) %>%
  mutate(concern_rank = row_number())

tibble(model = spf$type,
       sites = nrow(model_df),
       volume_exponent = round(beta, 3),
       overdispersion_theta = round(spf$theta, 2),
       AIC = round(AIC(spf$fit), 1)) %>%
  knitr::kable(caption = "Safety performance function")
Safety performance function
model sites volume_exponent overdispersion_theta AIC
Negative binomial 12 0.436 0.35 20.6 
count(model_df, concern_tier, name = "sites")
## # A tibble: 2 × 2
##   concern_tier      sites
##   <fct>             <int>
## 1 High                  2
## 2 Expected or below    10

The volume exponent is 0.44. Because it is below 1, crashes grow more slowly than foot traffic: per-pedestrian risk falls at busier sites (the ‘safety in numbers’ effect), so raw crash counts alone would over-flag busy locations.

tier_colors <- c("High" = "#b2182b", "Elevated" = "#ef8a62", "Expected or below" = "#67a9cf")

curve_df <- tibble(log_vol = seq(min(model_df$log_vol), max(model_df$log_vol), length.out = 200)) %>%
  mutate(ped_volume = exp(log_vol) - 1)
curve_df$predicted <- as.numeric(predict(spf$fit, newdata = curve_df, type = "response"))

ggplot(model_df, aes(ped_volume + 1, ped_crashes)) +
  geom_jitter(aes(color = concern_tier), width = 0, height = 0.15, alpha = 0.8, size = 2) +
  geom_line(data = curve_df, aes(y = predicted), linewidth = 1) +
  scale_x_log10(labels = scales::comma) +
  scale_color_manual(values = tier_colors) +
  labs(title = "Pedestrian crashes vs. pedestrian volume",
       subtitle = sprintf("Line = SPF prediction over %s year(s)", n_years),
       x = "Average daily pedestrians (log scale)", y = "Pedestrian crashes at site", color = "Concern tier") +
  theme_minimal()

drop_sev <- function(df) if (has_severity) df else select(df, -any_of(c("ka_crashes", "fatal_crashes", "severity_score")))

top_sites <- model_df %>% slice_head(n = 20)

top_sites %>%
  mutate(label = sprintf("%s - %s (%s)", Community, str_trunc(OnRoad, 30), site_id),
         label = factor(label, levels = rev(label))) %>%
  ggplot(aes(excess_crashes, label, fill = concern_tier)) +
  geom_col() +
  scale_fill_manual(values = tier_colors) +
  labs(title = "Top 20 count sites by excess expected pedestrian crashes",
       x = "Excess crashes (EB expected minus SPF prediction)", y = NULL, fill = "Concern tier") +
  theme_minimal()

top_sites %>%
  transmute(rank = concern_rank, site_id, Community, OnRoad,
            avg_daily_peds = round(ped_volume), days_counted, ped_crashes, ka_crashes,
            predicted = round(predicted, 2), eb_expected = round(eb_expected, 2),
            excess = round(excess_crashes, 2), rate_per_M_peds = round(crash_rate_per_M, 2),
            pct_high_speed = round(pct_high_speed), concern_tier) %>%
  drop_sev() %>%
  knitr::kable(caption = "Highest-concern count sites")
Highest-concern count sites
rank site_id Community OnRoad avg_daily_peds days_counted ped_crashes predicted eb_expected excess rate_per_M_peds pct_high_speed concern_tier
1 82-00058N Detroit Michigan Ave 12 3 2 0.29 1.06 0.77 88.86 0 High
2 63-00010N Pontiac M 59 btw Thorpe St & Lincoln St 6 3 1 0.22 0.52 0.30 86.52 100 High
3 82-00059N Detroit Michigan Ave 0 3 0 0.11 0.08 -0.02 0.00 NA Expected or below
4 63-00012N Pontiac M 59 btw Franklin Blvd & State St 1 3 0 0.12 0.09 -0.03 0.00 NA Expected or below
5 82-00065N Detroit US-12 (Michigan Ave) 2 3 0 0.14 0.10 -0.04 0.00 NA Expected or below
6 82-00060N Detroit Michigan Ave 7 3 0 0.23 0.14 -0.09 0.00 NA Expected or below
7 63-00009N Pontiac M 59 btw Dakota St & Waldo St 8 3 0 0.24 0.14 -0.10 0.00 NA Expected or below
8 63-00005N Orion Twp - Oakland Mid-Block Crossing on E Greenshield Rd 8 3 0 0.25 0.15 -0.10 0.00 NA Expected or below
9 63-00011N Pontiac M 59 btw Henderson St & Palmer St 13 3 0 0.30 0.16 -0.14 0.00 NA Expected or below
10 82-00057N Detroit Michigan Ave 16 3 0 0.32 0.17 -0.15 0.00 NA Expected or below
11 82-00004N Dearborn Heights Crosswalk on Ann Arbor Trail at Parkland 18 3 0 0.34 0.17 -0.16 0.00 NA Expected or below
12 63-00013N Novi Unmarked crossing on 11 Mile 61 3 0 0.56 0.22 -0.35 0.00 NA Expected or below

7. Crash hot spots (Getis-Ord Gi*)

cell <- CONFIG$hex_size_m

grid <- st_make_grid(crash_sf, cellsize = cell, square = FALSE)
grid <- st_sf(cell_id = seq_along(grid), geometry = grid)

# Study area for the test = cells within hex_context_m of at least one crash
keep <- lengths(st_is_within_distance(st_centroid(st_geometry(grid)), crash_sf,
                                      dist = CONFIG$hex_context_m)) > 0
grid <- grid[keep, ]

cell_hit <- st_intersects(crash_sf, grid)
crash_sf$cell_id <- vapply(cell_hit, function(i) if (length(i)) grid$cell_id[i[1]] else NA_integer_, integer(1))

cell_stats <- crash_sf %>%
  st_drop_geometry() %>%
  filter(!is.na(cell_id)) %>%
  group_by(cell_id) %>%
  summarise(crashes = n(), ka_crashes = sum(is_KA), severity_score = sum(sev_weight), .groups = "drop")

grid <- grid %>%
  left_join(cell_stats, by = "cell_id") %>%
  mutate(across(c(crashes, ka_crashes, severity_score), ~ coalesce(as.numeric(.x), 0)))

# Neighbors = the cell itself plus the 6 touching hexagons
nb <- st_is_within_distance(st_centroid(st_geometry(grid)), dist = cell * 1.5)

grid <- grid %>%
  mutate(gi_z    = gi_star(severity_score, nb),
         p_fdr   = p.adjust(pnorm(gi_z, lower.tail = FALSE), method = "BH"),
         hotspot = case_when(gi_z >= 2.576 ~ "Hot spot (99%)",
                             gi_z >= 1.960 ~ "Hot spot (95%)",
                             TRUE          ~ "Not significant"))

hot_cells <- filter(grid, gi_z >= 1.96)
count(st_drop_geometry(grid), hotspot, name = "cells")
##           hotspot cells
## 1  Hot spot (95%)   616
## 2  Hot spot (99%)   888
## 3 Not significant 19631
# Merge touching hot cells into zones and summarize each zone
if (nrow(hot_cells) > 0) {
  zones <- hot_cells %>%
    st_buffer(1) %>%           # closes hairline gaps between hexagons
    st_union() %>%
    st_cast("POLYGON") %>%
    st_sf(geometry = .) %>%
    mutate(zone_id = row_number())

  crash_zone <- st_intersects(crash_sf, zones)
  crash_sf$zone_id <- vapply(crash_zone, function(i) if (length(i)) zones$zone_id[i[1]] else NA_integer_, integer(1))

  zone_stats <- crash_sf %>%
    st_drop_geometry() %>%
    filter(!is.na(zone_id)) %>%
    group_by(zone_id) %>%
    summarise(ped_crashes      = n(),
              ka_crashes       = sum(is_KA),
              fatal_crashes    = sum(severity == "K"),
              severity_score   = sum(sev_weight),
              pct_intersection = round(100 * mean(is_intersection, na.rm = TRUE)),
              pct_high_speed   = round(100 * mean(high_speed)),
              median_speed_mph = median(speed_mph, na.rm = TRUE),
              city             = mode_of(city),
              .groups = "drop")

  zone_peak <- st_join(st_centroid(hot_cells), zones) %>%
    st_drop_geometry() %>%
    group_by(zone_id) %>%
    summarise(cells = n(), peak_gi_z = max(gi_z), .groups = "drop")

  zone_centers <- st_centroid(st_geometry(zones))
  near_site    <- st_nearest_feature(zone_centers, sites_sf)
  zone_ll      <- st_coordinates(st_transform(zone_centers, 4326))

  zones <- zones %>%
    left_join(zone_stats, by = "zone_id") %>%
    left_join(zone_peak,  by = "zone_id") %>%
    mutate(area_km2       = as.numeric(st_area(geometry)) / 1e6,
           center_lat     = zone_ll[, "Y"],
           center_lon     = zone_ll[, "X"],
           nearest_site   = sites_sf$site_id[near_site],
           nearest_road   = sites_sf$OnRoad[near_site],
           nearest_site_m = round(as.numeric(st_distance(zone_centers, sites_sf[near_site, ], by_element = TRUE)))) %>%
    mutate(across(c(ped_crashes, ka_crashes, fatal_crashes, severity_score), ~ coalesce(.x, 0))) %>%
    arrange(desc(severity_score), desc(ka_crashes), desc(pct_high_speed)) %>%
    mutate(zone_rank = row_number())

  # ---- Readable summary table ---------------------------------------------
  fmt_num  <- function(x, digits = 0) ifelse(is.na(x), "n/a", formatC(x, format = "f", digits = digits, big.mark = ","))
  fmt_pct  <- function(x) ifelse(is.na(x), "n/a", paste0(round(x), "%"))
  fmt_dist <- function(m) ifelse(m * 3.28084 < 1000,
                                 paste(fmt_num(m * 3.28084), "ft"),
                                 paste(fmt_num(m / 1609.34, 1), "mi"))

  zone_table <- zones %>%
    st_drop_geometry() %>%
    slice_head(n = 15) %>%
    transmute(
      zone_rank,
      city             = coalesce(city, "n/a"),
      ped_crashes      = fmt_num(ped_crashes),
      ka_crashes       = fmt_num(ka_crashes),
      fatal_crashes    = fmt_num(fatal_crashes),
      severity_score   = fmt_num(severity_score),
      pct_intersection = fmt_pct(pct_intersection),
      pct_high_speed   = fmt_pct(pct_high_speed),
      median_speed_mph = fmt_num(median_speed_mph),
      area_sq_mi       = fmt_num(area_km2 * 0.386102, 2),
      peak_gi_z        = fmt_num(peak_gi_z, 1),
      location         = sprintf("[%.4f, %.4f](https://www.google.com/maps/search/?api=1&query=%.6f,%.6f)",
                                 center_lat, center_lon, center_lat, center_lon),
      nearest_site     = sprintf("%s%s, %s away", nearest_site,
                                 ifelse(is.na(nearest_road), "", paste0(" (", str_trunc(nearest_road, 25), ")")),
                                 fmt_dist(nearest_site_m))
    ) %>%
    drop_sev()   # removes the severity columns when the crash file has no injury field

  zone_labels <- c(
    zone_rank        = "Rank",
    city             = "City",
    ped_crashes      = "Ped crashes",
    ka_crashes       = "Fatal + serious (K+A)",
    fatal_crashes    = "Fatal (K)",
    severity_score   = "Severity score",
    pct_intersection = "At intersections",
    pct_high_speed   = sprintf("On %d+ mph roads", CONFIG$high_speed_mph),
    median_speed_mph = "Median speed limit (mph)",
    area_sq_mi       = "Zone size (sq mi)",
    peak_gi_z        = "Hot spot strength (Gi* z)",
    location         = "Zone center (opens map)",
    nearest_site     = "Nearest count site"
  )

  knitr::kable(
    zone_table,
    col.names = unname(zone_labels[names(zone_table)]),
    align     = ifelse(names(zone_table) %in% c("city", "location", "nearest_site"), "l", "r"),
    caption   = sprintf(
      "Top %d pedestrian crash hot spot zones, ranked by %s. Hot spot strength is the peak Gi* z-score in the zone (1.96+ = 95%% confidence, 2.58+ = 99%%).",
      nrow(zone_table),
      if (has_severity) "severity-weighted crashes" else "number of pedestrian crashes")
  )
} else {
  zones <- NULL
  message("No statistically significant hot spot cells at this grid size. Try a larger CONFIG$hex_size_m.")
}
Top 15 pedestrian crash hot spot zones, ranked by number of pedestrian crashes. Hot spot strength is the peak Gi* z-score in the zone (1.96+ = 95% confidence, 2.58+ = 99%).
Rank City Ped crashes At intersections On 40+ mph roads Median speed limit (mph) Zone size (sq mi) Hot spot strength (Gi* z) Zone center (opens map) Nearest count site
1 Detroit 325 66% 2% 30 8.10 36.2 42.3537, -83.0600 82-00059N (Michigan Ave), 1.5 mi away
2 Detroit 283 48% 6% 35 9.71 11.4 42.4154, -83.2006 82-00004N (Crosswalk on Ann Arbor…), 6.1 mi away
3 Detroit 207 47% 9% 35 6.22 12.4 42.4275, -82.9698 82-00059N (Michigan Ave), 7.9 mi away
4 Detroit 166 45% 2% 30 6.33 9.1 42.3580, -83.2068 82-00004N (Crosswalk on Ann Arbor…), 2.9 mi away
5 Hamtramck 129 30% 5% 25 4.18 11.4 42.4175, -83.0675 82-00058N (Michigan Ave), 5.9 mi away
6 Detroit 93 41% 6% 30 3.22 9.6 42.3852, -83.1377 82-00065N (US-12 (Michigan Ave)), 4.7 mi away
7 Lincoln Park 55 78% 45% 35 2.20 6.3 42.2521, -83.1835 82-00004N (Crosswalk on Ann Arbor…), 6.8 mi away
8 Dearborn 47 55% 21% 35 2.42 5.8 42.3301, -83.1792 82-00004N (Crosswalk on Ann Arbor…), 3.8 mi away
9 Royal Oak 47 40% 6% 25 0.97 17.0 42.4840, -83.1447 82-00065N (US-12 (Michigan Ave)), 11.0 mi away
10 Detroit 39 44% 8% 35 1.56 7.2 42.3318, -83.1294 82-00065N (US-12 (Michigan Ave)), 2.4 mi away
11 Detroit 37 30% 22% 35 1.18 8.2 42.4401, -83.0204 82-00059N (Michigan Ave), 7.7 mi away
12 Detroit 35 34% 9% 35 1.02 10.0 42.3867, -83.0090 82-00059N (Michigan Ave), 4.5 mi away
13 Warren 32 12% 3% 35 1.13 6.8 42.4716, -83.0275 82-00059N (Michigan Ave), 9.8 mi away
14 Birmingham 30 40% 23% 25 0.64 12.4 42.5462, -83.2142 63-00012N (M 59 btw Franklin Blvd…), 7.7 mi away
15 Detroit 27 59% 4% 25 0.86 7.7 42.3121, -83.1220 82-00065N (US-12 (Michigan Ave)), 2.5 mi away

8. Maps

8.1 Static overview

p <- ggplot()
if (!is.null(zones)) p <- p + geom_sf(data = zones, fill = "#fdae61", color = "#d7301f", alpha = 0.6)
if (any(crash_sf$highlight)) p <- p + geom_sf(data = filter(crash_sf, highlight), color = "black", size = 0.5, alpha = 0.4)
p +
  geom_sf(data = left_join(select(sites_sf, site_id),
                           select(model_df, site_id, concern_tier), by = "site_id"),
          aes(color = concern_tier), size = 2.5) +
  scale_color_manual(values = tier_colors) +
  labs(title = "Pedestrian areas of concern",
       subtitle = sprintf("Orange = crash hot spot zones; black = %s; points = count sites", str_to_lower(highlight_label)),
       color = "Count-site tier") +
  theme_minimal()

8.2 Interactive map

library(leaflet)
to_ll <- function(x) st_transform(x, 4326)

map_sites <- sites_sf %>%
  select(site_id) %>%
  left_join(model_df, by = "site_id") %>%
  to_ll() %>%
  mutate(popup = sprintf(
    paste0("<b>%s</b><br>%s<br>Site %s (%s), legs: %s<br>",
           "Concern rank: %d of %d (<b>%s</b>)<br>Avg daily pedestrians: %s (%d day(s) counted)<br>",
           "Ped crashes within %d m: %d%s<br>Expected (SPF): %.2f | EB: %.2f | Excess: %.2f"),
    htmltools::htmlEscape(Community), htmltools::htmlEscape(str_trunc(coalesce(Description, OnRoad, ""), 90)),
    site_id, District, htmltools::htmlEscape(leg_ids),
    concern_rank, nrow(model_df), concern_tier,
    format(round(ped_volume), big.mark = ","), as.integer(days_counted),
    CONFIG$match_dist_m, as.integer(ped_crashes),
    if (has_severity) sprintf(" (fatal/serious: %d)", as.integer(ka_crashes)) else
      ifelse(ped_crashes > 0, sprintf(" (%d%% on %d+ mph roads)", round(pct_high_speed), CONFIG$high_speed_mph), ""),
    predicted, eb_expected, excess_crashes))

tier_pal <- colorFactor(unname(tier_colors), levels = names(tier_colors))

carto_key <- Sys.getenv("CARTO_API_KEY")

m <- leaflet()
m <- if (nzchar(carto_key)) {
  m %>% addTiles(
    urlTemplate = paste0("https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png?key=", carto_key),
    attribution = '&copy; <a href="https://www.openstreetmap.org/copyright">OpenStreetMap</a> contributors &copy; <a href="https://carto.com/attributions">CARTO</a>',
    options     = tileOptions(subdomains = "abcd", maxZoom = 20))
} else {
  m %>% addProviderTiles(providers$Esri.WorldGrayCanvas,
                         options = providerTileOptions(maxNativeZoom = 16, maxZoom = 19))
}
if (nrow(hot_cells) > 0) {
  hot_pal <- colorNumeric("YlOrRd", domain = range(hot_cells$gi_z))
  m <- m %>%
    addPolygons(data = to_ll(hot_cells), fillColor = ~ hot_pal(gi_z), fillOpacity = 0.55,
                color = "#7f0000", weight = 0.4, group = "Crash hot spots",
                label = ~ sprintf("Gi* z = %.1f | %d ped crashes", gi_z, as.integer(crashes))) %>%
    addLegend("bottomleft", pal = hot_pal, values = hot_cells$gi_z,
              title = "Hot spot Gi* z", group = "Crash hot spots")
}
if (any(!crash_sf$highlight)) {
  m <- m %>% addCircleMarkers(data = to_ll(filter(crash_sf, !highlight)), radius = 2, stroke = FALSE,
                              fillColor = "#636363", fillOpacity = 0.5, group = "Other ped crashes")
}
if (any(crash_sf$highlight)) {
  m <- m %>% addCircleMarkers(data = to_ll(filter(crash_sf, highlight)), radius = 3, stroke = FALSE,
                              fillColor = "black", fillOpacity = 0.8, group = highlight_label,
                              label = ~ sprintf("%s mph | %s", coalesce(as.character(speed_mph), "?"), coalesce(city, "")))
}

m %>%
  addCircleMarkers(data = map_sites, radius = ~ scales::rescale(sqrt(ped_volume), to = c(5, 13)),
                   color = "white", weight = 1, fillColor = ~ tier_pal(concern_tier), fillOpacity = 0.9,
                   popup = ~ popup, group = "Count sites (concern tier)") %>%
  addLegend("bottomright", pal = tier_pal, values = map_sites$concern_tier, title = "Count-site concern") %>%
  addLayersControl(overlayGroups = c("Crash hot spots", "Count sites (concern tier)",
                                     highlight_label, "Other ped crashes"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup("Other ped crashes")

9. Export results

site_out <- model_df %>%
  transmute(concern_rank, concern_tier, site_id, leg_ids, Community, CountyName, District, OnRoad, Description,
            Latitude, Longitude, avg_daily_peds = round(ped_volume, 1), days_counted,
            ped_crashes, ka_crashes, fatal_crashes, severity_score, intersection_crashes,
            high_speed_crashes, median_speed_mph,
            predicted = round(predicted, 3), eb_expected = round(eb_expected, 3),
            excess_crashes = round(excess_crashes, 3), crash_rate_per_M = round(crash_rate_per_M, 3)) %>%
  drop_sev()
write.csv(site_out, file.path(CONFIG$out_dir, "ped_site_concern_ranking.csv"), row.names = FALSE)

crash_sf %>%
  st_drop_geometry() %>%
  select(crash_id, crash_year, severity, speed_mph, is_intersection, city, Latitude, Longitude,
         assigned_site, dist_to_site_m, any_of("zone_id")) %>%
  write.csv(file.path(CONFIG$out_dir, "ped_crash_assignments.csv"), row.names = FALSE)

if (!is.null(zones)) {
  write.csv(st_drop_geometry(zones), file.path(CONFIG$out_dir, "ped_hotspot_zones.csv"), row.names = FALSE)
  st_write(st_transform(zones, 4326), file.path(CONFIG$out_dir, "ped_hotspot_zones.geojson"),
           delete_dsn = TRUE, quiet = TRUE)
  st_write(st_transform(hot_cells, 4326), file.path(CONFIG$out_dir, "ped_hotspot_cells.geojson"),
           delete_dsn = TRUE, quiet = TRUE)
}

list.files(CONFIG$out_dir)
## [1] "ped_crash_assignments.csv"    "ped_hotspot_cells.geojson"   
## [3] "ped_hotspot_zones.csv"        "ped_hotspot_zones.geojson"   
## [5] "ped_site_concern_ranking.csv"

Summary

  • 12 count sites inside the crash data coverage were screened against 4677 pedestrian crashes over 5 year(s); 3 crashes fell within 76 m of a count site.
  • 2 sites are High concern and 0 are Elevated (more crashes than expected for their pedestrian volume).
  • 127 crash hot spot zones were found; these include areas without counts.
  • The crash file has no severity field, so all crashes count equally. Re-export with ‘Worst Injury in Crash’ to weight fatal and serious crashes.

Limitations. Pedestrian volume is the only exposure variable; adding vehicle AADT, number of lanes, or land use to the SPF would sharpen the expected values. Most sites were counted on a single day, so volumes are noisy, especially without seasonal factors. Crash years and count years are assumed comparable.