Allowing for non-linearity

Note

I run the main regressions allowing for non-linear effects of temperature on domestic violence.

I estimate several regressions, with different binned temperature measures (tmean, tmax, tmin…)

  • I save plots for all the regs
  • For the main reg, with bins_tmean_one, I also
    • save a reg table for the main reg with tmean
    • compute summary statistics for the estimates
    • plot the support of bins and save it
library(tidyverse)
library(here) # path control
library(fixest) # estimations with fixed effects
library(modelsummary) # for the modelplot figure
library(patchwork) # assemble the plots
library(knitr) # for kable
library(tictoc)
tic("Total render time")

source("fun_dynamic_bins.R")

# Massive regressions, take ages to run,
# set memory limit v high, also uses swap
mem.maxVSize(vsize = 40000)
[1] 40000
data <-
  read_rds(here("..", "output", "data_reports.rds")) %>% 
  filter(covid_phase == "Pre-pandemic")

data <- data %>% 
  mutate(prec_quintile = ntile(prec, 5),
         rh_quintile = ntile(rh, 5),
         wsp_quintile = ntile(wsp, 5))

Reports with bins_tmean_one

reg_bins_reports <- data %>% 
  mutate(bins_tmean_one = create_dynamic_bins(tmean, width = 1, 0.01)) %>% 
  fepois(reports_dv ~ i(bins_tmean_one, ref = "[16") |
           prec_quintile + rh_quintile + wsp_quintile +
           ageb + year^month + day_of_week + day_of_year,
         cluster = ~ ageb)

I can plot the estimates:

Or also show a table, and save it as reg_bins_reports.tex

reg_bins_reports
Dependent Var.: Reports
Daily Mean Temp. Bin = <11 -0.2213*** (0.0415)
Daily Mean Temp. Bin = [11,12) -0.1432*** (0.0394)
Daily Mean Temp. Bin = [12,13) -0.0764** (0.0292)
Daily Mean Temp. Bin = [13,14) -0.0789*** (0.0240)
Daily Mean Temp. Bin = [14,15) -0.0602** (0.0188)
Daily Mean Temp. Bin = [15,16) -0.0155 (0.0153)
Daily Mean Temp. Bin = [17,18) 0.0151 (0.0139)
Daily Mean Temp. Bin = [18,19) 0.0409** (0.0155)
Daily Mean Temp. Bin = [19,20) 0.1018*** (0.0184)
Daily Mean Temp. Bin = [20,21) 0.0850*** (0.0233)
Daily Mean Temp. Bin = [21,22) 0.0981*** (0.0291)
Daily Mean Temp. Bin = >=22 0.1341*** (0.0371)
Fixed-Effects: ——————-
prec_quintile Yes
rh_quintile Yes
wsp_quintile Yes
ageb Yes
year-month Yes
day_of_week Yes
day_of_year Yes
______________________________ ___________________
S.E.: Clustered by: ageb
Observations 3,624,826
Squared Cor. 0.01334
Pseudo R2 0.05616
BIC 804,317.1

I want to compute some summary stats:

sum_stats <- iplot(reg_bins_reports, only.params = T)$prms %>% 
  as_tibble() %>% 
  mutate(ref = estimate==0) %>% 
  left_join(data %>%
              count(estimate_names = bins_tmean_one) %>% 
              mutate(share_in_pc = 100*n / sum(n))) 

sum_stats %>% 
  #remove first and last, because its not a one degree jump
  slice(2:(n() - 1)) %>% 
  select(estimate_names, estimate) %>% 
  mutate(value_diff = estimate - lag(estimate),
         mean_diff = mean(value_diff, na.rm = T)) %>%
  kable()
estimate_names estimate value_diff mean_diff
[11,12) -0.1432142 NA 0.024133
[12,13) -0.0764416 0.0667726 0.024133
[13,14) -0.0788947 -0.0024531 0.024133
[14,15) -0.0602501 0.0186447 0.024133
[15,16) -0.0155372 0.0447129 0.024133
[16,17) 0.0000000 0.0155372 0.024133
[17,18) 0.0151162 0.0151162 0.024133
[18,19) 0.0408825 0.0257663 0.024133
[19,20) 0.1017882 0.0609058 0.024133
[20,21) 0.0849635 -0.0168247 0.024133
[21,22) 0.0981159 0.0131524 0.024133

I also plot the support of the temp bins and save it as reg_bins_mean_temp_distrib.png

Reports with bins_tmax_one

reg_bins_reports_tmax <- data %>% 
  mutate(bins_tmax_one = create_dynamic_bins(tmax, width = 1, 0.01)) %>% 
  fepois(reports_dv ~ i(bins_tmax_one, ref = "[22") |
           prec_quintile + rh_quintile + wsp_quintile +
           ageb + year^month + day_of_week + day_of_year,
         cluster = ~ ageb)

Plot and save as reg_bins_reports_tmax_1C.png

Reports with bins_tmin_one

reg_bins_reports_tmin <- data %>% 
  mutate(bins_tmin_one = create_dynamic_bins(tmin, width = 1, 0.01)) %>% 
  fepois(reports_dv ~ i(bins_tmin_one, ref = "[13") |
           prec_quintile+ rh_quintile + wsp_quintile +
           ageb + year^month + day_of_week + day_of_year,
         cluster = ~ ageb)

Plot and save as reg_bins_reports_tmin_1C.png

Reports with bins_midnight_temperature

I extract the midnight temperature from the hourly data. It is simply the first of the list of the 24 hourly measures.

data <- data %>% 
  left_join(
    read_rds(here("..", "output", "temporary_data",
                  "tmp_hourly_ageb.rds"))
  ) %>% 
  mutate(midnight_temp = map_dbl(tmp_hourly, ~ .x[1]))

Run the reg:

reg_midnight <- data %>% 
  mutate(bins_midnight_temp = create_dynamic_bins(midnight_temp, width = 1, 0.01)) %>% 
  fepois(reports_dv ~ i(bins_midnight_temp, ref = "[13") |
           prec_quintile + rh_quintile + wsp_quintile +
           ageb + year^month + day_of_week + day_of_year,
         cluster = ~ ageb)

Plot and save as reg_bins_reports_midnight_temp.png

Reports with bins_tmin_one and bins_tmax_one estimated jointly

reg_bins_reports_tmin_tmax <- data %>% 
  mutate(bins_tmin_one = create_dynamic_bins(tmin, width = 1, 0.01),
         bins_tmax_one = create_dynamic_bins(tmax, width = 1, 0.01)) %>% 
  fepois(reports_dv ~ i(bins_tmin_one, ref = "[13") + i(bins_tmax_one, ref = "[22") |
           prec_quintile+ rh_quintile + wsp_quintile +
           ageb + year^month + day_of_week + day_of_year,
         cluster = ~ ageb)

Plot and save as reg_bins_reports_tmin_tmax_jointly.png and the horizontal version as reg_bins_reports_tmin_tmax_jointly_horizontal.png

Rendered on: 2025-04-03 19:20:30
Total render time: 1295.122 sec elapsed