start_time <- Sys.time()
library(tidyverse)
library(here) # path control
library(lubridate) # handle dates
library(knitr)
library(tictoc)
library(fixest) # estimations with fixed effects
library(modelsummary) # for the modelplot figure
tic("Total render time")Types of crime
Note
In this script, I run regs on different types of crime
data <-
read_rds(here("..", "output", "data_reports.rds")) %>%
filter(covid_phase == "Pre-pandemic")
# quintiles are computed last minute, on the selected period
data <- data %>%
mutate(prec_quintile = ntile(prec, 5),
rh_quintile = ntile(rh, 5),
wsp_quintile = ntile(wsp, 5))here are all the categories of crimes I track
data %>%
filter(covid_phase == "Pre-pandemic") %>%
select(reports_dv, reports_drugs:reports_other) %>%
summarise(across(everything(), ~sum(.x))) %>%
pivot_longer(everything(), names_to = "type", values_to = "total") %>%
mutate(proportion = total/sum(total),
prop_in_pc = round(100*proportion, 0)) %>%
kable()| type | total | proportion | prop_in_pc |
|---|---|---|---|
| reports_dv | 84217 | 0.0919063 | 9 |
| reports_drugs | 16969 | 0.0185183 | 2 |
| reports_feminicide | 214 | 0.0002335 | 0 |
| reports_fraud | 52885 | 0.0577136 | 6 |
| reports_homicide | 7027 | 0.0076686 | 1 |
| reports_rapes | 3784 | 0.0041295 | 0 |
| reports_sexual | 11548 | 0.0126024 | 1 |
| reports_suicide | 1789 | 0.0019523 | 0 |
| reports_theft | 431996 | 0.4714389 | 47 |
| reports_threats | 45885 | 0.0500745 | 5 |
| reports_trust | 15564 | 0.0169851 | 2 |
| reports_other | 244457 | 0.2667769 | 27 |
Run the regressions separately, I need to do it so for memory.
# this is the baseline
reports_dv <-
fepois(reports_dv ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_drugs <-
fepois(reports_drugs ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_fraud <-
fepois(reports_fraud ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_homicide <-
fepois(reports_homicide ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_rapes <-
fepois(reports_rapes ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_sexual <-
fepois(reports_sexual ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_suicide <-
fepois(reports_suicide ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_theft <-
fepois(reports_theft ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_threats <-
fepois(reports_threats ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_trust <-
fepois(reports_trust ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
reports_other <-
fepois(reports_other ~ tmean | prec_quintile + rh_quintile + wsp_quintile + ageb + year^month + day_of_week + day_of_year, cluster = ~ ageb, data = data)
etable(mget(ls(pattern = "reports_")), view= F) reports_drugs reports_dv reports_fraud
Dependent Var.: reports_drugs reports_dv reports_fraud
tmean -0.0009 (0.0083) 0.0274*** (0.0034) 0.0052 (0.0046)
Fixed-Effects: ---------------- ------------------ ---------------
prec_quintile Yes Yes Yes
rh_quintile Yes Yes Yes
wsp_quintile Yes Yes Yes
ageb Yes Yes Yes
year-month Yes Yes Yes
day_of_week Yes Yes Yes
day_of_year Yes Yes Yes
_______________ ________________ __________________ _______________
S.E.: Clustered by: ageb by: ageb by: ageb
Observations 3,058,670 3,624,826 3,454,739
Squared Cor. 0.04594 0.01333 0.09170
Pseudo R2 0.15917 0.05612 0.19726
BIC 214,333.0 804,221.2 487,853.5
reports_homicide reports_other reports_rapes
Dependent Var.: reports_homicide reports_other reports_rapes
tmean 0.0268* (0.0111) 0.0109*** (0.0020) 0.0240 (0.0166)
Fixed-Effects: ---------------- ------------------ ---------------
prec_quintile Yes Yes Yes
rh_quintile Yes Yes Yes
wsp_quintile Yes Yes Yes
ageb Yes Yes Yes
year-month Yes Yes Yes
day_of_week Yes Yes Yes
day_of_year Yes Yes Yes
_______________ ________________ __________________ _______________
S.E.: Clustered by: ageb by: ageb by: ageb
Observations 2,843,215 3,633,862 2,493,915
Squared Cor. 0.00285 0.09356 0.00170
Pseudo R2 0.05332 0.11949 0.05402
BIC 127,678.8 1,667,373.4 84,408.0
reports_sexual reports_suicide reports_theft
Dependent Var.: reports_sexual reports_suicide reports_theft
tmean 0.0049 (0.0091) 0.0136 (0.0221) 0.0023 (0.0015)
Fixed-Effects: --------------- --------------- ---------------
prec_quintile Yes Yes Yes
rh_quintile Yes Yes Yes
wsp_quintile Yes Yes Yes
ageb Yes Yes Yes
year-month Yes Yes Yes
day_of_week Yes Yes Yes
day_of_year Yes Yes Yes
_______________ _______________ _______________ _______________
S.E.: Clustered by: ageb by: ageb by: ageb
Observations 3,230,260 1,679,139 3,636,874
Squared Cor. 0.00558 0.00091 0.14487
Pseudo R2 0.07082 0.03776 0.14511
BIC 181,111.3 49,247.3 2,433,237.9
reports_threats reports_trust
Dependent Var.: reports_threats reports_trust
tmean 0.0322*** (0.0047) 0.0115 (0.0075)
Fixed-Effects: ------------------ ---------------
prec_quintile Yes Yes
rh_quintile Yes Yes
wsp_quintile Yes Yes
ageb Yes Yes
year-month Yes Yes
day_of_week Yes Yes
day_of_year Yes Yes
_______________ __________________ _______________
S.E.: Clustered by: ageb by: ageb
Observations 3,575,128 3,234,781
Squared Cor. 0.00945 0.01094
Pseudo R2 0.05819 0.08588
BIC 506,967.0 219,961.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1etable(reports_dv, #reg_linear_calls,
reports_homicide,
# reports_suicide,
# reports_rapes,
reports_theft,
reports_fraud,
# reports_drugs,
dict = c(reports_dv = "Domestic violence",
reports_theft = "Theft",
reports_fraud = "Fraud",
reports_drugs = "Drugs possession",
reports_homicide = "Homicide",
reports_rapes = "Rapes",
reports_suicide = "Suicide"),
fitstat = ~ n + my,
digits.stats = "s3", digits = "r4",
file = here("..", "output", "tables", "reg_other_types.tex"), replace = T,
view = T,
# fixef_sizes = T,
fixef.group=list("Prec, hum, wsp quintiles"="quintile",
"Date FEs"="Month|Day|Year"),
style.tex = style.tex("aer",
model.format = "(i)",
tpt = TRUE,
notes.tpt.intro = "\\footnotesize",
fixef.suffix = " FEs"),
adjustbox = TRUE, float = T,
title = "Temperature Effects on Reports for Other Types of Crime",
label = "reg_other_types",
notes = "Notes: Standard errors clustered by neighborhood in parentheses. * denotes significance at the 10% level, ** at the 5% level, and *** at the 1% level. Each column reports results from a separate regression of daily crime counts on temperature and weather controls. All regressions include neighborhood fixed effects, day-of-year fixed effects, year fixed effects, and quintile controls for precipitation, humidity, and wind speed. The dependent variable varies across columns and corresponds to daily counts per neighborhood for the indicated crime category."
)Run the semi-parametric on differnet crime types
I do not run it for the time being, it takes way too long.
source(here("codes", "my_functions.R"))
## regs for bins ----
reports_dv <- data %>%
filter(covid_phase == "Pre-pandemic") %>%
mutate(bins_tmean_one = create_dynamic_bins(tmean, width = 1, 0.01)) %>%
fepois(reports_dv ~ i(bins_tmean_one, ref = "[16") +
i(prec_quintile) + i(rh_quintile) + i(wsp_quintile) |
ageb + year^month + day_of_week + day_of_year,
cluster = ~ ageb)
reports_homicide <- data %>%
filter(covid_phase == "Pre-pandemic") %>%
mutate(bins_tmean_one = create_dynamic_bins(tmean, width = 1, 0.01)) %>%
fepois(reports_homicide ~ i(bins_tmean_one, ref = "[16") +
i(prec_quintile) + i(rh_quintile) + i(wsp_quintile) |
ageb + year^month + day_of_week + day_of_year,
cluster = ~ ageb)
reports_theft <- data %>%
filter(covid_phase == "Pre-pandemic") %>%
mutate(bins_tmean_one = create_dynamic_bins(tmean, width = 1, 0.01)) %>%
fepois(reports_theft ~ i(bins_tmean_one, ref = "[16") +
i(prec_quintile) + i(rh_quintile) + i(wsp_quintile) |
ageb + year^month + day_of_week + day_of_year,
cluster = ~ ageb)
reports_fraud <- data %>%
filter(covid_phase == "Pre-pandemic") %>%
mutate(bins_tmean_one = create_dynamic_bins(tmean, width = 1, 0.01)) %>%
fepois(reports_fraud ~ i(bins_tmean_one, ref = "[16") +
i(prec_quintile) + i(rh_quintile) + i(wsp_quintile) |
ageb + year^month + day_of_week + day_of_year,
cluster = ~ ageb)Plot all those:
library(patchwork)
breaks = seq(-0.3, 0.2, 0.1)
labels = scales::label_number(accuracy = 0.1)
plot_dv <- iplot(reports_dv, only.params = T)$prms %>%
ggplot(aes(x = estimate_names, y = estimate, ymin = ci_low, ymax = ci_high, group = 1)) +
geom_errorbar(width = .2) + geom_point(size = 2) + geom_hline(yintercept = 0, alpha = .5) +
theme_light() +
labs(title = "Reports for DV", x = "Daily mean temperature", y = "Estimate and 95% CI") +
scale_y_continuous(breaks = breaks, labels = labels) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
coord_cartesian(ylim = c(-0.3, 0.2))
plot_hom <- iplot(reports_homicide, only.params = T)$prms %>%
ggplot(aes(x = estimate_names, y = estimate, ymin = ci_low, ymax = ci_high, group = 1)) +
geom_errorbar(width = .2) + geom_point(size = 2) + geom_hline(yintercept = 0, alpha = .5) +
theme_light() +
labs(title = "Reports for homicides", x = "Daily mean temperature", y = "Estimate and 95% CI") +
scale_y_continuous(breaks = breaks, labels = labels) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
coord_cartesian(ylim = c(-0.3, 0.2))
plot_theft <- iplot(reports_theft, only.params = T)$prms %>%
ggplot(aes(x = estimate_names, y = estimate, ymin = ci_low, ymax = ci_high, group = 1)) +
geom_errorbar(width = .2) + geom_point(size = 2) + geom_hline(yintercept = 0, alpha = .5) +
theme_light() +
labs(title = "Reports for thefts", x = "Daily mean temperature", y = "Estimate and 95% CI") +
scale_y_continuous(breaks = breaks, labels = labels) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
coord_cartesian(ylim = c(-0.3, 0.2))
plot_fraud <- iplot(reports_fraud, only.params = T)$prms %>%
ggplot(aes(x = estimate_names, y = estimate, ymin = ci_low, ymax = ci_high, group = 1)) +
geom_errorbar(width = .2) + geom_point(size = 2) + geom_hline(yintercept = 0, alpha = .5) +
theme_light() +
labs(title = "Reports for frauds", x = "Daily mean temperature", y = "Estimate and 95% CI") +
scale_y_continuous(breaks = breaks, labels = labels) +
scale_x_discrete(guide = guide_axis(angle = 45)) +
coord_cartesian(ylim = c(-0.3, 0.2))Save the plots
(plot_dv | plot_hom) / (plot_theft | plot_fraud)
ggsave(path = here("output", "figures"),
filename = "reg_types_combined.png",
width = 12, height = 9)
# separately
ggsave(plot_dv, path = here("output", "figures"),
filename = "reg_types_dv.png",
width = 6, height = 4)
ggsave(plot_hom, path = here("output", "figures"),
filename = "reg_types_hom.png",
width = 6, height = 4)
ggsave(plot_theft, path = here("output", "figures"),
filename = "reg_types_theft.png",
width = 6, height = 4)
ggsave(plot_fraud, path = here("output", "figures"),
filename = "reg_types_fraud.png",
width = 6, height = 4)Rendered on: 2025-07-27 03:55:17Total render time: 477.206 sec elapsed