pkgs <- c("sf", "dplyr", "readr", "readxl", "stringr", "MASS", "sandwich", "lmtest")
missing <- pkgs[!pkgs %in% rownames(installed.packages())]
if (length(missing) > 0) install.packages(missing, repos = "https://cloud.r-project.org")
library(sf) # geometry, distances, overlay
## Warning: package 'sf' was built under R version 4.5.3
## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE
library(dplyr) # data wrangling
## Warning: package 'dplyr' was built under R version 4.5.3
##
## 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
library(readr) # CSV reading
## Warning: package 'readr' was built under R version 4.5.3
library(readxl) # Excel reading
## Warning: package 'readxl' was built under R version 4.5.3
library(stringr) # street-name cleaning
## Warning: package 'stringr' was built under R version 4.5.3
library(MASS) # glm.nb (negative binomial)
## Warning: package 'MASS' was built under R version 4.5.3
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
library(sandwich) # cluster-robust standard errors
## Warning: package 'sandwich' was built under R version 4.5.3
library(lmtest) # coeftest
## Warning: package 'lmtest' was built under R version 4.5.3
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
options(width = 130)
data_dir <- "C:/Users/jalex/OneDrive/Documents/JS PAD Captsone Files/excel data from COSA"
if (!dir.exists(data_dir)) stop("data_dir does not exist: ", data_dir)
out_dir <- file.path(data_dir, "R_outputs")
dir.create(out_dir, showWarnings = FALSE)
list.files(data_dir)
## [1] "CorridorPlans_7237623346671139028.csv" "Corridors_-5661148773361461175.csv"
## [3] "Major_Thoroughfare_Plan__MTP_3253228272832874683.csv" "my_table.csv"
## [5] "pwSevereBicyclistInjuryAreas.geojson" "pwSevereBicyclistInjuryAreas_8521414027376615119.csv"
## [7] "pwSeverePedestrianInjuryAreas.geojson" "pwSeverePedestrianInjuryAreas_7442757312003544754.csv"
## [9] "pwVisionZeroProjects_1393929762377528359 (1).csv" "R_outputs"
## [11] "Streets_649949320172161544.csv"
files <- c(
pedestrian = "pwSeverePedestrianInjuryAreas_7442757312003544754",
bicyclist = "pwSevereBicyclistInjuryAreas_8521414027376615119",
projects = "pwVisionZeroProjects_1393929762377528359 (1)",
streets = "Streets_649949320172161544",
mtp = "Major_Thoroughfare_Plan__MTP_3253228272832874683",
corridors = "Corridors_-5661148773361461175",
plans = "CorridorPlans_7237623346671139028",
cris = "my_table"
)
city_path <- function(stem) {
f <- file.path(data_dir, paste0(stem, ".xlsx"))
if (!file.exists(f)) f <- file.path(data_dir, paste0(stem, ".csv"))
if (!file.exists(f)) stop("Cannot find ", stem, " (.xlsx or .csv) in ", data_dir)
f
}
read_city <- function(stem) {
f <- city_path(stem)
if (tolower(tools::file_ext(f)) == "xlsx") {
d <- read_excel(f, col_types = "text", .name_repair = "minimal")
} else {
d <- read_csv(f, col_types = cols(.default = col_character()), show_col_types = FALSE, progress = FALSE)
}
names(d) <- str_squish(str_replace_all(names(d), "\\s+", " "))
d %>% mutate(across(everything(), ~ str_squish(as.character(.x))))
}
num <- function(x) suppressWarnings(as.numeric(x))
spia_csv <- read_city(files["pedestrian"])
sbia_csv <- read_city(files["bicyclist"])
proj_csv <- read_city(files["projects"])
mtp_csv <- read_city(files["mtp"])
streets <- read_city(files["streets"])
corr_csv <- read_city(files["corridors"])
plan_csv <- read_city(files["plans"])
c(pedestrian = nrow(spia_csv), bicyclist = nrow(sbia_csv), projects = nrow(proj_csv),
streets = nrow(streets), mtp = nrow(mtp_csv), corridors = nrow(corr_csv), plans = nrow(plan_csv))
## pedestrian bicyclist projects streets mtp corridors plans
## 166 14 365 92216 1405 20 12
# expect: 166, 14, 365, 92216, 1405, 20, 12
spia_csv <- spia_csv %>% mutate(across(c(OBJECTID, CorridorID, Incapacitated_Injuries, Fatal_Injuries, Total_Injuries, SPIA_Year), ~ as.integer(num(.x))), Shape__Length = num(Shape__Length))
sbia_csv <- sbia_csv %>% mutate(across(c(OBJECTID, Incapacitated_Injuries, Fatal_Injuries, Total_Injuries, SBIA_Year), ~ as.integer(num(.x))), Shape__Length = num(Shape__Length))
all(spia_csv$Total_Injuries == spia_csv$Fatal_Injuries + spia_csv$Incapacitated_Injuries)
## [1] TRUE
table(spia_csv$SPIA_Year)
##
## 2017 2020
## 76 90
spia_csv %>% group_by(SPIA_Year) %>%
summarise(segments = n(), miles = sum(Shape__Length) / 5280, killed = sum(Fatal_Injuries), incap = sum(Incapacitated_Injuries), ksi = sum(Total_Injuries))
## # A tibble: 2 × 6
## SPIA_Year segments miles killed incap ksi
## <int> <int> <dbl> <int> <int> <int>
## 1 2017 76 37.3 98 185 283
## 2 2020 90 53.4 115 252 367
c(bike_miles = sum(sbia_csv$Shape__Length) / 5280, bike_killed = sum(sbia_csv$Fatal_Injuries), bike_incap = sum(sbia_csv$Incapacitated_Injuries))
## bike_miles bike_killed bike_incap
## 12.82048 4.00000 31.00000
geo_url <- c(pedestrian = "https://services.arcgis.com/g1fRTDLeMgspWrYp/ArcGIS/rest/services/pwSeverePedestrianInjuryAreas/FeatureServer/0/query?where=1%3D1&outFields=*&f=geojson",
bicyclist = "https://services.arcgis.com/g1fRTDLeMgspWrYp/ArcGIS/rest/services/pwSevereBicyclistInjuryAreas/FeatureServer/0/query?where=1%3D1&outFields=*&f=geojson")
read_geom <- function(pattern, what, save_as) {
f <- list.files(data_dir, pattern = pattern, full.names = TRUE, ignore.case = TRUE)
if (length(f) == 0) {
dest <- file.path(data_dir, save_as)
message("No GeoJSON for the ", what, " layer in the data folder; downloading it to ", basename(dest))
ok <- suppressWarnings(tryCatch({ utils::download.file(geo_url[[what]], dest, mode = "wb", quiet = TRUE); TRUE }, error = function(e) FALSE))
if (!ok || !file.exists(dest) || file.size(dest) < 1000)
stop("Could not download the ", what, " layer. Copy spia_2278.geojson and sbia_2278.geojson into ", data_dir, " and run this step again.")
f <- dest
}
g <- st_read(f[1], quiet = TRUE) %>% mutate(across(-geometry, as.character))
names(g) <- str_replace_all(names(g), "\\s+", "")
bb <- st_bbox(g)
if (abs(bb["xmin"]) > 360 || abs(bb["ymin"]) > 360) {
suppressWarnings(st_crs(g) <- 2278)
} else {
if (is.na(st_crs(g))) st_crs(g) <- 4326
g <- st_transform(g, 2278)
}
message("Read ", basename(f[1]), ": ", nrow(g), " features")
g
}
spia_geo <- read_geom("(Pedestrian|spia).*\\.geojson$", "pedestrian", "pwSeverePedestrianInjuryAreas.geojson")
## Read pwSeverePedestrianInjuryAreas.geojson: 166 features
sbia_geo <- read_geom("(Bicyclist|sbia).*\\.geojson$", "bicyclist", "pwSevereBicyclistInjuryAreas.geojson")
## Read pwSevereBicyclistInjuryAreas.geojson: 14 features
spia_csv <- spia_csv %>% mutate(across(c(OBJECTID, CorridorID, Incapacitated_Injuries, Fatal_Injuries, Total_Injuries, SPIA_Year), ~ as.integer(num(.x))), Shape__Length = num(Shape__Length))
sbia_csv <- sbia_csv %>% mutate(across(c(OBJECTID, Incapacitated_Injuries, Fatal_Injuries, Total_Injuries, SBIA_Year), ~ as.integer(num(.x))), Shape__Length = num(Shape__Length))
spia_geo <- spia_geo %>% mutate(CorridorID = as.integer(num(CorridorID))) %>% dplyr::select(CorridorID, geometry)
spia <- spia_csv %>% inner_join(spia_geo, by = "CorridorID") %>% st_as_sf()
nrow(spia)
## [1] 166
sbia_geo <- sbia_geo %>% mutate(FROM = toupper(FromStreet), TO = toupper(ToStreet), glen = as.numeric(st_length(geometry)))
idx <- sapply(seq_len(nrow(sbia_csv)), function(i) {
cand <- which(sbia_geo$FROM == toupper(sbia_csv$FromStreet[i]) & sbia_geo$TO == toupper(sbia_csv$ToStreet[i]) & abs(sbia_geo$glen - sbia_csv$Shape__Length[i]) < 1)
if (length(cand) == 1) cand else NA_integer_ })
sum(!is.na(idx))
## [1] 14
sbia <- st_sf(sbia_csv, geometry = st_geometry(sbia_geo)[idx])
cor <- bind_rows(spia %>% mutate(uid = paste0("P", CorridorID), mode = "Pedestrian", year = SPIA_Year),
sbia %>% mutate(uid = paste0("B", OBJECTID), mode = "Bicyclist", year = SBIA_Year)) %>%
mutate(fatal = Fatal_Injuries, incap = Incapacitated_Injuries, ksi = Total_Injuries,
length_ft = Shape__Length, miles = length_ft / 5280, computed_ft = as.numeric(st_length(geometry)))
max(abs(cor$computed_ft - cor$length_ft))
## [1] 0.002647927
table(cor$mode, cor$year)
##
## 2017 2020
## Bicyclist 0 14
## Pedestrian 76 90
pts <- proj_csv %>% mutate(x = num(x), y = num(y), numbered = str_detect(ProjectNumber, "^VZ_\\d+")) %>%
st_as_sf(coords = c("x", "y"), crs = 2278, remove = FALSE)
c(points = nrow(pts), locations = nrow(distinct(st_drop_geometry(pts), x, y)), numbered = sum(pts$numbered))
## points locations numbered
## 365 193 108
dmat <- st_distance(pts, cor)
nearest <- as.numeric(apply(dmat, 1, min))
round(quantile(nearest, c(.1, .25, .5, .75, .9)))
## 10% 25% 50% 75% 90%
## 2 9 1281 2889 4681
for (b in c(100, 150, 250, 500, 1000)) cat(b, "ft:", sum(nearest <= b), "points (", round(100 * mean(nearest <= b), 1), "%)\n")
## 100 ft: 130 points ( 35.6 %)
## 150 ft: 130 points ( 35.6 %)
## 250 ft: 136 points ( 37.3 %)
## 500 ft: 147 points ( 40.3 %)
## 1000 ft: 165 points ( 45.2 %)
bins <- c(0, 10, 50, 150, 500, 1000, 2500, 5000, Inf)
table(cut(nearest, bins, right = TRUE, include.lowest = TRUE))
##
## [0,10] (10,50] (50,150] (150,500] (500,1e+03] (1e+03,2.5e+03] (2.5e+03,5e+03] (5e+03,Inf]
## 101 24 5 17 18 81 87 32
for (b in c(100, 150, 250)) {
within <- st_is_within_distance(cor, pts, dist = b)
cor[[paste0("pts_", b)]] <- lengths(within)
cor[[paste0("ptsnum_", b)]] <- sapply(within, function(i) sum(pts$numbered[i]))
cor[[paste0("treated_", b)]] <- as.integer(lengths(within) > 0)
cor[[paste0("dose_", b)]] <- lengths(within) / cor$miles
}
cor <- cor %>% mutate(treated = treated_150, dose = dose_150, later = as.integer(year == 2020), ksi_per_mi = ksi / miles, logmiles = log(miles))
cor %>% st_drop_geometry() %>% group_by(mode, year) %>% summarise(segments = n(), treated = sum(treated), points = sum(pts_150), .groups = "drop")
## # A tibble: 3 × 5
## mode year segments treated points
## <chr> <int> <int> <int> <int>
## 1 Bicyclist 2020 14 4 14
## 2 Pedestrian 2017 76 23 98
## 3 Pedestrian 2020 90 20 112
write.csv(st_drop_geometry(pts) %>% mutate(nearest_ft = nearest), file.path(out_dir, "project_points_nearest_corridor.csv"), row.names = FALSE)
street_key <- function(name) {
s <- toupper(name)
s <- str_replace_all(s, "[.']", "")
s <- str_replace_all(s, "\\b(SPUR|FM|LOOP|US|IH|SH|RM)\\s*\\d+\\b", "")
s <- str_replace_all(s, "\\bHIGHWAY\\b|\\bHWY\\b", "")
str_squish(s)
}
cor <- cor %>% mutate(street_key = street_key(StreetName))
ped_keys <- st_drop_geometry(cor) %>% filter(mode == "Pedestrian")
length(unique(ped_keys$street_key))
## [1] 70
both <- intersect(ped_keys$street_key[ped_keys$year == 2017], ped_keys$street_key[ped_keys$year == 2020]); length(both) #
## [1] 35
streets <- streets %>% mutate(MSAG_NAME = toupper(MSAG_NAME), LengthFeet = num(LengthFeet),
Owner = toupper(coalesce(Owner, "UNKNOWN")), CouncilDistrict = coalesce(CouncilDistrict, "Unknown"))
mtp <- mtp_csv %>% mutate(StreetName_u = toupper(StreetName), len = num(Shape__Length))
word_match <- function(x, key) str_detect(x, paste0("(^|\\s)", str_replace_all(key, "([\\W])", "\\\\\\1"), "(\\s|$)"))
ctrl <- lapply(sort(unique(cor$street_key)), function(key) {
s <- streets[word_match(streets$MSAG_NAME, key), ]
m <- mtp[word_match(mtp$StreetName_u, key), ]
by_d <- s %>% filter(CouncilDistrict != "0", CouncilDistrict != "Unknown") %>% group_by(CouncilDistrict) %>% summarise(L = sum(LengthFeet, na.rm = TRUE))
mt <- m %>% group_by(Type) %>% summarise(L = sum(len, na.rm = TRUE))
tibble(street_key = key, streets_rows = nrow(s),
txdot_share = if (nrow(s)) sum(s$LengthFeet[s$Owner == "TXDOT"], na.rm = TRUE) / sum(s$LengthFeet, na.rm = TRUE) else NA_real_,
council_district = if (nrow(by_d)) by_d$CouncilDistrict[which.max(by_d$L)] else "Unmatched",
mtp_type = if (nrow(mt)) mt$Type[which.max(mt$L)] else NA_character_)
}) %>% bind_rows() %>%
mutate(state_owned = ifelse(is.na(txdot_share), NA, as.integer(txdot_share >= 0.5)),
arterial_class = case_when(is.na(mtp_type) ~ NA_character_,
str_detect(toupper(mtp_type), "PRIMARY|SUPER|FREEWAY|EXPRESSWAY") ~ "Primary/super",
str_detect(toupper(mtp_type), "SECONDARY|ENHANCED") ~ "Secondary",
TRUE ~ "Other"))
cor <- cor %>% left_join(ctrl, by = "street_key")
c(streets_coverage = mean(cor$streets_rows > 0), mtp_coverage = mean(!is.na(cor$arterial_class)))
## streets_coverage mtp_coverage
## 0.8555556 0.8111111
table(cor$arterial_class, useNA = "ifany"); table(cor$state_owned, useNA = "ifany")
##
## Primary/super Secondary <NA>
## 62 84 34
##
## 0 1 <NA>
## 111 43 26
dat <- st_drop_geometry(cor)
ped <- dat %>% filter(mode == "Pedestrian")
bike <- dat %>% filter(mode == "Bicyclist")
write.csv(dat, file.path(out_dir, "corridors_analysis_table.csv"), row.names = FALSE)
saveRDS(dat, file.path(out_dir, "corridors_analysis_table.rds"))
st_write(cor, file.path(out_dir, "corridors_with_treatment.geojson"), delete_dsn = TRUE, quiet = TRUE)
desc <- function(d) d %>% summarise(segments = n(), median_length = median(miles), at_floor = sum(ksi == 2),
seg_median_rate = median(ksi_per_mi), seg_mean_rate = mean(ksi_per_mi), seg_sd_rate = sd(ksi_per_mi),
treated = sum(treated), points = sum(pts_150), killed = sum(fatal), incap = sum(incap),
ksi_per_mile = sum(ksi) / sum(miles), miles = sum(miles), ksi = sum(ksi))
table1 <- bind_rows(desc(filter(ped, year == 2017)) %>% mutate(group = "Pedestrian 2017"),
desc(filter(ped, year == 2020)) %>% mutate(group = "Pedestrian 2020"),
desc(ped) %>% mutate(group = "Pedestrian both"),
desc(bike) %>% mutate(group = "Bicyclist 2020")) %>% dplyr::select(group, everything())
print(as.data.frame(table1 %>% mutate(across(where(is.numeric), ~ round(.x, 2)))))
## group segments median_length at_floor seg_median_rate seg_mean_rate seg_sd_rate treated points killed incap
## 1 Pedestrian 2017 76 0.34 32 8.47 15.53 24.54 23 98 98 185
## 2 Pedestrian 2020 90 0.42 41 6.91 21.19 52.88 20 112 115 252
## 3 Pedestrian both 166 0.37 73 7.46 18.60 42.31 43 210 213 437
## 4 Bicyclist 2020 14 0.76 9 3.28 12.84 27.33 4 14 4 31
## ksi_per_mile miles ksi
## 1 7.59 37.30 283
## 2 6.88 53.37 367
## 3 7.17 90.67 650
## 4 2.73 12.82 35
c(mean = mean(ped$ksi), variance = var(ped$ksi))
## mean variance
## 3.915663 7.883753
write.csv(table1, file.path(out_dir, "table1_descriptives.csv"), row.names = FALSE)
cells <- ped %>% group_by(year, treated) %>%
summarise(segments = n(), miles = sum(miles), killed = sum(fatal), ksi = sum(ksi),
ksi_per_mile = sum(ksi) / sum(miles), deaths_per_mile = sum(fatal) / sum(miles), .groups = "drop")
print(as.data.frame(cells %>% mutate(across(c(miles, ksi_per_mile, deaths_per_mile), ~ round(.x, 2)))))
## year treated segments miles killed ksi ksi_per_mile deaths_per_mile
## 1 2017 0 53 19.21 52 160 8.33 2.71
## 2 2017 1 23 18.09 46 123 6.80 2.54
## 3 2020 0 70 33.41 73 242 7.24 2.18
## 4 2020 1 20 19.96 42 125 6.26 2.10
r <- function(y, t) cells$ksi_per_mile[cells$year == y & cells$treated == t]
raw_did <- (r(2020, 1) / r(2017, 1)) / (r(2020, 0) / r(2017, 0)); raw_did
## [1] 1.059774
write.csv(cells, file.path(out_dir, "table2_cells.csv"), row.names = FALSE)
panel <- ped %>% filter(street_key %in% both) %>% group_by(street_key, year) %>%
summarise(ksi = sum(ksi), miles = sum(miles), treated = max(treated), .groups = "drop") %>%
group_by(street_key) %>% mutate(treated_street = max(treated)) %>% ungroup() %>%
mutate(later = as.integer(year == 2020), logmiles = log(miles))
pooled <- panel %>% group_by(year) %>% summarise(ksi = sum(ksi), miles = sum(miles), rate = sum(ksi) / sum(miles)); pooled # expect: 7.25 -> 6.99
## # A tibble: 2 × 4
## year ksi miles rate
## <int> <int> <dbl> <dbl>
## 1 2017 231 31.9 7.25
## 2 2020 296 42.3 6.99
rr_p <- pooled$rate[2] / pooled$rate[1]; se_p <- sqrt(1 / pooled$ksi[1] + 1 / pooled$ksi[2])
round(c(ratio = rr_p, lower = exp(log(rr_p) - 1.96 * se_p), upper = exp(log(rr_p) + 1.96 * se_p)), 3) # expect: 0.965 (0.812, 1.146)
## ratio lower upper
## 0.965 0.812 1.146
panel %>% group_by(treated_street, year) %>% summarise(streets = n(), rate = round(sum(ksi) / sum(miles), 2), .groups = "drop")
## # A tibble: 4 × 4
## treated_street year streets rate
## <int> <int> <int> <dbl>
## 1 0 2017 18 8.24
## 2 0 2020 18 6.66
## 3 1 2017 17 6.89
## 4 1 2020 17 7.15
# expect: untreated 18 streets 8.24 -> 6.66 ; treated 17 streets 6.89 -> 7.15
write.csv(panel, file.path(out_dir, "street_panel_35_streets.csv"), row.names = FALSE)
rr_table <- function(fit, cluster, terms) {
V <- vcovCL(fit, cluster = cluster, type = "HC1"); ct <- coeftest(fit, vcov. = V)
data.frame(term = terms, rate_ratio = round(exp(ct[terms, "Estimate"]), 3),
lower = round(exp(ct[terms, "Estimate"] - 1.96 * ct[terms, "Std. Error"]), 3),
upper = round(exp(ct[terms, "Estimate"] + 1.96 * ct[terms, "Std. Error"]), 3),
p = round(ct[terms, "Pr(>|z|)"], 3), row.names = NULL)
}
nb <- glm.nb(ksi ~ treated * later + offset(logmiles), data = ped)
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in theta.ml(Y, mu, sum(w), w, limit = control$maxit, trace = control$trace > : iteration limit reached
## Warning in glm.nb(ksi ~ treated * later + offset(logmiles), data = ped): alternation limit reached
nb$theta
## [1] 1572.952
rr_table(nb, ped$street_key, c("treated", "later", "treated:later"))
## term rate_ratio lower upper p
## 1 treated 0.816 0.658 1.013 0.065
## 2 later 0.869 0.727 1.039 0.124
## 3 treated:later 1.060 0.819 1.371 0.660
m1 <- glm(ksi ~ treated * later + offset(logmiles), family = poisson, data = ped)
sum(residuals(m1, type = "pearson")^2) / m1$df.residual
## [1] 3.399414
table3 <- rr_table(m1, ped$street_key, c("treated", "later", "treated:later")); table3
## term rate_ratio lower upper p
## 1 treated 0.816 0.657 1.013 0.065
## 2 later 0.869 0.727 1.040 0.125
## 3 treated:later 1.060 0.819 1.372 0.659
se_did <- sqrt(vcovCL(m1, cluster = ped$street_key, type = "HC1")["treated:later", "treated:later"])
mde80 <- exp(-(1.96 + 0.84) * se_did); mde80
## [1] 0.6917744
write.csv(table3, file.path(out_dir, "table3_primary_model.csv"), row.names = FALSE)
m_len <- glm(ksi ~ treated * later + logmiles, family = poisson, data = ped)
table3_len <- rr_table(m_len, ped$street_key, c("treated", "later", "treated:later", "logmiles")); table3_len
## term rate_ratio lower upper p
## 1 treated 1.263 1.028 1.552 0.026
## 2 later 1.037 0.909 1.183 0.590
## 3 treated:later 0.991 0.827 1.188 0.921
## 4 logmiles 1.589 1.398 1.805 0.000
round(coef(m_len)["logmiles"], 3)
## logmiles
## 0.463
write.csv(table3_len, file.path(out_dir, "table3_length_as_covariate.csv"), row.names = FALSE)
checks <- list()
m_fe <- glm(ksi ~ treated * later + factor(street_key) + offset(logmiles), family = poisson, data = ped)
checks$fixed_effects <- rr_table(m_fe, ped$street_key, "treated:later")
checks$length_covariate <- table3_len[3, ]
pc <- ped %>% filter(!is.na(arterial_class), !is.na(state_owned)) %>% mutate(primary = as.integer(arterial_class == "Primary/super")); nrow(pc)
## [1] 146
m_ctl <- glm(ksi ~ treated * later + primary + state_owned + offset(logmiles), family = poisson, data = pc)
checks$controls <- rr_table(m_ctl, pc$street_key, "treated:later")
m_pan <- glm(ksi ~ treated_street * later + offset(logmiles), family = poisson, data = panel)
checks$street_panel <- rr_table(m_pan, panel$street_key, "treated_street:later")
for (b in c(100, 250)) { d <- ped; d$tr <- d[[paste0("treated_", b)]]
checks[[paste0("band_", b)]] <- rr_table(glm(ksi ~ tr * later + offset(logmiles), family = poisson, data = d), d$street_key, "tr:later") }
d <- ped %>% mutate(tr = as.integer(ptsnum_150 > 0)); sum(d$tr)
## [1] 19
checks$numbered_only <- rr_table(glm(ksi ~ tr * later + offset(logmiles), family = poisson, data = d), d$street_key, "tr:later")
spf <- glm(ksi ~ later + offset(logmiles), family = poisson, data = filter(ped, treated == 0))
eb <- bind_rows(lapply(c(2017, 2020), function(y) { tr <- filter(ped, treated == 1, year == y); e <- sum(predict(spf, newdata = tr, type = "response"))
data.frame(year = y, segments = nrow(tr), observed = sum(tr$ksi), expected = round(e, 1), ratio = round(sum(tr$ksi) / e, 3)) })); eb
## year segments observed expected ratio
## 1 2017 23 123 150.7 0.816
## 2 2020 20 125 144.5 0.865
prio <- c("FREDERICKSBURG", "WEST", "NEW BRAUNFELS", "COMMERCE", "GOLIAD", "MARBACH", "FIVE PALMS", "GENERAL MCMULLEN",
"MILITARY", "GRISSOM", "CULEBRA", "CALLAGHAN", "WURZBACH", "BABCOCK", "NACOGDOCHES", "THOUSAND OAKS")
unt <- ped %>% filter(treated == 0) %>% mutate(placebo = as.integer(street_key %in% prio)); sum(unt$placebo)
## [1] 34
checks$placebo <- rr_table(glm(ksi ~ placebo * later + offset(logmiles), family = poisson, data = unt), unt$street_key, "placebo:later")
set.seed(20261105)
streets_u <- unique(ped$street_key); groups <- split(ped, ped$street_key); boot <- numeric(0)
for (i in 1:1000) { samp <- sample(streets_u, length(streets_u), replace = TRUE)
d <- bind_rows(lapply(seq_along(samp), function(k) groups[[samp[k]]]))
if (sum(d$treated * d$later) == 0) next
boot <- c(boot, coef(glm(ksi ~ treated * later + offset(logmiles), family = poisson, data = d))["treated:later"]) }
round(exp(c(median = median(boot), quantile(boot, c(.025, .975)))), 3); mean(boot < 0)
## median 2.5% 97.5%
## 1.068 0.841 1.429
## [1] 0.314
checks$bootstrap <- data.frame(term = "treated:later", rate_ratio = round(exp(median(boot)), 3), lower = round(exp(unname(quantile(boot, .025))), 3), upper = round(exp(unname(quantile(boot, .975))), 3), p = NA)
robust <- bind_rows(lapply(names(checks), function(n) checks[[n]] %>% mutate(check = n))) %>% dplyr::select(check, everything()) %>% as.data.frame(); rownames(robust) <- NULL; robust
## check term rate_ratio lower upper p
## 1 fixed_effects treated:later 0.923 0.655 1.300 0.646
## 2 length_covariate treated:later 0.991 0.827 1.188 0.921
## 3 controls treated:later 1.048 0.789 1.394 0.745
## 4 street_panel treated_street:later 1.283 0.954 1.726 0.100
## 5 band_100 tr:later 1.060 0.819 1.372 0.659
## 6 band_250 tr:later 1.133 0.905 1.419 0.276
## 7 numbered_only tr:later 0.887 0.697 1.128 0.327
## 8 placebo placebo:later 1.088 0.717 1.650 0.692
## 9 bootstrap treated:later 1.068 0.841 1.429 NA
write.csv(robust, file.path(out_dir, "table_robustness_checks.csv"), row.names = FALSE)
write.csv(eb, file.path(out_dir, "table_empirical_bayes_check.csv"), row.names = FALSE)
later <- dat %>% filter(year == 2020) %>% mutate(bike = as.integer(mode == "Bicyclist")); nrow(later)
## [1] 104
t4 <- list()
t4$dose <- rr_table(glm(ksi ~ dose + bike + offset(logmiles), family = poisson, data = later), later$street_key, "dose")
round(c(summary(later$dose[later$treated == 1])[c("Median", "1st Qu.", "3rd Qu.")], max = max(later$dose)), 2)
## Median 1st Qu. 3rd Qu. max
## 5.66 3.03 8.25 20.13
lc <- later %>% filter(!is.na(arterial_class), !is.na(state_owned)) %>% mutate(primary = as.integer(arterial_class == "Primary/super")); nrow(lc)
## [1] 79
t4$dose_controls <- rr_table(glm(ksi ~ dose + primary + state_owned + offset(logmiles), family = poisson, data = lc), lc$street_key, c("dose", "primary", "state_owned"))
t4$cross_section <- rr_table(glm(ksi ~ treated + primary + state_owned + offset(logmiles), family = poisson, data = lc), lc$street_key, "treated")
table4 <- bind_rows(lapply(names(t4), function(n) t4[[n]] %>% mutate(test = n))) %>% dplyr::select(test, everything()); table4
## test term rate_ratio lower upper p
## 1 dose dose 1.006 0.979 1.034 0.646
## 2 dose_controls dose 1.016 0.988 1.045 0.260
## 3 dose_controls primary 0.822 0.651 1.038 0.100
## 4 dose_controls state_owned 1.075 0.835 1.384 0.574
## 5 cross_section treated 0.907 0.767 1.071 0.250
bike_tab <- bike %>% group_by(treated) %>% summarise(segments = n(), miles = round(sum(miles), 2), ksi = sum(ksi), ksi_per_mile = round(sum(ksi) / sum(miles), 2)); bike_tab
## # A tibble: 2 × 5
## treated segments miles ksi ksi_per_mile
## <int> <int> <dbl> <int> <dbl>
## 1 0 10 4.88 22 4.51
## 2 1 4 7.94 13 1.64
write.csv(table4, file.path(out_dir, "table4_dose_response_cross_section.csv"), row.names = FALSE)
write.csv(bike_tab, file.path(out_dir, "table4_bicyclist_corridors.csv"), row.names = FALSE)
cris_path <- city_path(files["cris"])
if (tolower(tools::file_ext(cris_path)) == "xlsx") {
raw <- read_excel(cris_path, col_names = FALSE, col_types = "text", .name_repair = "minimal")
hdr <- which(str_detect(coalesce(as.character(raw[[1]]), ""), "^Contributing Factors"))[1]
cris <- as.data.frame(raw[(hdr + 1):nrow(raw), ]); names(cris) <- as.character(unlist(raw[hdr, ]))
cris <- cris %>% mutate(across(-1, num))
} else {
lines <- readLines(cris_path, encoding = "UTF-8", warn = FALSE)
hdr <- grep("^\"?Contributing Factors", lines)[1]
cris <- read.csv(text = paste(lines[hdr:length(lines)], collapse = "\n"), check.names = FALSE, stringsAsFactors = FALSE)
}
names(cris)[1] <- "factors"; cris <- cris %>% filter(!is.na(factors), factors != "Total")
nrow(cris)
## [1] 742
grp <- function(f) { parts <- str_trim(str_split(toupper(f), ";")[[1]]); non <- parts[!str_detect(parts, "SPEED")]
if (length(non) == 0) return("Speed factors only"); j <- paste(non, collapse = " ")
if (str_detect(j, "INFLUENCE|ALCOHOL|DRINKING|DRUG")) return("Speed with impairment")
if (str_detect(j, "PEDESTRIAN")) return("Speed with a pedestrian factor")
if (str_detect(j, "INATTENTION|DISTRACT|CELL")) return("Speed with inattention or distraction")
if (str_detect(j, "DISREGARD|FAILED TO YIELD|STOP SIGN|SIGNAL")) return("Speed with disregard of control or failure to yield")
"Speed with other factors" }
cris$group <- sapply(cris$factors, grp)
table6 <- cris %>% group_by(group) %>%
summarise(crashes = sum(Total), fatal = sum(`K - FATAL INJURY`), serious = sum(`A - SUSPECTED SERIOUS INJURY`), .groups = "drop") %>%
mutate(fatal_or_serious = fatal + serious, pct = round(100 * fatal_or_serious / crashes, 1)) %>% arrange(desc(crashes))
print(as.data.frame(table6)); colSums(table6[, c("crashes", "fatal", "serious")])
## group crashes fatal serious fatal_or_serious pct
## 1 Speed with inattention or distraction 3652 13 58 71 1.9
## 2 Speed factors only 3563 13 40 53 1.5
## 3 Speed with other factors 1844 13 35 48 2.6
## 4 Speed with disregard of control or failure to yield 1047 12 44 56 5.3
## 5 Speed with impairment 513 19 27 46 9.0
## 6 Speed with a pedestrian factor 16 1 5 6 37.5
## crashes fatal serious
## 10635 71 209
write.csv(table6, file.path(out_dir, "table6_cris_groups.csv"), row.names = FALSE)
table7 <- ped %>% group_by(council_district) %>%
summarise(segments = n(), miles = round(sum(miles), 1), ksi = sum(ksi), killed = sum(fatal), ksi_per_mile = round(sum(ksi) / sum(miles), 1),
treated = sum(treated), treated_pct = round(100 * sum(treated) / n()), points = sum(pts_150), .groups = "drop")
print(as.data.frame(table7))
## council_district segments miles ksi killed ksi_per_mile treated treated_pct points
## 1 1 11 4.4 39 13 8.9 1 9 3
## 2 10 12 5.5 41 12 7.5 1 8 2
## 3 2 20 9.3 72 26 7.7 2 10 6
## 4 2_3 1 0.6 3 0 5.0 0 0 0
## 5 3 28 15.6 93 32 6.0 9 32 37
## 6 4 20 9.8 78 23 8.0 8 40 46
## 7 5 18 9.9 69 19 7.0 7 39 29
## 8 6 9 8.7 62 28 7.1 6 67 42
## 9 7 7 4.9 28 5 5.7 0 0 0
## 10 8 15 9.4 66 15 7.0 3 20 6
## 11 9 13 5.7 48 21 8.4 3 23 25
## 12 Unmatched 12 6.7 51 19 7.6 3 25 14
write.csv(table7, file.path(out_dir, "table7_by_district.csv"), row.names = FALSE)
orange <- "#D95F02"; blue <- "#1F5FAD"
png(file.path(out_dir, "Figure1_project_point_distances.png"), width = 1950, height = 1000, res = 300, family = "serif")
par(mar = c(4, 4.2, 0.8, 0.5))
cnt <- table(cut(nearest, bins, right = TRUE, include.lowest = TRUE))
bp <- barplot(cnt, col = c(rep(orange, 3), rep(blue, 5)), names.arg = c("\u226410", "10-50", "50-150", "150-500", "500-1k", "1k-2.5k", "2.5k-5k", ">5k"),
ylab = "Vision Zero project points", xlab = "Distance to the nearest designated corridor (feet)", border = NA, cex.names = 0.8, ylim = c(0, 135), las = 1)
text(bp, cnt, labels = cnt, pos = 3, cex = 0.8)
legend("top", c("Within the 150-ft band", "Outside the band"), fill = c(orange, blue), border = NA, bty = "n", horiz = TRUE, cex = 0.8)
invisible(dev.off())
png(file.path(out_dir, "Figure2_ksi_per_mile_slope.png"), width = 1950, height = 1100, res = 300, family = "serif")
par(mar = c(3.2, 4.2, 0.8, 0.5))
plot(c(1, 2), c(r(2017, 0), r(2020, 0)), type = "b", pch = 19, col = blue, ylim = c(0, 10), xlim = c(0.7, 2.3), xaxt = "n", lwd = 2, las = 1,
xlab = "", ylab = "KSI per corridor mile", main = "")
lines(c(1, 2), c(r(2017, 1), r(2020, 1)), type = "b", pch = 19, col = orange, lwd = 2)
text(c(1, 2), c(r(2017, 0), r(2020, 0)), sprintf("%.1f", c(r(2017, 0), r(2020, 0))), pos = 3, cex = 0.8)
text(c(1, 2), c(r(2017, 1), r(2020, 1)), sprintf("%.1f", c(r(2017, 1), r(2020, 1))), pos = 1, cex = 0.8)
axis(1, at = c(1, 2), labels = c("2017 designation (crashes 2011-2015)", "2020 designation (crashes 2014-2018)"), cex.axis = 0.8)
legend("bottomleft", c("Corridors without a Vision Zero project", "Corridors with a Vision Zero project"), col = c(blue, orange), lwd = 2, pch = 19, bty = "n", cex = 0.8)
invisible(dev.off())
png(file.path(out_dir, "Figure3_rate_ratio_forest.png"), width = 1950, height = 1400, res = 300, family = "serif")
labels_fp <- c(primary = "Primary model (negative binomial, 150-ft band, clustered by street)", poisson = "Poisson with clustered standard errors",
fixed_effects = "Street fixed effects", length_covariate = "Segment length as a covariate instead of an offset",
controls = "With arterial class and state ownership", street_panel = "Street-level panel (35 streets in both windows)",
band_100 = "Treatment band 100 feet", band_250 = "Treatment band 250 feet", numbered_only = "Numbered projects only",
bootstrap = "Street bootstrap (1,000 resamples)", placebo = "Placebo: 2024 priority corridors without a project")
fp <- bind_rows(table3[3, ] %>% mutate(check = "primary"), table3[3, ] %>% mutate(check = "poisson"), robust) %>%
mutate(check = factor(check, levels = names(labels_fp))) %>% arrange(check) %>% mutate(y = rev(seq_len(n())))
par(mar = c(3.6, 17.5, 1.2, 7.5), mgp = c(2.2, 0.6, 0))
col_fp <- ifelse(fp$check == "primary", orange, blue)
plot(fp$rate_ratio, fp$y, log = "x", xlim = c(0.6, 2.0), ylim = c(0.5, nrow(fp) + 0.5), pch = 19, col = col_fp, cex = 1.1,
yaxt = "n", xaxt = "n", ylab = "", xlab = "Rate ratio for treated-and-later (log scale)", bty = "l", cex.lab = 0.85)
axis(1, at = c(0.6, 0.8, 1.0, 1.5, 2.0), labels = c("0.6", "0.8", "1.0", "1.5", "2.0"), cex.axis = 0.8, gap.axis = -1)
segments(fp$lower, fp$y, fp$upper, fp$y, col = col_fp, lwd = 2.2)
abline(v = 1, lty = 2, col = "grey40")
axis(2, at = fp$y, labels = labels_fp[as.character(fp$check)], las = 1, cex.axis = 0.72, tick = FALSE, line = -0.4)
text(2.9, fp$y, sprintf("%.2f (%.2f\u2013%.2f)", fp$rate_ratio, fp$lower, fp$upper), xpd = TRUE, adj = 0, cex = 0.72)
text(2.9, nrow(fp) + 1.1, "Rate ratio (95% CI)", xpd = TRUE, adj = 0, cex = 0.72, col = "grey30")
invisible(dev.off())
write.csv(fp %>% dplyr::select(check, rate_ratio, lower, upper, p), file.path(out_dir, "figure3_forest_data.csv"), row.names = FALSE)
out_txt <- capture.output({
cat("Primary model (Poisson, offset log miles), cluster-robust by street\n"); print(table3)
cat("\nLength as covariate\n"); print(table3_len); cat("\nMinimum detectable rate ratio at 80% power:", round(mde80, 3), "\n")
cat("\nRobustness checks\n"); print(robust); cat("\nEmpirical Bayes check\n"); print(eb)
cat("\nDose-response and cross-section\n"); print(table4); cat("\nBicyclist corridors\n"); print(as.data.frame(bike_tab))
cat("\nTable 6 (CRIS)\n"); print(as.data.frame(table6)); cat("\nTable 7 (districts)\n"); print(as.data.frame(table7))
cat("\nsessionInfo()\n"); print(sessionInfo())
})
writeLines(out_txt, file.path(out_dir, "model_output.txt"))
save.image(file.path(out_dir, "vision_zero_results.RData"))
list.files(out_dir)
## [1] "corridors_analysis_table.csv" "corridors_analysis_table.rds" "corridors_with_treatment.geojson"
## [4] "Figure1_project_point_distances.png" "Figure2_ksi_per_mile_slope.png" "figure3_forest_data.csv"
## [7] "Figure3_rate_ratio_forest.png" "model_output.txt" "project_points_nearest_corridor.csv"
## [10] "street_panel_35_streets.csv" "table_empirical_bayes_check.csv" "table_robustness_checks.csv"
## [13] "table1_descriptives.csv" "table2_cells.csv" "table3_length_as_covariate.csv"
## [16] "table3_primary_model.csv" "table4_bicyclist_corridors.csv" "table4_dose_response_cross_section.csv"
## [19] "table6_cris_groups.csv" "table7_by_district.csv" "vision_zero_results.RData"