library(tidyverse)
library(here) # path control
library(lubridate) # handle dates
library(fixest) # estimations with fixed effects
library(modelsummary) # for the modelplot figure
data <-
readRDS(here("output",
"final_data_for_analysis.rds")) %>%
filter(covid_phase == "Pre-pandemic")
# QUINTILES!!!!
# all the 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))
# this is the baseline
reports_dv <-
fepois(reports_dv ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_drugs <-
fepois(reports_drugs ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_fraud <-
fepois(reports_fraud ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_homicide <-
fepois(reports_homicide ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_rapes <-
fepois(reports_rapes ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_sexual <-
fepois(reports_sexual ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_suicide <-
fepois(reports_suicide ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_theft <-
fepois(reports_theft ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_threats <-
fepois(reports_threats ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_trust <-
fepois(reports_trust ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
reports_other <-
fepois(reports_other ~ tmean | prec_quintile + rh_quintile + wsp_quintile + col_id + year^month + day_of_week + day_of_year, cluster = ~ col_id, data = data)
etable(mget(ls(pattern = "reports_")), view= F)
etable(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 = "Effect of temperature on different types of crime",
label = "reg_other_types"
)
# Try the semi-parametric on differnet crime types -----
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) |
col_id + year^month + day_of_week + day_of_year,
cluster = ~ col_id)
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) |
col_id + year^month + day_of_week + day_of_year,
cluster = ~ col_id)
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) |
col_id + year^month + day_of_week + day_of_year,
cluster = ~ col_id)
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) |
col_id + year^month + day_of_week + day_of_year,
cluster = ~ col_id)
## plot ----
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))
(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)Types of crime
Note
I run regs on different types of crime