1 Study Overview

This reproducible R Markdown workflow analyzes public feedback and reported incidents involving autonomous vehicles in Austin. The analysis is organized into four stages:

  1. descriptive assessment of incident patterns;
  2. structural topic modeling and word co-occurrence analysis;
  3. Policy Analysis 1, which develops an interpretable civic-risk screening and prioritization framework; and
  4. Policy Analysis 2, which evaluates advanced hybrid prediction, calibration, conformal triage, policy simulation, and spatial risk patterns.

2 1. Setup and Libraries

All required packages are loaded once at the beginning of the document. Duplicate package calls from the original script have been consolidated.

knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE,
  alert = FALSE,
  fig.align = "center",
  fig.width = 14,
  fig.height = 6,
  dpi = 300
)

set.seed(1234)

set.seed(1234)

required_packages <- c(
  "httr",
  "jsonlite",
  "dplyr",
  "purrr",
  "writexl",
  "readxl",
  "tidyverse",
  "quanteda",
  "stm",
  "tidytext",
  "textstem",
  "ggplot2",
  "ggrepel",
  "lubridate",
  "quanteda.textplots",
  "quanteda.textstats",
  "igraph",
  "ggraph",
  "widyr",
  "wordcloud2",
  "RColorBrewer",
  "text2vec",
  "uwot",
  "glmnet",
  "caret",
  "pROC",
  "tidymodels",
  "PRROC",
  "Matrix",
  "yardstick",
  "patchwork",
  "scales",
  "forcats",
  "sf",
  "units"
)

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0) {
  stop(
    "Install the following packages before knitting: ",
    paste(missing_packages, collapse = ", ")
  )
}

invisible(
  lapply(required_packages, library, character.only = TRUE)
)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attaching package: 'purrr'
## The following object is masked from 'package:jsonlite':
## 
##     flatten
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ readr     2.2.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()  masks stats::filter()
## ✖ purrr::flatten() masks jsonlite::flatten()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Package version: 4.4
## Unicode version: 15.1
## ICU version: 74.1
## 
## Parallel computing: 56 of 56 threads used.
## 
## See https://quanteda.io for tutorials and examples.
## 
## stm v1.3.8 successfully loaded. See ?stm for help. 
##  Papers, resources, and other materials at structuraltopicmodel.com
## 
## Loading required package: koRpus.lang.en
## 
## Loading required package: koRpus
## 
## Loading required package: sylly
## 
## For information on available language packages for 'koRpus', run
## 
##   available.koRpus.lang()
## 
## and see ?install.koRpus.lang()
## 
## 
## 
## Attaching package: 'koRpus'
## 
## 
## The following objects are masked from 'package:quanteda':
## 
##     tokens, types
## 
## 
## The following object is masked from 'package:readr':
## 
##     tokenize
## 
## 
## 
## Attaching package: 'igraph'
## 
## 
## The following object is masked from 'package:quanteda.textplots':
## 
##     as.igraph
## 
## 
## The following objects are masked from 'package:lubridate':
## 
##     %--%, union
## 
## 
## The following object is masked from 'package:tidyr':
## 
##     crossing
## 
## 
## The following object is masked from 'package:tibble':
## 
##     as_data_frame
## 
## 
## The following objects are masked from 'package:purrr':
## 
##     compose, simplify
## 
## 
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## 
## 
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## 
## 
## The following object is masked from 'package:base':
## 
##     union
## 
## 
## 
## Attaching package: 'text2vec'
## 
## 
## The following object is masked from 'package:igraph':
## 
##     normalize
## 
## 
## Loading required package: Matrix
## 
## 
## Attaching package: 'Matrix'
## 
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## Loaded glmnet 5.0
## 
## Loading required package: lattice
## 
## 
## Attaching package: 'lattice'
## 
## 
## The following object is masked from 'package:stm':
## 
##     cloud
## 
## 
## 
## Attaching package: 'caret'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     lift
## 
## 
## The following object is masked from 'package:httr':
## 
##     progress
## 
## 
## Type 'citation("pROC")' for a citation.
## 
## 
## Attaching package: 'pROC'
## 
## 
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var
## 
## 
## ── Attaching packages ────────────────────────────────────── tidymodels 1.5.0 ──
## 
## ✔ broom        1.0.12     ✔ rsample      1.3.2 
## ✔ dials        1.4.4      ✔ tailor       0.1.0 
## ✔ infer        1.1.0      ✔ tune         2.1.0 
## ✔ modeldata    1.5.1      ✔ workflows    1.3.0 
## ✔ parsnip      1.6.0      ✔ workflowsets 1.1.1 
## ✔ recipes      1.3.3      ✔ yardstick    1.4.0 
## 
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ igraph::as_data_frame()  masks tibble::as_data_frame(), dplyr::as_data_frame()
## ✖ rsample::calibration()   masks caret::calibration()
## ✖ igraph::compose()        masks purrr::compose()
## ✖ igraph::crossing()       masks tidyr::crossing()
## ✖ dials::degree()          masks igraph::degree()
## ✖ scales::discard()        masks purrr::discard()
## ✖ Matrix::expand()         masks tidyr::expand()
## ✖ dplyr::filter()          masks stats::filter()
## ✖ workflows::fit()         masks tailor::fit(), parsnip::fit(), infer::fit(), text2vec::fit()
## ✖ recipes::fixed()         masks stringr::fixed()
## ✖ purrr::flatten()         masks jsonlite::flatten()
## ✖ dplyr::lag()             masks stats::lag()
## ✖ caret::lift()            masks purrr::lift()
## ✖ dials::neighbors()       masks igraph::neighbors()
## ✖ Matrix::pack()           masks tidyr::pack()
## ✖ yardstick::precision()   masks caret::precision()
## ✖ yardstick::recall()      masks caret::recall()
## ✖ yardstick::sensitivity() masks caret::sensitivity()
## ✖ igraph::simplify()       masks purrr::simplify()
## ✖ yardstick::spec()        masks readr::spec()
## ✖ yardstick::specificity() masks caret::specificity()
## ✖ recipes::step()          masks stats::step()
## ✖ Matrix::unpack()         masks tidyr::unpack()
## ✖ recipes::update()        masks Matrix::update(), stats::update()
## 
## Loading required package: rlang
## 
## 
## Attaching package: 'rlang'
## 
## 
## The following object is masked from 'package:igraph':
## 
##     is_named
## 
## 
## The following objects are masked from 'package:purrr':
## 
##     flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
##     flatten_raw, invoke, splice
## 
## 
## The following objects are masked from 'package:jsonlite':
## 
##     flatten, unbox
## 
## 
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
## 
## udunits database from C:/Users/mvx13/AppData/Local/R/win-library/4.6/units/share/udunits/udunits2.xml
data_path <- file.path(params$data_directory, params$data_file)

if (!file.exists(data_path)) {
  stop("Input file not found: ", data_path)
}

3 2. Data Import and Preparation

The source spreadsheet is imported, school-bus stop-arm violations are excluded, and the retained incident records are used throughout the subsequent analyses.

### READ DATA


dat_geocoded <- readxl::read_excel(data_path)
dim(dat_geocoded)
## [1] 295  15
dat_geocoded= subset(dat_geocoded, IncidentType!="School bus stop arm violation")
dat= dat_geocoded
table(dat$IncidentType, useNA = "ifany")
## 
##      Blocking Traffic             Collision  Ignore APD Direction 
##                    46                    13                    26 
##             Near Miss              Nuisance        Safety Concern 
##                    37                    43                    97 
## Safety Issue - School 
##                     5
dim(dat)
## [1] 267  15

4 3. Descriptive Analysis

4.1 3.1 Distribution of Incident Types

incident_plot <- dat %>%
  count(IncidentType, name = "n") %>%
  mutate(
    IncidentType = fct_reorder(IncidentType, n),
    percent = 100 * n / sum(n)
  )

ggplot(incident_plot, aes(x = IncidentType, y = n)) +
  geom_col(width = 0.72, fill = "#798E87", alpha = 0.92) +
  geom_text(
    aes(label = paste0(n, " (", round(percent, 1), "%)")),
    hjust = -0.10,
    size = 4.5
  ) +
  coord_flip() +
  scale_y_continuous(expand = expansion(mult = c(0, 0.18))) +
  labs(
    x = NULL,
    y = "Number of feedback records"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.y = element_text(size = 13),
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14, face = "bold"),
    plot.margin = margin(10, 25, 10, 10)
  )

4.2 3.2 Company-by-Year Distribution

company_year_dot <- dat %>%
  count(Company, Year, name = "n") %>%
  complete(
    Company,
    Year = 2023:2026,
    fill = list(n = 0)
  ) %>%
  mutate(
    Year = factor(Year),
    Company = fct_reorder(Company, n, .fun = sum)
  )

ggplot(company_year_dot, aes(x = Year, y = Company)) +
  geom_point(
    aes(size = n, fill = n),
    shape = 21,
    color = "white",
    stroke = 1.1,
    alpha = 0.95
  ) +
  geom_text(
    aes(label = ifelse(n > 0, n, "")),
    size = 4.4,
    fontface = "bold",
    color = "white"
  ) +
  scale_size_continuous(
    range = c(3, 18),
    breaks = c(1, 5, 20, 50, 100),
    name = "Records"
  ) +
  scale_fill_gradient(
    low = "#C7B19C",
    high = "#798E87",
    name = "Records"
  ) +
  labs(
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.major = element_line(linewidth = 0.25, color = "grey88"),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13, face = "bold"),
    axis.text.y = element_text(size = 13, face = "bold"),
    legend.position = "right",
    plot.margin = margin(10, 20, 10, 10)
  )

5 4. Topic and Text-Network Analysis

5.1 4.1 Structural Topic Model

Feedback narratives are normalized, tokenized, filtered, converted to a document-feature matrix, and modeled using a structural topic model with incident type as a prevalence covariate.

# Clean data
dat_clean <- dat %>%
  mutate(
    id = row_number(),
    Feedback = as.character(Feedback),
    IncidentType = as.factor(IncidentType),
    text_clean = Feedback %>%
      str_to_lower() %>%
      str_replace_all("[^a-z0-9\\s]", " ") %>%
      str_squish()
  ) %>%
  filter(!is.na(text_clean), text_clean != "") %>%
  filter(!is.na(IncidentType))

# Create quanteda corpus with metadata
corp <- quanteda::corpus(
  dat_clean,
  text_field = "text_clean"
)

# Tokenize
toks <- quanteda::tokens(
  corp,
  remove_punct = TRUE,
  remove_numbers = TRUE,
  remove_symbols = TRUE
)

toks <- quanteda::tokens_remove(
  toks,
  pattern = c(
    quanteda::stopwords("en"),
    "austin", "tx", "st", "street", "rd", "road",
    "vehicle", "vehicles", "car", "cars",
    "said", "stated", "called", "report", "reported",
    "issue", "resident"
  )
)

# Create dfm
dfm_mat <- quanteda::dfm(toks)

# Trim sparse/common words
dfm_mat <- quanteda::dfm_trim(
  dfm_mat,
  min_docfreq = 3,
  max_docfreq = nrow(dat_clean) * 0.80,
  docfreq_type = "count"
)

# IMPORTANT: remove documents with zero remaining words
dfm_mat <- dfm_mat[quanteda::ntoken(dfm_mat) > 0, ]

# Convert to STM
stm_input <- quanteda::convert(dfm_mat, to = "stm")

docs  <- stm_input$documents
vocab <- stm_input$vocab
meta  <- stm_input$meta

# Check dimensions before fitting
length(docs)
## [1] 267
nrow(meta)
## [1] 267
table(meta$IncidentType, useNA = "ifany")
## 
##      Blocking Traffic             Collision  Ignore APD Direction 
##                    46                    13                    26 
##             Near Miss              Nuisance        Safety Concern 
##                    37                    43                    97 
## Safety Issue - School 
##                     5
set.seed(1234)

K_final <- 8

stm_fit <- stm::stm(
  documents = docs,
  vocab = vocab,
  K = K_final,
  prevalence = ~ IncidentType,
  data = meta,
  init.type = "Spectral",
  max.em.its = 100
)
## Beginning Spectral Initialization 
##   Calculating the gram matrix...
##   Finding anchor words...
##      ........
##   Recovering initialization...
##      ..........
## Initialization complete.
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 1 (approx. per word bound = -6.339) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 2 (approx. per word bound = -6.017, relative change = 5.075e-02) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 3 (approx. per word bound = -5.940, relative change = 1.283e-02) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 4 (approx. per word bound = -5.922, relative change = 3.084e-03) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 5 (approx. per word bound = -5.914, relative change = 1.385e-03) 
## Topic 1: traffic, waymo, blocking, lane, stopped 
##  Topic 2: cruise, us, safety, waymo, night 
##  Topic 3: waymo, parking, one, lot, drive 
##  Topic 4: waymo, lot, t, area, driving 
##  Topic 5: stop, waymo, way, sign, intersection 
##  Topic 6: waymo, lane, traffic, turn, left 
##  Topic 7: waymo, traffic, lane, turn, left 
##  Topic 8: traffic, guadalupe, waymo, one, also 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 6 (approx. per word bound = -5.909, relative change = 7.288e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 7 (approx. per word bound = -5.907, relative change = 4.474e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 8 (approx. per word bound = -5.905, relative change = 2.918e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 9 (approx. per word bound = -5.904, relative change = 2.210e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 10 (approx. per word bound = -5.902, relative change = 2.017e-04) 
## Topic 1: traffic, waymo, blocking, lane, stopped 
##  Topic 2: cruise, us, safety, waymo, night 
##  Topic 3: waymo, parking, lot, one, drive 
##  Topic 4: waymo, lot, t, neighborhood, parking 
##  Topic 5: stop, waymo, way, sign, intersection 
##  Topic 6: waymo, lane, traffic, turn, left 
##  Topic 7: waymo, traffic, lane, left, around 
##  Topic 8: traffic, 6th, guadalupe, waymo, one 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 11 (approx. per word bound = -5.901, relative change = 1.915e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 12 (approx. per word bound = -5.900, relative change = 1.845e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 13 (approx. per word bound = -5.899, relative change = 1.705e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 14 (approx. per word bound = -5.898, relative change = 1.423e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 15 (approx. per word bound = -5.898, relative change = 1.040e-04) 
## Topic 1: traffic, waymo, blocking, lane, stopped 
##  Topic 2: cruise, us, safety, waymo, night 
##  Topic 3: waymo, parking, lot, one, blocking 
##  Topic 4: waymo, lot, neighborhood, t, parking 
##  Topic 5: stop, waymo, way, sign, intersection 
##  Topic 6: waymo, lane, traffic, turn, left 
##  Topic 7: waymo, traffic, lane, left, around 
##  Topic 8: traffic, 6th, one, guadalupe, waymo 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 16 (approx. per word bound = -5.897, relative change = 1.026e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 17 (approx. per word bound = -5.896, relative change = 1.336e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 18 (approx. per word bound = -5.896, relative change = 1.231e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 19 (approx. per word bound = -5.895, relative change = 8.050e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 20 (approx. per word bound = -5.895, relative change = 8.271e-05) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, us, safety, waymo, city 
##  Topic 3: waymo, lot, parking, blocking, one 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, sign, intersection 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, traffic, lane, left, around 
##  Topic 8: traffic, 6th, one, guadalupe, cruise 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 21 (approx. per word bound = -5.894, relative change = 9.530e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 22 (approx. per word bound = -5.894, relative change = 1.071e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 23 (approx. per word bound = -5.893, relative change = 1.190e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 24 (approx. per word bound = -5.892, relative change = 1.207e-04) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 25 (approx. per word bound = -5.891, relative change = 1.171e-04) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, us, safety, city, waymo 
##  Topic 3: waymo, lot, parking, blocking, one 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, intersection, sign 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, lane, traffic, left, around 
##  Topic 8: traffic, cruise, one, 6th, guadalupe 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 26 (approx. per word bound = -5.891, relative change = 9.802e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 27 (approx. per word bound = -5.890, relative change = 8.623e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 28 (approx. per word bound = -5.890, relative change = 9.090e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 29 (approx. per word bound = -5.889, relative change = 8.676e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 30 (approx. per word bound = -5.889, relative change = 6.944e-05) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, city, us, safety, waymo 
##  Topic 3: waymo, lot, parking, blocking, intersection 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, intersection, driving 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, lane, traffic, around, left 
##  Topic 8: traffic, cruise, one, 6th, guadalupe 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 31 (approx. per word bound = -5.889, relative change = 4.461e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 32 (approx. per word bound = -5.888, relative change = 3.308e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 33 (approx. per word bound = -5.888, relative change = 3.341e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 34 (approx. per word bound = -5.888, relative change = 4.566e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 35 (approx. per word bound = -5.888, relative change = 4.731e-05) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, city, safety, us, waymo 
##  Topic 3: waymo, lot, parking, intersection, blocking 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, intersection, driving 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, lane, traffic, around, left 
##  Topic 8: traffic, cruise, one, 6th, guadalupe 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 36 (approx. per word bound = -5.888, relative change = 2.852e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 37 (approx. per word bound = -5.887, relative change = 1.392e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 38 (approx. per word bound = -5.887, relative change = 1.055e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 39 (approx. per word bound = -5.887, relative change = 1.463e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 40 (approx. per word bound = -5.887, relative change = 1.890e-05) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, city, safety, us, waymo 
##  Topic 3: waymo, lot, parking, intersection, blocking 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, intersection, driving 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, lane, traffic, around, left 
##  Topic 8: cruise, traffic, one, 6th, guadalupe 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 41 (approx. per word bound = -5.887, relative change = 2.086e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 42 (approx. per word bound = -5.887, relative change = 1.914e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 43 (approx. per word bound = -5.887, relative change = 1.969e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 44 (approx. per word bound = -5.887, relative change = 2.426e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 45 (approx. per word bound = -5.887, relative change = 2.556e-05) 
## Topic 1: traffic, waymo, blocking, minutes, lane 
##  Topic 2: cruise, city, safety, waymo, us 
##  Topic 3: waymo, lot, parking, intersection, blocking 
##  Topic 4: waymo, neighborhood, lot, t, parking 
##  Topic 5: stop, waymo, way, intersection, driving 
##  Topic 6: waymo, lane, turn, traffic, left 
##  Topic 7: waymo, lane, traffic, around, left 
##  Topic 8: cruise, traffic, one, 6th, guadalupe 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 46 (approx. per word bound = -5.886, relative change = 1.622e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Completing Iteration 47 (approx. per word bound = -5.886, relative change = 1.632e-05) 
## .....................................................................................................................................
## Completed E-Step (0 seconds). 
## Completed M-Step. 
## Model Converged
stm::labelTopics(stm_fit, n = 10)
## Topic 1 Top Words:
##       Highest Prob: traffic, waymo, blocking, minutes, lane, stopped, get, officer, t, driving 
##       FREX: game, minutes, chavez, support, blocking, cesar, officers, 15th, river, impeding 
##       Lift: ridiculous, taxi, cannon, football, game, major, site, william, impeding, eb 
##       Score: football, officer, game, officers, impeding, blocking, river, minutes, control, technician 
## Topic 2 Top Words:
##       Highest Prob: cruise, city, safety, waymo, us, incident, parked, night, front, close 
##       FREX: safety, city, curb, speed, response, cruise, close, afd, night, limit 
##       Lift: according, council, difficult, management, states, dogs, present, appropriately, row, state 
##       Score: present, cruise, delays, city, issued, curb, dogs, management, small, seriously 
## Topic 3 Top Words:
##       Highest Prob: waymo, lot, parking, intersection, blocking, scene, one, drive, photo, shows 
##       FREX: scene, police, shows, photo, driveway, parking, responded, lot, signals, arrived 
##       Lift: alley, became, damage, escorted, later, motorcycle, runners, inappropriate, closing, procession 
##       Score: became, police, signals, procession, scene, parking, arrived, lot, damage, blocking 
## Topic 4 Top Words:
##       Highest Prob: waymo, neighborhood, lot, t, parking, can, area, s, driving, around 
##       FREX: neighborhood, quiet, valet, day, residential, park, picking, kids, property, lot 
##       Lift: additionally, advise, express, loud, requesting, homes, kids, local, main, presence 
##       Score: kids, neighborhood, residential, property, local, picking, school, quiet, using, dropping 
## Topic 5 Top Words:
##       Highest Prob: stop, waymo, way, intersection, driving, sign, s, crossing, zoox, just 
##       FREX: zoox, crossing, sign, crosswalk, stop, concern, mph, guard, neighborhoods, bike 
##       Lift: creates, xfx, unpredictable, involved, riding, zoox, neighborhoods, crosswalk, bicycle, central 
##       Score: creates, sign, zoox, guard, neighborhoods, mph, concern, stop, crosswalk, crossing 
## Topic 6 Top Words:
##       Highest Prob: waymo, lane, turn, traffic, left, right, light, stopped, green, turning 
##       FREX: lamar, green, turning, turn, left, oncoming, right, eastbound, bound, behind 
##       Lift: law, rather, 38th, cutting, poles, woodward, n, corner, observe, october 
##       Score: cutting, turn, lamar, officer, oncoming, bound, light, left, 38th, turning 
## Topic 7 Top Words:
##       Highest Prob: waymo, lane, traffic, around, left, turn, stopped, red, stop, light 
##       FREX: burnet, went, northbound, flashing, red, onto, new, lights, go, truck 
##       Lift: b, remember, clyde, littlefield, abruptly, burnet, confirmed, avoided, warning, new 
##       Score: remember, burnet, barrels, new, clyde, littlefield, knocked, turn, confirmed, cuts 
## Topic 8 Top Words:
##       Highest Prob: cruise, traffic, one, 6th, guadalupe, intersection, waymo, streets, closure, san 
##       FREX: guadalupe, 6th, closure, campus, allowed, confused, san, 11th, jacinto, streets 
##       Lift: 11th, guadalupe, campus, 6th, 12th, closure, escalated, caller, waterloo, confused 
##       Score: 11th, guadalupe, closure, 6th, campus, cruise, confused, operator, remote, lt

5.2 4.2 Text Normalization and Feature Construction

dat_clean <- dat %>%
  mutate(
    id = row_number(),
    Feedback = as.character(Feedback),
    IncidentType = as.factor(IncidentType),
    text_clean = Feedback %>%
      str_to_lower() %>%
      str_replace_all("cruise|waymo|google self driving|driverless car|driverless cars|autonomous vehicle|autonomous vehicles", "av") %>%
      str_replace_all("afd|fire department|fire engine", "fire_response") %>%
      str_replace_all("apd|police department|officer", "police_response") %>%
      str_replace_all("ems|ambulance", "ems_response") %>%
      str_replace_all("[^a-z0-9\\s]", " ") %>%
      str_squish()
  ) %>%
  filter(!is.na(text_clean), text_clean != "")


custom_stop <- c(
  stopwords("en"),
  "austin", "tx", "st", "street", "rd", "road", "ave", "blvd",
  "vehicle", "vehicles", "car", "cars", "av",
  "said", "stated", "called", "report", "reported", "issue",
  "resident", "caller", "would", "could", "also", "one", "two"
)

corp <- quanteda::corpus(dat_clean, text_field = "text_clean")

toks <- quanteda::tokens(
  corp,
  remove_punct = TRUE,
  remove_numbers = TRUE,
  remove_symbols = TRUE
) %>%
  quanteda::tokens_remove(custom_stop) %>%
  quanteda::tokens_wordstem()

dfm_av <- quanteda::dfm(toks)

dfm_av <- quanteda::dfm_trim(
  dfm_av,
  min_docfreq = 3,
  docfreq_type = "count"
)

dfm_av
## Document-feature matrix of: 267 documents, 899 features (95.95% sparse) and 16 docvars.
##        features
## docs    time 6th let sever s area make question move wit
##   text1    1   2   1     2 1    1    1        1    2   1
##   text2    1   0   0     1 1    0    2        0    0   0
##   text3    0   0   0     0 0    0    0        0    0   0
##   text4    0   0   0     0 0    0    0        0    0   0
##   text5    2   0   0     1 0    0    0        0    0   0
##   text6    0   0   0     0 0    0    1        0    0   0
## [ reached max_ndoc ... 261 more documents, reached max_nfeat ... 889 more features ]
# Get top words from dfm, not from fcm
top_words <- names(sort(quanteda::colSums(dfm_av), decreasing = TRUE))[1:60]

# Create feature co-occurrence matrix
fcm_av <- quanteda::fcm(
  toks,
  context = "window",
  window = 5,
  count = "frequency"
)

# Select those top words in the fcm
fcm_trim <- quanteda::fcm_select(
  fcm_av,
  pattern = top_words,
  selection = "keep"
)

5.3 4.3 Quick Word Co-occurrence Network

# Get top words from dfm instead of fcm
top_words <- names(sort(quanteda::colSums(dfm_av), decreasing = TRUE))[1:60]

# Trim fcm to top words
fcm_trim <- quanteda::fcm_select(
  fcm_av,
  pattern = top_words,
  selection = "keep"
)

# Plot network
quanteda.textplots::textplot_network(
  fcm_trim,
  min_freq = 5,
  vertex_color = "#798E87",
  vertex_labelcolor = "black",
  edge_color = "#9C964A",
  edge_alpha = 0.6,
  edge_size = 2,
  vertex_labelsize = 6)

6 5. Policy Analysis 1: Interpretable Civic-Risk Screening

Policy Analysis 1 translates textual feedback into transparent risk domains, criticality classes, review priorities, and recommended policy responses.

6.1 5.1 Initial Policy-Risk Flags and Priority Scores

policy_cols <- c(
  "#D8A499",
  "#798E87",
  "#CCC591",
  "#CDC08C",
  "#C27D38"
)

dat_policy <- dat_clean %>%
  mutate(
    text_lower = str_to_lower(text_clean),

    emergency_flag = str_detect(text_lower, "fire|afd|ems|ambulance|police|apd|officer|emergency|code 3"),
    blockage_flag  = str_detect(text_lower, "block|blocking|stopped|lane|traffic|closure|congest"),
    VRU_flag       = str_detect(text_lower, "pedestrian|cyclist|bike|bicycle|child|children|stroller|walking"),
    residential_flag = str_detect(text_lower, "neighborhood|home|porch|sleep|night|noise|family|daughter"),
    privacy_flag   = str_detect(text_lower, "camera|spyware|privacy|surveillance|consent|experiment"),
    collision_flag = IncidentType == "Collision",
    near_miss_flag = IncidentType == "Near Miss",

    policy_priority_score =
      5 * collision_flag +
      4 * near_miss_flag +
      4 * emergency_flag +
      3 * VRU_flag +
      3 * blockage_flag +
      2 * residential_flag +
      1 * privacy_flag
  )

risk_long <- dat_policy %>%
  select(
    IncidentType,
    emergency_flag,
    blockage_flag,
    VRU_flag,
    residential_flag,
    privacy_flag,
    collision_flag,
    near_miss_flag
  ) %>%
  pivot_longer(
    cols = -IncidentType,
    names_to = "Risk_Domain",
    values_to = "Flag"
  ) %>%
  filter(Flag == TRUE) %>%
  count(IncidentType, Risk_Domain, name = "n") %>%
  mutate(
    Risk_Domain = Risk_Domain %>%
      str_remove("_flag") %>%
      str_replace_all("_", " ") %>%
      str_to_title(),
    IncidentType = fct_reorder(IncidentType, n, .fun = sum)
  )

ggplot(risk_long, aes(x = Risk_Domain, y = IncidentType, fill = n)) +
  geom_tile(color = "white", linewidth = 0.6) +
  geom_text(aes(label = n), size = 3.8, color = "black") +
  scale_fill_gradientn(
    colors = policy_cols
  ) +
  labs(
    title = "",
    subtitle = "",
    x = NULL,
    y = NULL,
    fill = "Records"
  ) +
  theme_minimal(base_size = 18) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    axis.text.x = element_text(angle = 35, hjust = 1),
    panel.grid = element_blank(),
    legend.position = "right"
  )

#### bubble

ggplot(risk_long, aes(x = Risk_Domain, y = IncidentType)) +
  geom_point(
    aes(size = n, fill = n),
    shape = 21,
    color = "gray25",
    alpha = 0.90,
    stroke = 0.35
  ) +
  scale_size_continuous(
    range = c(3, 13)
  ) +
  scale_fill_gradientn(
    colors = policy_cols
  ) +
  labs(
    title = "",
    subtitle = "",
    x = NULL,
    y = NULL,
    size = "Records",
    fill = "Records"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 11),
    axis.text.x = element_text(angle = 35, hjust = 1),
    panel.grid.major = element_line(color = "gray90"),
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

priority_summary <- dat_policy %>%
  group_by(IncidentType) %>%
  summarise(
    mean_priority = mean(policy_priority_score, na.rm = TRUE),
    max_priority = max(policy_priority_score, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_priority))

ggplot(priority_summary,
       aes(x = reorder(IncidentType, mean_priority),
           y = mean_priority)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Average Policy Priority Score by Incident Type",
    x = NULL,
    y = "Mean policy priority score"
  ) +
  theme_minimal(base_size = 13)

6.2 5.2 Expanded Civic-Risk Taxonomy

set.seed(1234)


dat <- dat %>%
  mutate(
    id = row_number(),
    Feedback = as.character(Feedback),
    IncidentType = as.factor(IncidentType)
  ) %>%
  filter(!is.na(Feedback), Feedback != "")

table(dat$IncidentType, useNA = "ifany")
## 
##      Blocking Traffic             Collision  Ignore APD Direction 
##                    46                    13                    26 
##             Near Miss              Nuisance        Safety Concern 
##                    37                    43                    97 
## Safety Issue - School 
##                     5
dat_clean <- dat %>%
  mutate(
    text_clean = Feedback %>%
      str_to_lower() %>%
      str_replace_all("cruise|waymo|google self-driving|google self driving", "av_operator") %>%
      str_replace_all("driverless cars|driverless car|self-driving cars|self driving cars|autonomous vehicles|autonomous vehicle", "av") %>%
      str_replace_all("afd|fire department|fire engine|fire truck|engine", "fire_response") %>%
      str_replace_all("apd|police department|police officer|officer", "police_response") %>%
      str_replace_all("ems|ambulance", "ems_response") %>%
      str_replace_all("moody center", "moody_center") %>%
      str_replace_all("near miss", "near_miss") %>%
      str_replace_all("blocking traffic", "blocking_traffic") %>%
      str_replace_all("code 3", "code_3") %>%
      str_replace_all("[^a-z0-9_\\s]", " ") %>%
      str_squish()
  )


dat_risk <- dat_clean %>%
  mutate(
    emergency_response = str_detect(
      text_clean,
      "fire_response|ems_response|police_response|ambulance|emergency|code_3|lights|sirens|apparatus|traffic accident"
    ),

    traffic_obstruction = str_detect(
      text_clean,
      "block|blocking|blocked|stopped|stop|lane|traffic|congest|closure|bunch|middle lane|hazards|flashers"
    ),

    vru_risk = str_detect(
      text_clean,
      "pedestrian|cyclist|bicyclist|bicycle|bike|child|children|stroller|walking|walk|crossing"
    ),

    residential_disruption = str_detect(
      text_clean,
      "neighborhood|home|porch|sleep|night|noise|quiet|family|daughter|residential|street i live"
    ),

    privacy_consent = str_detect(
      text_clean,
      "camera|spyware|privacy|surveillance|consent|experiment|training data|private company"
    ),

    manual_control_failure = str_detect(
      text_clean,
      "directing traffic|signals|commands|against what|traffic control|officer|police_response|hand signal|flag"
    ),

    collision_recovery = str_detect(
      text_clean,
      "collision|collided|crash|crashed|hit|struck|recovery mode|tow truck|unable to move|damaged"
    ),

    nighttime_operation = str_detect(
      text_clean,
      "night|2 am|3 am|0200|late|evening|6p|6 pm|6 a|6a"
    )
  )


dat_risk <- dat_risk %>%
  mutate(
    policy_priority_score =
      6 * (IncidentType == "Collision") +
      5 * (IncidentType == "Near Miss") +
      5 * emergency_response +
      4 * manual_control_failure +
      4 * collision_recovery +
      3 * traffic_obstruction +
      3 * vru_risk +
      2 * residential_disruption +
      2 * nighttime_operation +
      1 * privacy_consent,

    criticality_class = case_when(
      policy_priority_score >= 10 ~ "Urgent",
      policy_priority_score >= 6  ~ "High",
      policy_priority_score >= 3  ~ "Moderate",
      TRUE ~ "Low"
    )
  )

dat_risk %>%
  count(criticality_class, sort = TRUE)
## # A tibble: 4 × 2
##   criticality_class     n
##   <chr>             <int>
## 1 Urgent              115
## 2 Moderate             66
## 3 High                 65
## 4 Low                  21
dat_risk <- dat_risk %>%
  mutate(
    dominant_risk_domain = case_when(
      collision_recovery ~ "Collision / recovery failure",
      emergency_response ~ "Emergency-response interference",
      manual_control_failure ~ "Manual traffic-control failure",
      vru_risk ~ "Vulnerable road-user risk",
      traffic_obstruction ~ "Traffic obstruction",
      residential_disruption ~ "Residential disruption",
      privacy_consent ~ "Privacy / consent concern",
      TRUE ~ "General AV concern"
    )
  )

dat_risk %>%
  count(dominant_risk_domain, sort = TRUE)
## # A tibble: 7 × 2
##   dominant_risk_domain                n
##   <chr>                           <int>
## 1 Emergency-response interference    82
## 2 Traffic obstruction                64
## 3 Collision / recovery failure       60
## 4 Vulnerable road-user risk          28
## 5 General AV concern                 22
## 6 Residential disruption              6
## 7 Manual traffic-control failure      5
policy_actions <- tibble(
  dominant_risk_domain = c(
    "Collision / recovery failure",
    "Emergency-response interference",
    "Manual traffic-control failure",
    "Vulnerable road-user risk",
    "Traffic obstruction",
    "Residential disruption",
    "Privacy / consent concern",
    "General AV concern"
  ),
  recommended_policy_action = c(
    "Require rapid removal protocol, remote-operator response standard, and responder-access procedure.",
    "Require emergency vehicle interaction protocol and responder override/communication pathway.",
    "Restrict AV operations during manually directed traffic, events, and temporary traffic-control conditions.",
    "Apply VRU-sensitive ODD limits near pedestrian, bicycle, school, campus, and nightlife zones.",
    "Require minimum-risk maneuver standards that prevent lane blockage and unsafe stopping.",
    "Apply nighttime residential routing limits, notification rules, and neighborhood geofencing.",
    "Require public transparency, data-use disclosure, and camera/privacy governance.",
    "Monitor and classify for future escalation."
  )
)

dat_policy <- dat_risk %>%
  left_join(policy_actions, by = "dominant_risk_domain")


policy_actions <- tibble(
  dominant_risk_domain = c(
    "Collision / recovery failure",
    "Emergency-response interference",
    "Manual traffic-control failure",
    "Vulnerable road-user risk",
    "Traffic obstruction",
    "Residential disruption",
    "Privacy / consent concern",
    "General AV concern"
  ),
  recommended_policy_action = c(
    "Require rapid removal protocol, remote-operator response standard, and responder-access procedure.",
    "Require emergency vehicle interaction protocol and responder override/communication pathway.",
    "Restrict AV operations during manually directed traffic, events, and temporary traffic-control conditions.",
    "Apply VRU-sensitive ODD limits near pedestrian, bicycle, school, campus, and nightlife zones.",
    "Require minimum-risk maneuver standards that prevent lane blockage and unsafe stopping.",
    "Apply nighttime residential routing limits, notification rules, and neighborhood geofencing.",
    "Require public transparency, data-use disclosure, and camera/privacy governance.",
    "Monitor and classify for future escalation."
  )
)

dat_policy <- dat_risk %>%
  left_join(policy_actions, by = "dominant_risk_domain")

6.3 5.3 Latent Semantic Representation and UMAP Visualization

tokens <- dat_policy$text_clean %>%
  word_tokenizer()

it <- itoken(
  tokens,
  ids = dat_policy$id,
  progressbar = FALSE
)

vocab <- create_vocabulary(it) %>%
  prune_vocabulary(
    term_count_min = 3,
    doc_proportion_max = 0.80
  )

vectorizer <- vocab_vectorizer(vocab)

dtm <- create_dtm(it, vectorizer)

tfidf <- TfIdf$new()
dtm_tfidf <- fit_transform(dtm, tfidf)

lsa <- LSA$new(n_topics = 30)
doc_embeddings <- fit_transform(dtm_tfidf, lsa)


umap_out <- uwot::umap(
  doc_embeddings,
  n_neighbors = 15,
  min_dist = 0.12,
  metric = "cosine",
  scale = TRUE
)

umap_df <- dat_policy %>%
  bind_cols(
    tibble(
      UMAP1 = umap_out[, 1],
      UMAP2 = umap_out[, 2]
    )
  )

ggplot(
  umap_df,
  aes(
    x = UMAP1,
    y = UMAP2,
    color = dominant_risk_domain,
    shape = criticality_class
  )
) +
  geom_point(size = 3, alpha = 0.85) +
  labs(
    title = "",
    subtitle = "",
    x = "UMAP dimension 1",
    y = "UMAP dimension 2",
    color = "Risk domain",
    shape = "Criticality"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    legend.position = "right",
    plot.title = element_text(face = "bold")
  )

6.4 5.4 Priority Retrieval and Risk-Domain Visualization

priority_curve <- dat_policy %>%
  arrange(desc(policy_priority_score)) %>%
  mutate(
    review_rank = row_number(),
    high_risk = criticality_class %in% c("Urgent", "High"),
    cumulative_high_risk = cumsum(high_risk),
    total_high_risk = sum(high_risk),
    recall_high_risk = cumulative_high_risk / total_high_risk,
    review_share = review_rank / n()
  )

ggplot(priority_curve, aes(x = review_share, y = recall_high_risk)) +
  geom_line(linewidth = 1.2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  scale_x_continuous(labels = scales::percent) +
  scale_y_continuous(labels = scales::percent) +
  labs(
    title = "",
    subtitle = "Higher curve indicates faster retrieval of urgent and high-risk cases",
    x = "Share of feedback reviewed",
    y = "Share of urgent/high-risk feedback captured"
  ) +
  theme_minimal(base_size = 16)

risk_heat <- dat_policy %>%
  count(IncidentType, dominant_risk_domain) %>%
  group_by(IncidentType) %>%
  mutate(percent = 100 * n / sum(n)) %>%
  ungroup()

ggplot(
  risk_heat,
  aes(
    x = IncidentType,
    y = dominant_risk_domain,
    fill = percent
  )
) +
  geom_tile(color = "white") +
  geom_text(aes(label = paste0(round(percent, 1), "%")), size = 3.5) +
  labs(
    title = "Operational Risk Domains by AV Incident Type",
    x = "Incident type",
    y = "Operational risk domain",
    fill = "Percent"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1),
    plot.title = element_text(face = "bold")
  )

risk_cols <- c(
  "#D8A499",
  "#CCC591",
  "#CDC08C",
  "#798E87",
  "#C27D38"
)

risk_heat <- risk_heat %>%
  mutate(
    IncidentType = forcats::fct_reorder(IncidentType, percent, .fun = max),
    dominant_risk_domain = forcats::fct_reorder(dominant_risk_domain, percent, .fun = max)
  )

ggplot(risk_heat, aes(x = IncidentType, y = dominant_risk_domain, fill = percent)) +
  geom_tile(color = "white", linewidth = 0.7) +
  geom_text(aes(label = paste0(round(percent, 1), "%")), size = 3.5) +
  scale_fill_gradientn(
    colors = risk_cols,
    labels = function(x) paste0(x, "%")
  ) +
  labs(
    title = "",
    subtitle = "",
    x = NULL,
    y = NULL,
    fill = "Percent"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    axis.text.x = element_text(angle = 35, hjust = 1),
    plot.title = element_text(face = "bold", size = 15),
    panel.grid = element_blank()
  )

6.5 5.5 High-Risk Feedback Classification and Review Triage

model_df <- dat_policy %>%
  mutate(
    high_risk_label = if_else(
      criticality_class %in% c("Urgent", "High"),
      1,
      0
    )
  )

x <- as.matrix(doc_embeddings)
y <- model_df$high_risk_label

set.seed(1234)

train_index <- createDataPartition(y, p = 0.75, list = FALSE)

x_train <- x[train_index, ]
x_test  <- x[-train_index, ]

y_train <- y[train_index]
y_test  <- y[-train_index]

cv_fit <- cv.glmnet(
  x_train,
  y_train,
  family = "binomial",
  alpha = 0.5,
  type.measure = "auc"
)

pred_prob <- as.numeric(
  predict(cv_fit, newx = x_test, s = "lambda.min", type = "response")
)

roc_obj <- roc(y_test, pred_prob)

auc(roc_obj)
## Area under the curve: 0.7073
plot(
  roc_obj,
  main = "High-Risk AV Feedback Classifier ROC Curve"
)

all_pred_prob <- as.numeric(
  predict(cv_fit, newx = x, s = "lambda.min", type = "response")
)

dat_policy <- dat_policy %>%
  mutate(
    model_high_risk_prob = all_pred_prob,
    model_uncertainty = abs(model_high_risk_prob - 0.50),
    review_priority = case_when(
      model_high_risk_prob >= 0.80 ~ "Auto-escalate",
      model_high_risk_prob >= 0.60 ~ "Priority human review",
      model_uncertainty <= 0.10 ~ "Uncertain, needs review",
      TRUE ~ "Routine monitoring"
    )
  )

dat_policy %>%
  count(review_priority, sort = TRUE)
## # A tibble: 4 × 2
##   review_priority             n
##   <chr>                   <int>
## 1 Auto-escalate             116
## 2 Priority human review      77
## 3 Routine monitoring         47
## 4 Uncertain, needs review    27
review_cols <- c(
  "#D8A499",
  "#798E87",
  "#CCC591",
  "#CDC08C",
  "#C27D38"
)

ggplot(
  dat_policy,
  aes(
    x = model_high_risk_prob,
    fill = review_priority
  )
) +
  geom_histogram(
    bins = 25,
    alpha = 0.85,
    color = "white",
    linewidth = 0.25
  ) +
  scale_fill_manual(
    values = review_cols
  ) +
  labs(
    title = "",
    subtitle = "",
    x = "Predicted probability of high-risk feedback",
    y = "Number of feedback records",
    fill = "Review category"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    plot.title = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )

ggsave(
  "Fig_uncertainty_triage.png",
  width = 9,
  height = 6,
  dpi = 400
)

6.6 5.6 Policy Summaries and Excel Export

top_escalation <- dat_policy %>%
  arrange(desc(policy_priority_score), desc(model_high_risk_prob)) %>%
  select(
    id,
    IncidentType,
    criticality_class,
    policy_priority_score,
    model_high_risk_prob,
    review_priority,
    dominant_risk_domain,
    recommended_policy_action,
    Feedback
  )

risk_summary <- dat_policy %>%
  group_by(dominant_risk_domain, recommended_policy_action) %>%
  summarise(
    n_records = n(),
    mean_priority = mean(policy_priority_score, na.rm = TRUE),
    urgent_high_count = sum(criticality_class %in% c("Urgent", "High")),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_priority), desc(n_records))

incident_summary <- dat_policy %>%
  group_by(IncidentType) %>%
  summarise(
    n_records = n(),
    mean_priority = mean(policy_priority_score, na.rm = TRUE),
    urgent_high_count = sum(criticality_class %in% c("Urgent", "High")),
    mean_model_high_risk_prob = mean(model_high_risk_prob, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_priority))

write_xlsx(
  list(
    policy_ready_data = dat_policy,
    top_escalation_cases = top_escalation,
    risk_domain_summary = risk_summary,
    incident_type_summary = incident_summary,
    policy_actions = policy_actions
  ),
  "AV_CivicRisk_policy_outputs.xlsx"
)

6.7 5.7 Spatial Distribution of Policy Risk

if (all(c("Latitude", "Longitude") %in% names(dat_policy))) {

  spatial_risk <- dat_policy %>%
    filter(!is.na(Latitude), !is.na(Longitude))

  ggplot(
    spatial_risk,
    aes(
      x = Longitude,
      y = Latitude,
      color = dominant_risk_domain,
      size = policy_priority_score
    )
  ) +
    geom_point(alpha = 0.75) +
    labs(
      title = "Spatial Distribution of AV Operational Risk Reports",
      x = "Longitude",
      y = "Latitude",
      color = "Risk domain",
      size = "Priority score"
    ) +
    theme_minimal(base_size = 13)

  ggsave(
    "Fig_spatial_AV_risk_points.png",
    width = 9,
    height = 7,
    dpi = 400
  )
}


risk_cols <- c(
  "#D8A499",
  "#798E87",
  "#CCC591",
  "#CDC08C",
  "#C27D38",
  "#C7B19C",
  "#A2A475"
)

risk_levels <- sort(unique(spatial_risk$dominant_risk_domain))

risk_cols_named <- setNames(
  risk_cols[seq_along(risk_levels)],
  risk_levels
)

ggplot(
  spatial_risk,
  aes(
    x = Longitude,
    y = Latitude,
    color = dominant_risk_domain,
    size = policy_priority_score
  )
) +
  geom_point(alpha = 0.75) +
  scale_color_manual(values = risk_cols_named) +
  scale_size_continuous(range = c(2, 8)) +
  labs(
    x = "Longitude",
    y = "Latitude",
    color = "Risk domain",
    size = "Priority score"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

7 6. Policy Analysis 2: Advanced Hybrid Modeling and Decision Support

Policy Analysis 2 compares rule-based, metadata-based, text-embedding, and hybrid models; evaluates discrimination and calibration; implements conformal triage; simulates policy coverage; and summarizes spatial hotspots.

7.1 6.1 Colors, Outcomes, and Structured Features

# ------------------------------------------------------------
# 0. Colors
# ------------------------------------------------------------

risk_cols <- c(
  "#D8A499",
  "#798E87",
  "#CCC591",
  "#CDC08C",
  "#C27D38",
  "#C7B19C",
  "#A2A475"
)

risk_levels <- sort(unique(dat_policy$dominant_risk_domain))

risk_cols_named <- setNames(
  risk_cols[seq_along(risk_levels)],
  risk_levels
)

review_cols <- c(
  "Auto-escalate" = "#C27D38",
  "Priority human review" = "#CDC08C",
  "Uncertain, needs review" = "#D8A499",
  "Routine monitoring" = "#798E87"
)

# ------------------------------------------------------------
# 1. Create high-risk label and structured features
# ------------------------------------------------------------

model_df <- dat_policy %>%
  mutate(
    high_risk_label = if_else(
      criticality_class %in% c("Urgent", "High"),
      1L, 0L
    ),
    high_risk_factor = factor(
      high_risk_label,
      levels = c(0, 1),
      labels = c("LowModerate", "UrgentHigh")
    ),
    IncidentType = as.factor(IncidentType),
    dominant_risk_domain = as.factor(dominant_risk_domain)
  )

# LSA embeddings should already exist from your earlier code:
# doc_embeddings <- fit_transform(dtm_tfidf, lsa)

embed_df <- as_tibble(doc_embeddings) %>%
  setNames(paste0("LSA_", seq_len(ncol(doc_embeddings))))

model_all <- bind_cols(
  model_df,
  embed_df
)

7.2 6.2 Train-Test Split and Evaluation Function

# ------------------------------------------------------------
# 2. Train/test split
# ------------------------------------------------------------

set.seed(1234)

split_obj <- initial_split(
  model_all,
  prop = 0.75,
  strata = high_risk_factor
)

train_dat <- training(split_obj)
test_dat  <- testing(split_obj)

# ------------------------------------------------------------
# 3. Helper function for model evaluation
# ------------------------------------------------------------

get_binary_metrics <- function(truth, prob) {

  truth_num <- ifelse(truth == "UrgentHigh", 1, 0)

  roc_obj <- pROC::roc(
    response = truth_num,
    predictor = prob,
    quiet = TRUE
  )

  pr_obj <- PRROC::pr.curve(
    scores.class0 = prob[truth_num == 1],
    scores.class1 = prob[truth_num == 0],
    curve = TRUE
  )

  pred_class <- ifelse(prob >= 0.50, 1, 0)

  tibble(
    AUC = as.numeric(pROC::auc(roc_obj)),
    PRAUC = as.numeric(pr_obj$auc.integral),
    Brier = mean((prob - truth_num)^2),
    Accuracy = mean(pred_class == truth_num),
    Precision = ifelse(sum(pred_class == 1) == 0, NA_real_,
                       sum(pred_class == 1 & truth_num == 1) / sum(pred_class == 1)),
    Recall = ifelse(sum(truth_num == 1) == 0, NA_real_,
                    sum(pred_class == 1 & truth_num == 1) / sum(truth_num == 1))
  )
}

7.3 6.3 Benchmark and Hybrid Models

# ------------------------------------------------------------
# 4. Model 1: Rule-score only
# ------------------------------------------------------------

rule_fit <- glm(
  high_risk_factor ~ policy_priority_score,
  data = train_dat,
  family = binomial()
)

rule_prob <- predict(
  rule_fit,
  newdata = test_dat,
  type = "response"
)

rule_metrics <- get_binary_metrics(
  truth = test_dat$high_risk_factor,
  prob = rule_prob
) %>%
  mutate(Model = "Rule score only")

# ------------------------------------------------------------
# 5. Model 2: Metadata only
# ------------------------------------------------------------

meta_fit <- glm(
  high_risk_factor ~ IncidentType + dominant_risk_domain,
  data = train_dat,
  family = binomial()
)

meta_prob <- predict(
  meta_fit,
  newdata = test_dat,
  type = "response"
)

meta_metrics <- get_binary_metrics(
  truth = test_dat$high_risk_factor,
  prob = meta_prob
) %>%
  mutate(Model = "Metadata only")

# ------------------------------------------------------------
# 6. Model 3: Text embedding only, elastic net
# ------------------------------------------------------------

lsa_cols <- grep("^LSA_", names(model_all), value = TRUE)

x_train_lsa <- as.matrix(train_dat[, lsa_cols])
x_test_lsa  <- as.matrix(test_dat[, lsa_cols])

y_train <- ifelse(train_dat$high_risk_factor == "UrgentHigh", 1, 0)
y_test  <- ifelse(test_dat$high_risk_factor == "UrgentHigh", 1, 0)

set.seed(1234)

cv_lsa <- cv.glmnet(
  x = x_train_lsa,
  y = y_train,
  family = "binomial",
  alpha = 0.5,
  type.measure = "auc"
)

lsa_prob <- as.numeric(
  predict(cv_lsa, newx = x_test_lsa, s = "lambda.min", type = "response")
)

lsa_metrics <- get_binary_metrics(
  truth = test_dat$high_risk_factor,
  prob = lsa_prob
) %>%
  mutate(Model = "Text embeddings only")

# ------------------------------------------------------------
# 7. Model 4: Hybrid model
# Text embeddings + metadata + rule score
# ------------------------------------------------------------

hybrid_train <- train_dat %>%
  select(
    high_risk_factor,
    policy_priority_score,
    IncidentType,
    dominant_risk_domain,
    all_of(lsa_cols)
  )

hybrid_test <- test_dat %>%
  select(
    high_risk_factor,
    policy_priority_score,
    IncidentType,
    dominant_risk_domain,
    all_of(lsa_cols)
  )

x_train_hybrid <- model.matrix(
  high_risk_factor ~ .,
  data = hybrid_train
)[, -1]

x_test_hybrid <- model.matrix(
  high_risk_factor ~ .,
  data = hybrid_test
)[, -1]

set.seed(1234)

cv_hybrid <- cv.glmnet(
  x = x_train_hybrid,
  y = y_train,
  family = "binomial",
  alpha = 0.5,
  type.measure = "auc"
)

hybrid_prob <- as.numeric(
  predict(cv_hybrid, newx = x_test_hybrid, s = "lambda.min", type = "response")
)

hybrid_metrics <- get_binary_metrics(
  truth = test_dat$high_risk_factor,
  prob = hybrid_prob
) %>%
  mutate(Model = "Hybrid policy-text model")

7.4 6.4 Model Comparison

# ------------------------------------------------------------
# 8. Model comparison table
# ------------------------------------------------------------

model_comparison <- bind_rows(
  rule_metrics,
  meta_metrics,
  lsa_metrics,
  hybrid_metrics
) %>%
  select(Model, AUC, PRAUC, Brier, Accuracy, Precision, Recall) %>%
  arrange(desc(AUC))

model_comparison
## # A tibble: 4 × 7
##   Model                      AUC PRAUC    Brier Accuracy Precision Recall
##   <chr>                    <dbl> <dbl>    <dbl>    <dbl>     <dbl>  <dbl>
## 1 Rule score only          1     1     3.25e-19    1         1      1    
## 2 Hybrid policy-text model 1     1     2.12e- 1    0.672     0.672  1    
## 3 Metadata only            0.971 0.988 4.14e- 2    0.955     0.977  0.956
## 4 Text embeddings only     0.685 0.755 2.01e- 1    0.716     0.717  0.956
model_comparison_long <- model_comparison %>%
  pivot_longer(
    cols = c(AUC, PRAUC, Accuracy, Precision, Recall),
    names_to = "Metric",
    values_to = "Value"
  )

ggplot(
  model_comparison_long,
  aes(
    x = Value,
    y = fct_reorder(Model, Value),
    fill = Metric
  )
) +
  geom_col(
    position = position_dodge(width = 0.75),
    width = 0.65
  ) +
  scale_fill_manual(
    values = c(
      "AUC" = "#C27D38",
      "PRAUC" = "#798E87",
      "Accuracy" = "#CCC591",
      "Precision" = "#CDC08C",
      "Recall" = "#D8A499"
    )
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = percent_format(accuracy = 1)
  ) +
  labs(
    x = "Metric value",
    y = NULL,
    fill = "Metric"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    legend.position = "bottom"
  )

7.5 6.5 Calibration and Reliability Assessment

calibration_df <- tibble(
  truth = y_test,
  pred_prob = hybrid_prob
) %>%
  mutate(
    bin = cut(
      pred_prob,
      breaks = seq(0, 1, by = 0.10),
      include.lowest = TRUE
    )
  ) %>%
  group_by(bin) %>%
  summarise(
    mean_pred = mean(pred_prob, na.rm = TRUE),
    observed_rate = mean(truth, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  filter(!is.na(mean_pred), !is.na(observed_rate))

ECE <- calibration_df %>%
  mutate(weight = n / sum(n)) %>%
  summarise(
    ECE = sum(weight * abs(mean_pred - observed_rate))
  ) %>%
  pull(ECE)

Brier_hybrid <- mean((hybrid_prob - y_test)^2)

ECE
## [1] 0.01193531
Brier_hybrid
## [1] 0.2117
ggplot(
  calibration_df,
  aes(
    x = mean_pred,
    y = observed_rate,
    size = n
  )
) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed",
    color = "gray45"
  ) +
  geom_point(
    shape = 21,
    fill = "#C27D38",
    color = "gray20",
    alpha = 0.85
  ) +
  geom_line(
    linewidth = 1.1,
    color = "#798E87"
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = percent_format(accuracy = 1)
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = percent_format(accuracy = 1)
  ) +
  scale_size_continuous(range = c(3, 10)) +
  labs(
    x = "Mean predicted high-risk probability",
    y = "Observed high-risk rate",
    size = "Records"
  ) +
  annotate(
    "text",
    x = 0.65,
    y = 0.15,
    label = paste0(
      "ECE = ", round(ECE, 3),
      "\nBrier = ", round(Brier_hybrid, 3)
    ),
    size = 5,
    hjust = 0
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

ggplot(
  calibration_df,
  aes(
    x = mean_pred,
    y = observed_rate,
    size = n
  )
) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed",
    color = "gray45"
  ) +
  geom_point(
    shape = 21,
    fill = "#C27D38",
    color = "gray20",
    alpha = 0.85
  ) +
  geom_line(
    linewidth = 1.1,
    color = "#798E87"
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = percent_format(accuracy = 1)
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = percent_format(accuracy = 1)
  ) +
  scale_size_continuous(range = c(3, 10)) +
  labs(
    x = "Mean predicted high-risk probability",
    y = "Observed high-risk rate",
    size = "Records"
  ) +
  annotate(
    "text",
    x = 0.65,
    y = 0.15,
    label = paste0(
      "ECE = ", round(ECE, 3),
      "\nBrier = ", round(Brier_hybrid, 3)
    ),
    size = 5,
    hjust = 0
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

7.6 6.6 Policy-Triage Gain Curves

make_gain_curve <- function(prob, truth, model_name) {

  tibble(
    prob = prob,
    truth = truth
  ) %>%
    arrange(desc(prob)) %>%
    mutate(
      review_rank = row_number(),
      high_risk = truth == 1,
      cumulative_high_risk = cumsum(high_risk),
      total_high_risk = sum(high_risk),
      recall_high_risk = cumulative_high_risk / total_high_risk,
      review_share = review_rank / n(),
      Model = model_name
    )
}

gain_curves <- bind_rows(
  make_gain_curve(rule_prob, y_test, "Rule score only"),
  make_gain_curve(meta_prob, y_test, "Metadata only"),
  make_gain_curve(lsa_prob, y_test, "Text embeddings only"),
  make_gain_curve(hybrid_prob, y_test, "Hybrid policy-text model")
)



ggplot(
  gain_curves,
  aes(
    x = review_share,
    y = recall_high_risk,
    color = Model
  )
) +
  geom_line(linewidth = 1.2) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed",
    color = "gray55"
  ) +
  scale_color_manual(
    values = c(
      "Rule score only" = "#D8A499",
      "Metadata only" = "#CCC591",
      "Text embeddings only" = "#798E87",
      "Hybrid policy-text model" = "#C27D38"
    )
  ) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  labs(
    x = "Share of feedback reviewed",
    y = "Share of urgent/high-risk records captured",
    color = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "bottom"
  )

7.7 6.7 Conformal Triage

# Split training into proper training and calibration set
set.seed(2026)

cal_split <- initial_split(
  train_dat,
  prop = 0.75,
  strata = high_risk_factor
)

proper_train <- training(cal_split)
cal_dat      <- testing(cal_split)

x_proper <- model.matrix(
  high_risk_factor ~ policy_priority_score + IncidentType + dominant_risk_domain + .,
  data = proper_train %>%
    select(
      high_risk_factor,
      policy_priority_score,
      IncidentType,
      dominant_risk_domain,
      all_of(lsa_cols)
    )
)[, -1]

x_cal <- model.matrix(
  high_risk_factor ~ policy_priority_score + IncidentType + dominant_risk_domain + .,
  data = cal_dat %>%
    select(
      high_risk_factor,
      policy_priority_score,
      IncidentType,
      dominant_risk_domain,
      all_of(lsa_cols)
    )
)[, -1]

y_proper <- ifelse(proper_train$high_risk_factor == "UrgentHigh", 1, 0)
y_cal    <- ifelse(cal_dat$high_risk_factor == "UrgentHigh", 1, 0)

set.seed(2026)

cv_conf <- cv.glmnet(
  x = x_proper,
  y = y_proper,
  family = "binomial",
  alpha = 0.5,
  type.measure = "auc"
)

cal_prob <- as.numeric(
  predict(cv_conf, newx = x_cal, s = "lambda.min", type = "response")
)

# Nonconformity score for true class
cal_nonconf <- ifelse(y_cal == 1, 1 - cal_prob, cal_prob)

alpha <- 0.10
qhat <- quantile(
  cal_nonconf,
  probs = ceiling((length(cal_nonconf) + 1) * (1 - alpha)) / length(cal_nonconf),
  na.rm = TRUE
)

# Apply to test set
conf_prob_test <- as.numeric(
  predict(cv_conf, newx = x_test_hybrid, s = "lambda.min", type = "response")
)

conformal_test <- test_dat %>%
  mutate(
    conformal_prob = conf_prob_test,
    include_lowmoderate = conformal_prob <= qhat,
    include_urgenthigh = (1 - conformal_prob) <= qhat,
    conformal_set = case_when(
      include_lowmoderate & include_urgenthigh ~ "{Low/Moderate, Urgent/High}",
      include_urgenthigh ~ "{Urgent/High}",
      include_lowmoderate ~ "{Low/Moderate}",
      TRUE ~ "{Review required}"
    )
  )

7.8 6.8 Policy Coverage Simulation

policy_sim <- dat_policy %>%
  mutate(
    high_risk = criticality_class %in% c("Urgent", "High")
  ) %>%
  group_by(dominant_risk_domain, recommended_policy_action) %>%
  summarise(
    total_records = n(),
    urgent_high_records = sum(high_risk),
    mean_priority = mean(policy_priority_score, na.rm = TRUE),
    max_priority = max(policy_priority_score, na.rm = TRUE),
    coverage_share = 100 * urgent_high_records / sum(dat_policy$criticality_class %in% c("Urgent", "High")),
    .groups = "drop"
  ) %>%
  arrange(desc(urgent_high_records), desc(mean_priority))

policy_sim
## # A tibble: 7 × 7
##   dominant_risk_domain  recommended_policy_a…¹ total_records urgent_high_records
##   <chr>                 <chr>                          <int>               <int>
## 1 Emergency-response i… Require emergency veh…            82                  81
## 2 Collision / recovery… Require rapid removal…            60                  59
## 3 Vulnerable road-user… Apply VRU-sensitive O…            28                  23
## 4 Traffic obstruction   Require minimum-risk …            64                  10
## 5 Manual traffic-contr… Restrict AV operation…             5                   5
## 6 General AV concern    Monitor and classify …            22                   2
## 7 Residential disrupti… Apply nighttime resid…             6                   0
## # ℹ abbreviated name: ¹​recommended_policy_action
## # ℹ 3 more variables: mean_priority <dbl>, max_priority <dbl>,
## #   coverage_share <dbl>
ggplot(
  policy_sim,
  aes(
    x = urgent_high_records,
    y = fct_reorder(dominant_risk_domain, urgent_high_records),
    fill = mean_priority
  )
) +
  geom_col(width = 0.72) +
  geom_text(
    aes(
      label = paste0(
        urgent_high_records,
        " cases; ",
        round(coverage_share, 1),
        "%"
      )
    ),
    hjust = -0.05,
    size = 4.6
  ) +
  scale_fill_gradientn(
    colors = risk_cols,
    name = "Mean priority"
  ) +
  scale_x_continuous(
    expand = expansion(mult = c(0, 0.25))
  ) +
  labs(
    x = "Urgent/high-risk records covered",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.major.y = element_blank(),
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

7.9 6.9 Representative High-Priority Cases

representative_cases <- dat_policy %>%
  group_by(dominant_risk_domain) %>%
  arrange(
    desc(policy_priority_score),
    desc(model_high_risk_prob),
    .by_group = TRUE
  ) %>%
  slice_head(n = 3) %>%
  ungroup() %>%
  select(
    dominant_risk_domain,
    IncidentType,
    criticality_class,
    policy_priority_score,
    model_high_risk_prob,
    recommended_policy_action,
    Feedback
  )

representative_cases
## # A tibble: 21 × 7
##    dominant_risk_domain     IncidentType criticality_class policy_priority_score
##    <chr>                    <fct>        <chr>                             <dbl>
##  1 Collision / recovery fa… Near Miss    Urgent                               26
##  2 Collision / recovery fa… Near Miss    Urgent                               24
##  3 Collision / recovery fa… Collision    Urgent                               22
##  4 Emergency-response inte… Safety Conc… Urgent                               19
##  5 Emergency-response inte… Near Miss    Urgent                               18
##  6 Emergency-response inte… Blocking Tr… Urgent                               17
##  7 General AV concern       Collision    High                                  6
##  8 General AV concern       Collision    High                                  6
##  9 General AV concern       Nuisance     Low                                   2
## 10 Manual traffic-control … Safety Conc… Urgent                               10
## # ℹ 11 more rows
## # ℹ 3 more variables: model_high_risk_prob <dbl>,
## #   recommended_policy_action <chr>, Feedback <chr>
write_xlsx(
  representative_cases,
  "Representative_high_priority_AV_feedback_cases.xlsx"
)

7.10 6.10 Spatial Hotspot Analysis

if (all(c("Latitude", "Longitude") %in% names(dat_policy))) {

  spatial_sf <- dat_policy %>%
    filter(!is.na(Latitude), !is.na(Longitude)) %>%
    st_as_sf(
      coords = c("Longitude", "Latitude"),
      crs = 4326,
      remove = FALSE
    ) %>%
    st_transform(3857)

  # Create 500-meter grid
  grid_500m <- st_make_grid(
    spatial_sf,
    cellsize = 500,
    square = TRUE
  ) %>%
    st_as_sf() %>%
    mutate(grid_id = row_number())

  spatial_joined <- st_join(
    spatial_sf,
    grid_500m,
    join = st_within
  )

  grid_risk <- spatial_joined %>%
    st_drop_geometry() %>%
    group_by(grid_id) %>%
    summarise(
      n_reports = n(),
      mean_priority = mean(policy_priority_score, na.rm = TRUE),
      urgent_high = sum(criticality_class %in% c("Urgent", "High")),
      .groups = "drop"
    ) %>%
    right_join(grid_500m, by = "grid_id") %>%
    st_as_sf() %>%
    mutate(
      n_reports = replace_na(n_reports, 0),
      mean_priority = replace_na(mean_priority, 0),
      urgent_high = replace_na(urgent_high, 0)
    ) %>%
    filter(n_reports > 0)

  grid_risk_4326 <- grid_risk %>%
    st_transform(4326)
}


ggplot() +
  geom_sf(
    data = grid_risk_4326,
    aes(fill = mean_priority),
    color = "white",
    linewidth = 0.15,
    alpha = 0.90
  ) +
  geom_point(
    data = dat_policy %>%
      filter(!is.na(Latitude), !is.na(Longitude)),
    aes(
      x = Longitude,
      y = Latitude,
      size = policy_priority_score
    ),
    color = "gray20",
    alpha = 0.35
  ) +
  scale_fill_gradientn(
    colors = risk_cols,
    name = "Mean priority"
  ) +
  scale_size_continuous(
    range = c(1.5, 5),
    name = "Case priority"
  ) +
  labs(
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

7.11 6.11 Risk-Probability Summary and Final Export

risk_prob_summary <- dat_policy %>%
  group_by(dominant_risk_domain) %>%
  summarise(
    n = n(),
    mean_prob = mean(model_high_risk_prob, na.rm = TRUE),
    median_prob = median(model_high_risk_prob, na.rm = TRUE),
    mean_priority = mean(policy_priority_score, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(desc(mean_prob))

ggplot(
  risk_prob_summary,
  aes(
    x = mean_prob,
    y = fct_reorder(dominant_risk_domain, mean_prob),
    size = n,
    fill = mean_priority
  )
) +
  geom_point(
    shape = 21,
    color = "gray25",
    alpha = 0.90,
    stroke = 0.35
  ) +
  scale_fill_gradientn(
    colors = risk_cols,
    name = "Mean priority"
  ) +
  scale_size_continuous(
    range = c(4, 14),
    labels = comma,
    name = "Records"
  ) +
  scale_x_continuous(
    labels = percent_format(accuracy = 1),
    limits = c(0, 1)
  ) +
  labs(
    x = "Mean predicted high-risk probability",
    y = NULL
  ) +
  theme_minimal(base_size = 16) +
  theme(
    panel.grid.minor = element_blank(),
    legend.position = "right"
  )

ggsave(
  "Fig_risk_domain_probability_bubble.png",
  width = 10,
  height = 6.5,
  dpi = 400
)


advanced_outputs <- list(
  model_comparison = model_comparison,
  calibration = calibration_df,
  policy_simulation = policy_sim,
  representative_cases = representative_cases,
  conformal_sets = conformal_test %>%
    st_drop_geometry() %>%
    select(
      id,
      IncidentType,
      criticality_class,
      policy_priority_score,
      conformal_prob,
      conformal_set,
      dominant_risk_domain,
      Feedback
    ),
  risk_probability_summary = risk_prob_summary
)

write_xlsx(
  advanced_outputs,
  "CivicRisk_AV_advanced_outputs.xlsx"
)

8 7. Reproducibility Notes