Step 1. Weather variables. We defined the dates by year and availability (whether the event took place or not), days of the week, time windows, and pre and post-event spans, which we then linked to nearby meteorological stations. Also we linked with Collisions data.
# remove objects and memory
rm(list=ls());gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 836662 44.7 1628492 87 1179400 63.0
Vcells 1692015 13.0 8388608 64 3370825 25.8#remove images
while(!dev.cur())dev.off()
cat("\014")load(paste0(getwd(),"/_data/step1.RData"))Particularly, the weathercan library, which contains information on climate change and environmental variables from stations across Canada1.
#borrar caché
#system("fc-cache -f -v")
#check R version
if(Sys.info()["sysname"]=="Windows"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
if(Sys.info()["sysname"]=="Linux"){
if (getRversion() != "4.4.1") { stop("Requiere versión de R 4.4.1. Actual: ", getRversion()) }
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
# install.packages(c("dplyr", #for data
# "tidyr", #for data
# "lubridate", #for dates
# "openxlsx", #for excel files
# "rio", #for importing and exporting data
# "purrr", #for iterating in databases
# "devtools", #for external packages
# "DiagrammeR", #para visualizar DAG
# "dagitty",
# "ggdag",
# "ggplot2" #for graphics
# "kableExtra", #pretty tables
# "quarto", #for documents
# "geosphere" #for coordinates and classifying
# "geepack", #for regression
# "glmmTMB", #For GLMMs
# "DHARMa", #For residual diagnostics
# "car", #For hypothesis testing
# "brms", #Bayesian model
# "bayesplot",
# "loo",
# "Synth", #for synthetic control method
# "weathercan",#for weather data
# "sandwich", #cluster robust intervals
# "emmeans", #for predictions
# "gnm", #Conditional Poisson models
# "splines", #nonlinearity
# "geeM", #negative binomial and more flexible GEE models
# "PanelMatch", #Matching technique with panel data
# "scpi" #control sintético
# "nixtlar", #for time series analysis and prediction
# "CausalImpact" #for time series causal impact
# "forecast" #for time series analysis prediction and decomposition
# ))
library(dplyr); library(lubridate); library(tidyverse); library(openxlsx); library(rio); library(purrr); library(dagitty); library(ggdag); library(kableExtra); library(geosphere); library(geepack); library(lme4); library(glmmTMB); library(DHARMa); library(car); library(brms); library(bayesplot); library(loo); library(Synth); library(weathercan); library(sandwich); library(emmeans); library(gnm); library(splines); library(geeM); library(plm); library(PanelMatch); library(scpi); library(fect); library(nixtlar); library(CausalImpact); library(forecast); library(webshot); library(tableone)
Adjuntando el paquete: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Adjuntando el paquete: 'lubridate'The following objects are masked from 'package:base':
date, intersect, setdiff, unionWarning: package 'readr' was built under R version 4.4.3── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ purrr 1.0.4 ✔ tidyr 1.3.1
✔ readr 2.1.5 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Adjuntando el paquete: 'dagitty'
The following object is masked from 'package:rio':
convert
Adjuntando el paquete: 'ggdag'
The following object is masked from 'package:stats':
filter
Adjuntando el paquete: 'kableExtra'
The following object is masked from 'package:dplyr':
group_rows
Cargando paquete requerido: Matrix
Adjuntando el paquete: 'Matrix'
The following objects are masked from 'package:tidyr':
expand, pack, unpack
Adjuntando el paquete: 'lme4'
The following object is masked from 'package:rio':
factorize
This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
Cargando paquete requerido: carData
Adjuntando el paquete: 'car'
The following object is masked from 'package:purrr':
some
The following object is masked from 'package:dplyr':
recode
Cargando paquete requerido: Rcpp
Loading 'brms' package (version 2.22.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Adjuntando el paquete: 'brms'
The following object is masked from 'package:glmmTMB':
lognormal
The following object is masked from 'package:lme4':
ngrps
The following object is masked from 'package:stats':
ar
This is bayesplot version 1.12.0
- Online documentation and vignettes at mc-stan.org/bayesplot
- bayesplot theme set to bayesplot::theme_default()
* Does _not_ affect other ggplot2 plots
* See ?bayesplot_theme_set for details on theme setting
Adjuntando el paquete: 'bayesplot'
The following object is masked from 'package:brms':
rhat
This is loo version 2.8.0
- Online documentation and vignettes at mc-stan.org/loo
- As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
- Windows 10 users: loo may be very slow if 'mc.cores' is set in your .Rprofile file (see https://github.com/stan-dev/loo/issues/94).
##
## Synth Package: Implements Synthetic Control Methods.
## See https://web.stanford.edu/~jhain/synthpage.html for additional information.Warning: package 'weathercan' was built under R version 4.4.3As of v0.7.2, the `normals` column in `stations()` reflects whether or not there
are *any* normals available (not just the most recent).
Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'
Registered S3 method overwritten by 'lfe':
method from
nobs.felm broom
Adjuntando el paquete: 'plm'
The following objects are masked from 'package:dplyr':
between, lag, lead
Adjuntando el paquete: 'PanelMatch'
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:stats':
weights
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Cargando paquete requerido: bsts
Cargando paquete requerido: BoomSpikeSlab
Cargando paquete requerido: Boom
Adjuntando el paquete: 'Boom'
The following objects are masked from 'package:brms':
ddirichlet, rdirichlet
The following object is masked from 'package:stats':
rWishart
Adjuntando el paquete: 'BoomSpikeSlab'
The following object is masked from 'package:stats':
knots
Cargando paquete requerido: zoo
Adjuntando el paquete: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Cargando paquete requerido: xts
######################### Warning from 'xts' package ##########################
# #
# The dplyr lag() function breaks how base R's lag() function is supposed to #
# work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
# source() into this session won't work correctly. #
# #
# Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
# conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
# dplyr from breaking base R's lag() function. #
# #
# Code in packages is not affected. It's protected by R's namespace mechanism #
# Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
# #
###############################################################################
Adjuntando el paquete: 'xts'
The following objects are masked from 'package:dplyr':
first, last
Adjuntando el paquete: 'bsts'
The following object is masked from 'package:BoomSpikeSlab':
SuggestBurn
Registered S3 method overwritten by 'quantmod':
method from
as.zoo.data.frame zoo
Adjuntando el paquete: 'forecast'
The following object is masked from 'package:brms':
ma#special repository indicated or the package
if(!require(weathercan)){
install.packages("weathercan",
repos = c("https://ropensci.r-universe.dev", "https://cloud.r-project.org")); library(weathercan)
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
if(!require(bpmn)){devtools::install_github("bergant/bpmn")}Cargando paquete requerido: bpmn#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
sum_dates <- function(x){
cbind.data.frame(
min= as.Date(min(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01"),
p001= as.Date(quantile(unclass(as.Date(x)), .001, na.rm=T), origin = "1970-01-01"),
p005= as.Date(quantile(unclass(as.Date(x)), .005, na.rm=T), origin = "1970-01-01"),
p025= as.Date(quantile(unclass(as.Date(x)), .025, na.rm=T), origin = "1970-01-01"),
p25= as.Date(quantile(unclass(as.Date(x)), .25, na.rm=T), origin = "1970-01-01"),
p50= as.Date(quantile(unclass(as.Date(x)), .5, na.rm=T), origin = "1970-01-01"),
p75= as.Date(quantile(unclass(as.Date(x)), .75, na.rm=T), origin = "1970-01-01"),
p975= as.Date(quantile(unclass(as.Date(x)), .975, na.rm=T), origin = "1970-01-01"),
p995= as.Date(quantile(unclass(as.Date(x)), .995, na.rm=T), origin = "1970-01-01"),
p999= as.Date(quantile(unclass(as.Date(x)), .999, na.rm=T), origin = "1970-01-01"),
max= as.Date(max(unclass(as.Date(x)), na.rm=T), origin = "1970-01-01")
)
}
smd_bin <- function(x,y){
z <- x*(1-x)
t <- y*(1-y)
k <- sum(z,t)
l <- k/2
return((x-y)/sqrt(l))
}
theme_custom_sjplot2 <- function(base_size = 12, base_family = "") {
theme_minimal(base_size = base_size, base_family = base_family) +
theme(
# Text elements
text = element_text(size = base_size, family = base_family),
plot.title = element_text(face = "bold", hjust = 0.5, size = base_size * 1.2),
plot.subtitle = element_text(hjust = 0.5, margin = margin(b = 10)),
axis.title = element_text(size = base_size, face = "bold"),
axis.text = element_text(size = base_size * 0.8),
axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.5),
axis.text.y = element_text(angle = 0, hjust = 1, vjust = 0.5),
axis.title.x = element_text(margin = margin(t = 10)),
axis.title.y = element_text(margin = margin(r = 10)),
# Plot layout
plot.margin = margin(t = 20, r = 20, b = 20, l = 20),
panel.grid.major = element_line(color = "grey80"),
panel.grid.minor = element_blank(),
legend.position = "right",
legend.text = element_text(size = base_size * 0.8),
legend.title = element_text(size = base_size, face = "bold"),
legend.background = element_rect(fill = "white", colour = NA),
legend.box.background = element_rect(colour = "grey80", linetype = "solid"),
legend.key = element_rect(fill = "white", colour = "white")
)
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
num_cores <- parallel::detectCores() -1
data.table::setDTthreads(threads = num_cores)#restore_after_fork = NULL, throttle = NULL)
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#CONFIG #######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
options(scipen=2) #display numbers rather scientific numberWe collapsed data by aggregating it at the day–year level for treatment and control areas. For each cell, we summed the number of high-speed tickets and documented vehicles and license holders, and averaged the minimum temperature, maximum temperature, while the 2-day-lagged precipitation and total daily precipitation were summarized by the median. We restricted post-outcome period to 7 days.
ordered<-
paste(rep(c("control","treatment"), each = length(2000:2022)),rep(2000:2022, times = 2), sep=".")
ordered_with_placebo_years<-
c("control.2000", "control.2001", "control.2002", "control.2003",
"control.2004", "control.2005", "control.2006", "control.2007", "control.2008",
"control.2009", "control.2010", "control.2011", "control.2012", "control.2013",
"control.2014", "control.2015", "control.2016", "control.2017", "control.2018",
"control.2019", "control.2020", "control.2021", "control.2022", "treatment.2009",
"treatment.2020", "treatment.2021", "treatment.2000", "treatment.2001", "treatment.2002",
"treatment.2003", "treatment.2004", "treatment.2005", "treatment.2006",
"treatment.2007", "treatment.2008", "treatment.2010", "treatment.2011",
"treatment.2012", "treatment.2013", "treatment.2014", "treatment.2015",
"treatment.2016", "treatment.2017", "treatment.2018", "treatment.2019",
"treatment.2022")
collisions_weather_corr_rect_synth <-
collisions_weather_corr_rect |>
filter(date_num <= unclass(race_date_rec) + 7) |>
group_by(tr_contr_sens, year, yday_corr, date) |>
summarise(
sum_velocidad = sum(velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(velocidad_lead1_post_conv_imp, na.rm = T),
sum_velocidad_lead2= sum(velocidad_lead2_post_conv_imp, na.rm = T),
exp_veh = sum(vehicles_use_type, na.rm = TRUE), # Nº de vehículos expuestos
exp_lic = sum(license_holders_sex_age, na.rm = TRUE), # Nº de conductores
median_lag_2_prec_median_imp = median(lag_2_prec_median_lin_imp, na.rm = TRUE),
mean_min_temp_mean_lin = mean(min_temp_mean_lin, na.rm = TRUE),
mean_max_temp_mean_lin = mean(max_temp_mean_lin, na.rm = TRUE),
median_total_precip_median_lin = median(total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) |>
mutate( # ── aux variables
rate_veh = (sum_velocidad / exp_veh) * 1e6, # not to model
rate_lic = (sum_velocidad / exp_lic) * 1e6, # not to model
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh), # offset
off_lic = log(exp_lic)
) |>
left_join(races, by = c("date" = "race_date")) |>
mutate(
race_day = if_else(!is.na(year.y), 1L, 0L),
id = as.numeric(factor(paste0(tr_contr_sens, ".",year.x), levels= ordered_with_placebo_years)),
id_lab = factor(paste0(tr_contr_sens, ".",year.x))
) |>
select(-year.y) |>
mutate(id_lab= case_when(id_lab== "treatment.2009"~ "placebo.2009",
id_lab== "treatment.2020"~ "placebo.2020",
id_lab== "treatment.2021"~ "placebo.2021",
T~id_lab)) |>
mutate(tr_contrl_corr= ifelse(grepl("treatment",id_lab),1,0))
if(identical(collisions_weather_corr_rect_synth$median_total_precip_median_lin, collisions_weather_corr_rect_synth$median_lag_2_prec_median_imp)){warning("Identical lagged with not lagged variable")}
cat("Day of the race \n")
collisions_weather_corr_rect_synth |>
filter(race_day==1) |> reframe(min=yday_corr, max=yday_corr) |> print(n=46)
############################################################
## 1. “Only-treated, one series” ─────────────────────────
## (collisions_weather_corr_rect_synth_tr)
############################################################
cat("Generate a summary of counts of every year in treated only\n")
collisions_weather_corr_rect_synth_tr <-
collisions_weather_corr_rect_synth %>% # start from the long table
filter(grepl("treatment", id_lab)) %>% # keep treated rows only
group_by(yday_corr) %>% # collapse across MRC × years
summarise(
sum_velocidad = sum(sum_velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(sum_velocidad_lead1, na.rm = T),
sum_velocidad_lead2= sum(sum_velocidad_lead2, na.rm = T),
exp_veh = sum(exp_veh, na.rm = TRUE),
exp_lic = sum(exp_lic, na.rm = TRUE),
mean_median_lag_2_prec_median_imp = mean(median_lag_2_prec_median_imp, na.rm = TRUE),
mean_min_temp_mean_lin = mean(mean_min_temp_mean_lin, na.rm = TRUE),
mean_max_temp_mean_lin = mean(mean_max_temp_mean_lin, na.rm = TRUE),
mean_median_total_precip_median_lin = mean(median_total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
rate_veh = (sum_velocidad / exp_veh) * 1e6,
rate_lic = (sum_velocidad / exp_lic) * 1e6,
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh),
off_lic = log(exp_lic),
year.x = 2000, # dummy year / date just to keep structure
date = as.Date("2019-05-19"),
race_day = 0L,
id = 46 # single ID for the synthetic treated series
)
############################################################
## 2. “One treated - many controls” ─────────────────────
## (collisions_weather_corr_rect_synth_one_tr_many_cntrls)
############################################################
# controls left at their original resolution (one ID per MRC–year)
controls_long <- collisions_weather_corr_rect_synth %>%
filter(!grepl("treatment", id_lab)) %>%
rename(
mean_median_lag_2_prec_median_imp = median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin = median_total_precip_median_lin
)
collisions_weather_corr_rect_synth_one_tr_many_cntrls <- bind_rows(
controls_long,
collisions_weather_corr_rect_synth_tr # the single treated series from step 1
)
############################################################
## 3. “Only-controls, one series” ───────────────────────
## (collisions_weather_corr_rect_synth_cntr)
############################################################
collisions_weather_corr_rect_synth_cntr <-
collisions_weather_corr_rect_synth %>%
filter(!grepl("treatment", id_lab)) %>% # keep treated rows only
group_by(yday_corr) %>% # collapse across MRC × years
summarise(
sum_velocidad = sum(sum_velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(sum_velocidad_lead1, na.rm = T),
sum_velocidad_lead2= sum(sum_velocidad_lead2, na.rm = T),
exp_veh = sum(exp_veh, na.rm = TRUE),
exp_lic = sum(exp_lic, na.rm = TRUE),
# distributional summaries to keep a sense of variability
p25_lag_2_prec_median_imp = quantile(median_lag_2_prec_median_imp, .25, na.rm = TRUE),
p75_lag_2_prec_median_imp = quantile(median_lag_2_prec_median_imp, .75, na.rm = TRUE),
mean_median_lag_2_prec_median_imp = mean(median_lag_2_prec_median_imp, na.rm = TRUE),
p25_min_temp_mean_lin = quantile(mean_min_temp_mean_lin, .25, na.rm = TRUE),
p75_min_temp_mean_lin = quantile(mean_min_temp_mean_lin, .75, na.rm = TRUE),
mean_min_temp_mean_lin = mean(mean_min_temp_mean_lin, na.rm = TRUE),
p25_max_temp_mean_lin = quantile(mean_max_temp_mean_lin, .25, na.rm = TRUE),
p75_max_temp_mean_lin = quantile(mean_max_temp_mean_lin, .75, na.rm = TRUE),
mean_max_temp_mean_lin = mean(mean_max_temp_mean_lin, na.rm = TRUE),
p25_total_precip_median_lin = quantile(median_total_precip_median_lin, .25, na.rm = TRUE),
p75_total_precip_median_lin = quantile(median_total_precip_median_lin, .75, na.rm = TRUE),
mean_median_total_precip_median_lin = mean(median_total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
rate_veh = (sum_velocidad / exp_veh) * 1e6,
rate_lic = (sum_velocidad / exp_lic) * 1e6,
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh),
off_lic = log(exp_lic),
year.x = 2000,
date = as.Date("2019-05-19"),
race_day = 0L,
id = 1 # single ID for the synthetic control
) %>%
# rename the median variables so their names match the treated table
rename(
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin = mean_median_total_precip_median_lin
)
############################################################
## Quick sanity checks
############################################################
# glimpse(collisions_weather_corr_rect_synth_tr)
# glimpse(collisions_weather_corr_rect_synth_one_tr_many_cntrls)
# glimpse(collisions_weather_corr_rect_synth_cntr)
cat("Check imbalances\n")
collisions_weather_corr_rect_synth_one_tr_many_cntrls |> group_by(id) |> reframe(max=max(yday_corr), min=min(yday_corr), n=n()) |> tail()
invisible("Consolidate one treated and one control")
collisions_weather_corr_rect_synth_cntr_tr <-
rbind.data.frame(collisions_weather_corr_rect_synth_cntr,collisions_weather_corr_rect_synth_tr |>
# rename(collisions_weather_corr_rect_synth_tr, "median_lag_2_prec_mean_imp"= "mean_lag_2_prec_mean_imp", "median_total_precip_mean_lin"= "mean_total_precip_mean_lin")|>
mutate(p25_lag_2_prec_median_imp=0, p75_lag_2_prec_median_imp=0, p25_min_temp_mean_lin=0, p75_min_temp_mean_lin=0, p25_max_temp_mean_lin=0, p75_max_temp_mean_lin=0, p25_total_precip_median_lin=0, p75_total_precip_median_lin=0))|> data.frame() |>
mutate(treated_cntrl= ifelse(id==46, "treated","control"))Day of the race
# A tibble: 40 × 2
min max
<int> <int>
1 38 38
2 38 38
3 38 38
4 38 38
5 38 38
6 38 38
7 38 38
8 38 38
9 38 38
10 38 38
11 38 38
12 38 38
13 38 38
14 38 38
15 38 38
16 38 38
17 38 38
18 38 38
19 38 38
20 38 38
21 38 38
22 38 38
23 38 38
24 38 38
25 38 38
26 38 38
27 38 38
28 38 38
29 38 38
30 38 38
31 38 38
32 38 38
33 38 38
34 38 38
35 38 38
36 38 38
37 38 38
38 38 38
39 38 38
40 38 38
Generate a summary of counts of every year in treated only
Check imbalances
# A tibble: 6 × 4
id max min n
<dbl> <int> <int> <int>
1 22 45 1 45
2 23 45 1 45
3 24 45 1 45
4 25 45 1 45
5 26 45 1 45
6 46 45 1 45collisions_weather_corr_rect_quebec_synth <-
collisions_weather_corr_rect |>
filter(date_num <= unclass(race_date_rec) + 7) |>
group_by(tr_contr, year, yday_corr, date) |>
summarise(
sum_velocidad = sum(velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(velocidad_lead1_post_conv_imp, na.rm = T),
sum_velocidad_lead2= sum(velocidad_lead2_post_conv_imp, na.rm = T),
exp_veh = sum(vehicles_use_type, na.rm = TRUE), # Nº de vehículos expuestos
exp_lic = sum(license_holders_sex_age, na.rm = TRUE), # Nº de conductores
median_lag_2_prec_median_imp = median(lag_2_prec_median_lin_imp, na.rm = TRUE),
mean_min_temp_mean_lin = mean(min_temp_mean_lin, na.rm = TRUE),
mean_max_temp_mean_lin = mean(max_temp_mean_lin, na.rm = TRUE),
median_total_precip_median_lin = median(total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) |>
mutate( # ── aux variables
rate_veh = (sum_velocidad / exp_veh) * 1e6, # not to model
rate_lic = (sum_velocidad / exp_lic) * 1e6, # not to model
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh), # offset
off_lic = log(exp_lic)
) |>
left_join(races, by = c("date" = "race_date"))|>
mutate(
race_day = if_else(!is.na(year.y), 1L, 0L),
id = as.numeric(factor(paste0(tr_contr, ".",year.x), levels= ordered_with_placebo_years)),
id_lab = factor(paste0(tr_contr, ".",year.x))
) |>
select(-year.y) |>
mutate(id_lab= case_when(id_lab== "treatment.2009"~ "placebo.2009",
id_lab== "treatment.2020"~ "placebo.2020",
id_lab== "treatment.2021"~ "placebo.2021",
T~id_lab)) |>
mutate(tr_contrl_corr= ifelse(grepl("treatment",id_lab),1,0))
collisions_weather_corr_rect_quebec_synth_tr <-
collisions_weather_corr_rect_quebec_synth %>% # start from the long table
filter(grepl("treatment", id_lab)) %>% # keep treated rows only
group_by(yday_corr) %>% # collapse across MRC × years
summarise(
sum_velocidad = sum(sum_velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(sum_velocidad_lead1, na.rm = T),
sum_velocidad_lead2= sum(sum_velocidad_lead2, na.rm = T),
exp_veh = sum(exp_veh, na.rm = TRUE),
exp_lic = sum(exp_lic, na.rm = TRUE),
mean_median_lag_2_prec_median_imp = mean(median_lag_2_prec_median_imp, na.rm = TRUE),
mean_min_temp_mean_lin = mean(mean_min_temp_mean_lin, na.rm = TRUE),
mean_max_temp_mean_lin = mean(mean_max_temp_mean_lin, na.rm = TRUE),
mean_median_total_precip_median_lin = mean(median_total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
rate_veh = (sum_velocidad / exp_veh) * 1e6,
rate_lic = (sum_velocidad / exp_lic) * 1e6,
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh),
off_lic = log(exp_lic),
year.x = 2000, # dummy year / date just to keep structure
date = as.Date("2019-05-19"),
race_day = 0L,
id = 46 # single ID for the synthetic treated series
)
controls_quebec_long <- collisions_weather_corr_rect_quebec_synth %>%
filter(!grepl("treatment", id_lab)) %>%
rename(
mean_median_lag_2_prec_median_imp = median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin = median_total_precip_median_lin
)
collisions_weather_corr_rect_quebec_synth_one_tr_many_cntrls <- bind_rows(
controls_quebec_long,
collisions_weather_corr_rect_quebec_synth_tr # the single treated series from step 1
)
collisions_weather_corr_rect_quebec_synth_cntr <-
collisions_weather_corr_rect_quebec_synth %>%
filter(!grepl("treatment", id_lab)) %>% # keep treated rows only
group_by(yday_corr) %>% # collapse across all control IDs
summarise(
sum_velocidad = sum(sum_velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(sum_velocidad_lead1, na.rm = T),
sum_velocidad_lead2= sum(sum_velocidad_lead2, na.rm = T),
exp_veh = sum(exp_veh, na.rm = TRUE),
exp_lic = sum(exp_lic, na.rm = TRUE),
# distributional summaries to keep a sense of variability
p25_lag_2_prec_median_imp = quantile(median_lag_2_prec_median_imp, .25, na.rm = TRUE),
p75_lag_2_prec_median_imp = quantile(median_lag_2_prec_median_imp, .75, na.rm = TRUE),
mean_median_lag_2_prec_median_imp = mean(median_lag_2_prec_median_imp, na.rm = TRUE),
p25_min_temp_mean_lin = quantile(mean_min_temp_mean_lin, .25, na.rm = TRUE),
p75_min_temp_mean_lin = quantile(mean_min_temp_mean_lin, .75, na.rm = TRUE),
mean_min_temp_mean_lin = mean(mean_min_temp_mean_lin, na.rm = TRUE),
p25_max_temp_mean_lin = quantile(mean_max_temp_mean_lin, .25, na.rm = TRUE),
p75_max_temp_mean_lin = quantile(mean_max_temp_mean_lin, .75, na.rm = TRUE),
mean_max_temp_mean_lin = mean(mean_max_temp_mean_lin, na.rm = TRUE),
p25_total_precip_median_lin = quantile(median_total_precip_median_lin, .25, na.rm = TRUE),
p75_total_precip_median_lin = quantile(median_total_precip_median_lin, .75, na.rm = TRUE),
mean_median_total_precip_median_lin = mean(median_total_precip_median_lin, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
rate_veh = (sum_velocidad / exp_veh) * 1e6,
rate_lic = (sum_velocidad / exp_lic) * 1e6,
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh),
off_lic = log(exp_lic),
year.x = 2000,
date = as.Date("2019-05-19"),
race_day = 0L,
id = 1 # single ID for the synthetic control
) %>%
# rename the median variables so their names match the treated table
rename(
mean_median_lag_2_prec_median_imp = mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin = mean_median_total_precip_median_lin
)
collisions_weather_corr_rect_quebec_synth_cntr_tr <-
rbind.data.frame(collisions_weather_corr_rect_quebec_synth_cntr,collisions_weather_corr_rect_quebec_synth_tr |>
mutate(p25_lag_2_prec_median_imp=0, p75_lag_2_prec_median_imp=0, p25_min_temp_mean_lin=0, p75_min_temp_mean_lin=0,
p25_max_temp_mean_lin=0, p75_max_temp_mean_lin=0, p25_total_precip_median_lin=0, p75_total_precip_median_lin=0))|> data.frame() |>
mutate(treated_cntrl= ifelse(id==46, "treated","control"))After the second meeting with Ignacio, we decided to compare exclusivelly one week with the race. So we selected thursday to sunday dates.
collisions_weather_corr_rect_pre_post <-
collisions_weather_corr_rect |>
#4 days , 2 weeks before
mutate(
pre_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - days((14+3)),
race_date_rec - days((14-0)))), 0L),
exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - days(3),
race_date_rec + days(0))), 0L),
#4 days , 2 weeks after
neg_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec + days(14-3),
race_date_rec + days(14+0))), 0L),
exposure_window= case_when(pre_exposure_post_conv==1~ "Pre-exposure",
exposure_post_conv==1~ "Exposure",
neg_exposure_post_conv==1~ "Post-exposure"))|>
group_by(cluster, exposure_window)|>
summarise(
sum_velocidad = sum(velocidad, na.rm = TRUE),
sum_velocidad_lead1= sum(velocidad_lead1_post_conv_imp, na.rm = T),
sum_velocidad_lead2= sum(velocidad_lead2_post_conv_imp, na.rm = T),
exp_veh = sum(vehicles_use_type, na.rm = TRUE), # Nº de vehículos expuestos
exp_lic = sum(license_holders_sex_age, na.rm = TRUE), # Nº de conductores
median_lag_2_prec_median_imp = median(lag_2_prec_median_lin_imp, na.rm = TRUE),
mean_min_temp_mean_lin = mean(min_temp_mean_lin, na.rm = TRUE),
sd_min_temp_mean_lin = sd(min_temp_mean_lin, na.rm = TRUE),
mean_max_temp_mean_lin = mean(max_temp_mean_lin, na.rm = TRUE),
sd_max_temp_mean_lin = sd(max_temp_mean_lin, na.rm = TRUE),
median_total_precip_median_lin = median(total_precip_median_lin, na.rm = TRUE),
q1_total_precip_median_lin = quantile(total_precip_median_lin, .25, na.rm = TRUE),
q3_total_precip_median_lin = quantile(total_precip_median_lin, .75, na.rm = TRUE),
min_date= as.Date(min(unclass(date))),
max_date= as.Date(max(unclass(date))),
.groups = "drop"
) |>
mutate( # ── aux variables
rate_veh = (sum_velocidad / exp_veh) * 1e6, # not to model
rate_lic = (sum_velocidad / exp_lic) * 1e6, # not to model
rate_lead1_veh = (sum_velocidad_lead1 / exp_veh) * 1e6, # not to model
rate_lead1_lic = (sum_velocidad_lead1 / exp_lic) * 1e6, # not to model
rate_lead2_veh = (sum_velocidad_lead2 / exp_veh) * 1e6, # not to model
rate_lead2_lic = (sum_velocidad_lead2 / exp_lic) * 1e6, # not to model
off_veh = log(exp_veh), # offset
off_lic = log(exp_lic)
) |>
filter(!is.na(exposure_window)) |>
mutate(
treated = !grepl("\\.(23)$", cluster) &
grepl(paste0("^(", paste(as.character(setdiff(2000:2022, c(2009, 2020, 2021))), collapse="|"), ")"), cluster)
) |>
mutate(
treated_sens = !grepl("\\.(23|43)$", cluster) &
grepl(paste0("^(", paste(as.character(setdiff(2000:2022, c(2009, 2020, 2021))), collapse="|"), ")"), cluster),
D= paste0(exposure_window, ".", treated),
D_sens= paste0(exposure_window, ".", treated_sens)
)
# collisions_weather_corr_rect_pre_post |>
# filter(pre_exposure_post_conv==1) |>
# group_by(cluster) |>
# summarise(min=min(date), max= max(date))
collisions_weather_corr_rect_pre_post_full <-
collisions_weather_corr_rect |>
#4 days , 2 weeks before
mutate(
pre_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - days((14+3)),
race_date_rec - days((14-0)))), 0L),
exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - days(3),
race_date_rec + days(0))), 0L),
#4 days , 2 weeks after
neg_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec + days(14-3),
race_date_rec + days(14+0))), 0L),
exposure_window= case_when(pre_exposure_post_conv==1~ "Pre-exposure",
exposure_post_conv==1~ "Exposure",
neg_exposure_post_conv==1~ "Post-exposure"))|>
filter(!is.na(exposure_window))|>
mutate(
treated = !grepl("\\.(23)$", cluster) &
grepl(paste0("^(", paste(as.character(setdiff(2000:2022, c(2009, 2020, 2021))), collapse="|"), ")"), cluster)
) |>
mutate(
treated_sens = !grepl("\\.(23|43)$", cluster) &
grepl(paste0("^(", paste(as.character(setdiff(2000:2022, c(2009, 2020, 2021))), collapse="|"), ")"), cluster),
D= paste0(exposure_window, ".", treated),
D_sens= paste0(exposure_window, ".", treated_sens)
)
rbind.data.frame(
head(collisions_weather_corr_rect_pre_post_full |> dplyr::select(yday_corr, date, race_date_rec, pre_exposure_post_conv, exposure_post_conv, neg_exposure_post_conv) |> mutate(wd=weekdays(date)), 13),
tail(collisions_weather_corr_rect_pre_post_full |> dplyr::select(yday_corr, date, race_date_rec, pre_exposure_post_conv, exposure_post_conv, neg_exposure_post_conv) |> mutate(wd=weekdays(date)), 13)
)|> knitr::kable("markdown", caption= "First and last 13 days of the database with different exposure windows")
# ── 1. load package ------------------------------------------
library(tableone)
# ── 2. identify the variables by distribution -----------------
nonnormal <- c("rate_veh",
"rate_lic",
"rate_lead1_veh",
"rate_lead1_lic",
"rate_lead2_veh",
"rate_lead2_lic",
"median_lag_2_prec_median_imp",
"median_total_precip_median_lin",
"q1_total_precip_median_lin",
"q3_total_precip_median_lin")
normal <- c("mean_min_temp_mean_lin",
"sd_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"sd_max_temp_mean_lin")
vars <- c(nonnormal, normal)
# ── 4. create Table 1 -----------------------------------------
#Exposure.FALSE Exposure.TRUE Post-exposure.FALSE Post-exposure.TRUE Pre-exposure.FALSE Pre-exposure.TRUE
tab1 <- CreateTableOne(vars = vars,
strata = "D",
data = subset(collisions_weather_corr_rect_pre_post, grepl("Pre",exposure_window)),
test = TRUE)
# ── 5. print / export -----------------------------------------
print(tab1, smd = TRUE, nonnormal= nonnormal, quote = TRUE, noSpaces = TRUE) # console
# write.csv(print(tab1), "table1.csv") # if you want a CSV
tab1_df <- print(
x = tab1,
smd = TRUE,
nonnormal= nonnormal,
quote = FALSE, # sin comillas
noSpaces = TRUE, # sin espacios extra
printToggle = FALSE # devuelve lo que imprime en vez de cat()
) |>
as.data.frame(stringsAsFactors = FALSE) |>
rownames_to_column(var = "Term") # pasa rownames a columna
names(tab1_df[c(1:4,6)]) <- c("Term", "Non-host MRC-years", "Host MRCs-years", "Sig.", "SMD")
tab1_df[,4:5]<-NULL
tab1_sens <- CreateTableOne(vars = vars,
strata = "D_sens",
data = subset(collisions_weather_corr_rect_pre_post, grepl("Pre",exposure_window)),
test = TRUE)
tab1_sens_df <- print(
x = tab1_sens,
smd = TRUE,
nonnormal= nonnormal,
quote = FALSE, # sin comillas
noSpaces = TRUE, # sin espacios extra
printToggle = FALSE # devuelve lo que imprime en vez de cat()
) |>
as.data.frame(stringsAsFactors = FALSE) |>
rownames_to_column(var = "Term") # pasa rownames a columna
names(tab1_sens_df[c(1:4,6)]) <- c("Term", "Non-host MRC-years", "Host MRCs-years", "Sig.", "SMD")
tab1_sens_df[,4:5]<-NULL
cbind.data.frame(type=c("Controls: Quebec only", rep("", 14),"Controls: Quebec+Sherbrooke", rep("", 14)),
rbind.data.frame(
tab1_df,
tab1_sens_df
)
) |>
knitr::kable("markdown", caption="Pre-exposure, MRCs-years affected by the race vs. non-host or in years 2009,2020, 2021")
# cbind.data.frame(type=c("Controls: Quebec only", rep("", 14),"Controls: Quebec+Sherbrooke", rep("", 14)),
# rbind.data.frame(
# tab1_df,
# tab1_sens_df
# )
# ) |> filter(!grepl("_veh",Term)) |> rio::export("clipboard")| yday_corr | date | race_date_rec | pre_exposure_post_conv | exposure_post_conv | neg_exposure_post_conv | wd |
|---|---|---|---|---|---|---|
| 21 | 2000-06-01 | 2000-06-18 | 1 | 0 | 0 | jueves |
| 22 | 2000-06-02 | 2000-06-18 | 1 | 0 | 0 | viernes |
| 23 | 2000-06-03 | 2000-06-18 | 1 | 0 | 0 | sábado |
| 24 | 2000-06-04 | 2000-06-18 | 1 | 0 | 0 | domingo |
| 35 | 2000-06-15 | 2000-06-18 | 0 | 1 | 0 | jueves |
| 36 | 2000-06-16 | 2000-06-18 | 0 | 1 | 0 | viernes |
| 37 | 2000-06-17 | 2000-06-18 | 0 | 1 | 0 | sábado |
| 38 | 2000-06-18 | 2000-06-18 | 0 | 1 | 0 | domingo |
| 49 | 2000-06-29 | 2000-06-18 | 0 | 0 | 1 | jueves |
| 50 | 2000-06-30 | 2000-06-18 | 0 | 0 | 1 | viernes |
| 51 | 2000-07-01 | 2000-06-18 | 0 | 0 | 1 | sábado |
| 52 | 2000-07-02 | 2000-06-18 | 0 | 0 | 1 | domingo |
| 21 | 2001-05-24 | 2001-06-10 | 1 | 0 | 0 | jueves |
| 52 | 2021-06-27 | 2021-06-13 | 0 | 0 | 1 | domingo |
| 21 | 2022-06-02 | 2022-06-19 | 1 | 0 | 0 | jueves |
| 22 | 2022-06-03 | 2022-06-19 | 1 | 0 | 0 | viernes |
| 23 | 2022-06-04 | 2022-06-19 | 1 | 0 | 0 | sábado |
| 24 | 2022-06-05 | 2022-06-19 | 1 | 0 | 0 | domingo |
| 35 | 2022-06-16 | 2022-06-19 | 0 | 1 | 0 | jueves |
| 36 | 2022-06-17 | 2022-06-19 | 0 | 1 | 0 | viernes |
| 37 | 2022-06-18 | 2022-06-19 | 0 | 1 | 0 | sábado |
| 38 | 2022-06-19 | 2022-06-19 | 0 | 1 | 0 | domingo |
| 49 | 2022-06-30 | 2022-06-19 | 0 | 0 | 1 | jueves |
| 50 | 2022-07-01 | 2022-06-19 | 0 | 0 | 1 | viernes |
| 51 | 2022-07-02 | 2022-06-19 | 0 | 0 | 1 | sábado |
| 52 | 2022-07-03 | 2022-06-19 | 0 | 0 | 1 | domingo |
"Stratified by D"
"" "Pre-exposure.FALSE"
"n" "35"
"rate_veh (median [IQR])" "0.76 [0.39, 1.66]"
"rate_lic (median [IQR])" "0.70 [0.47, 1.62]"
"rate_lead1_veh (median [IQR])" "0.54 [0.41, 1.83]"
"rate_lead1_lic (median [IQR])" "0.70 [0.49, 1.95]"
"rate_lead2_veh (median [IQR])" "0.62 [0.40, 1.81]"
"rate_lead2_lic (median [IQR])" "0.74 [0.48, 1.79]"
"median_lag_2_prec_median_imp (median [IQR])" "0.25 [0.01, 3.62]"
"median_total_precip_median_lin (median [IQR])" "0.60 [0.00, 4.09]"
"q1_total_precip_median_lin (median [IQR])" "0.00 [0.00, 1.98]"
"q3_total_precip_median_lin (median [IQR])" "2.22 [1.54, 7.74]"
"mean_min_temp_mean_lin (mean (SD))" "8.28 (3.47)"
"sd_min_temp_mean_lin (mean (SD))" "3.08 (2.29)"
"mean_max_temp_mean_lin (mean (SD))" "18.70 (3.48)"
"sd_max_temp_mean_lin (mean (SD))" "3.99 (2.51)"
"Stratified by D"
"" "Pre-exposure.TRUE" "p"
"n" "80" ""
"rate_veh (median [IQR])" "0.92 [0.00, 2.04]" "0.548"
"rate_lic (median [IQR])" "1.03 [0.00, 2.08]" "0.681"
"rate_lead1_veh (median [IQR])" "0.88 [0.15, 1.83]" "0.956"
"rate_lead1_lic (median [IQR])" "0.96 [0.18, 1.77]" "0.840"
"rate_lead2_veh (median [IQR])" "0.77 [0.00, 1.80]" "0.740"
"rate_lead2_lic (median [IQR])" "0.80 [0.00, 1.84]" "0.645"
"median_lag_2_prec_median_imp (median [IQR])" "0.53 [0.05, 2.00]" "0.958"
"median_total_precip_median_lin (median [IQR])" "0.66 [0.22, 2.45]" "0.552"
"q1_total_precip_median_lin (median [IQR])" "0.08 [0.00, 1.02]" "0.493"
"q3_total_precip_median_lin (median [IQR])" "2.97 [0.83, 5.07]" "0.582"
"mean_min_temp_mean_lin (mean (SD))" "10.11 (2.75)" "0.003"
"sd_min_temp_mean_lin (mean (SD))" "2.51 (1.18)" "0.083"
"mean_max_temp_mean_lin (mean (SD))" "20.78 (3.48)" "0.004"
"sd_max_temp_mean_lin (mean (SD))" "3.56 (1.42)" "0.247"
"Stratified by D"
"" "test" "SMD"
"n" "" ""
"rate_veh (median [IQR])" "nonnorm" "0.197"
"rate_lic (median [IQR])" "nonnorm" "0.130"
"rate_lead1_veh (median [IQR])" "nonnorm" "0.032"
"rate_lead1_lic (median [IQR])" "nonnorm" "0.040"
"rate_lead2_veh (median [IQR])" "nonnorm" "0.008"
"rate_lead2_lic (median [IQR])" "nonnorm" "0.057"
"median_lag_2_prec_median_imp (median [IQR])" "nonnorm" "0.119"
"median_total_precip_median_lin (median [IQR])" "nonnorm" "0.032"
"q1_total_precip_median_lin (median [IQR])" "nonnorm" "0.006"
"q3_total_precip_median_lin (median [IQR])" "nonnorm" "0.026"
"mean_min_temp_mean_lin (mean (SD))" "" "0.585"
"sd_min_temp_mean_lin (mean (SD))" "" "0.311"
"mean_max_temp_mean_lin (mean (SD))" "" "0.597"
"sd_max_temp_mean_lin (mean (SD))" "" "0.210"| type | Term | Pre-exposure.FALSE | Pre-exposure.TRUE | SMD |
|---|---|---|---|---|
| Controls: Quebec only | n | 35 | 80 | |
| rate_veh (median [IQR]) | 0.76 [0.39, 1.66] | 0.92 [0.00, 2.04] | 0.197 | |
| rate_lic (median [IQR]) | 0.70 [0.47, 1.62] | 1.03 [0.00, 2.08] | 0.130 | |
| rate_lead1_veh (median [IQR]) | 0.54 [0.41, 1.83] | 0.88 [0.15, 1.83] | 0.032 | |
| rate_lead1_lic (median [IQR]) | 0.70 [0.49, 1.95] | 0.96 [0.18, 1.77] | 0.040 | |
| rate_lead2_veh (median [IQR]) | 0.62 [0.40, 1.81] | 0.77 [0.00, 1.80] | 0.008 | |
| rate_lead2_lic (median [IQR]) | 0.74 [0.48, 1.79] | 0.80 [0.00, 1.84] | 0.057 | |
| median_lag_2_prec_median_imp (median [IQR]) | 0.25 [0.01, 3.62] | 0.53 [0.05, 2.00] | 0.119 | |
| median_total_precip_median_lin (median [IQR]) | 0.60 [0.00, 4.09] | 0.66 [0.22, 2.45] | 0.032 | |
| q1_total_precip_median_lin (median [IQR]) | 0.00 [0.00, 1.98] | 0.08 [0.00, 1.02] | 0.006 | |
| q3_total_precip_median_lin (median [IQR]) | 2.22 [1.54, 7.74] | 2.97 [0.83, 5.07] | 0.026 | |
| mean_min_temp_mean_lin (mean (SD)) | 8.28 (3.47) | 10.11 (2.75) | 0.585 | |
| sd_min_temp_mean_lin (mean (SD)) | 3.08 (2.29) | 2.51 (1.18) | 0.311 | |
| mean_max_temp_mean_lin (mean (SD)) | 18.70 (3.48) | 20.78 (3.48) | 0.597 | |
| sd_max_temp_mean_lin (mean (SD)) | 3.99 (2.51) | 3.56 (1.42) | 0.210 | |
| Controls: Quebec+Sherbrooke | n | 55 | 60 | |
| rate_veh (median [IQR]) | 0.42 [0.00, 1.14] | 1.34 [0.51, 2.45] | 0.688 | |
| rate_lic (median [IQR]) | 0.50 [0.00, 1.18] | 1.23 [0.54, 2.24] | 0.599 | |
| rate_lead1_veh (median [IQR]) | 0.44 [0.00, 1.25] | 1.01 [0.51, 2.08] | 0.502 | |
| rate_lead1_lic (median [IQR]) | 0.51 [0.00, 1.29] | 1.03 [0.54, 2.10] | 0.404 | |
| rate_lead2_veh (median [IQR]) | 0.42 [0.00, 1.55] | 1.09 [0.29, 1.96] | 0.463 | |
| rate_lead2_lic (median [IQR]) | 0.50 [0.00, 1.55] | 1.03 [0.29, 1.89] | 0.373 | |
| median_lag_2_prec_median_imp (median [IQR]) | 0.53 [0.01, 4.03] | 0.38 [0.05, 1.38] | 0.320 | |
| median_total_precip_median_lin (median [IQR]) | 0.80 [0.04, 3.88] | 0.64 [0.12, 2.02] | 0.237 | |
| q1_total_precip_median_lin (median [IQR]) | 0.11 [0.00, 2.21] | 0.05 [0.00, 0.74] | 0.224 | |
| q3_total_precip_median_lin (median [IQR]) | 3.18 [1.54, 7.74] | 2.85 [0.76, 5.03] | 0.294 | |
| mean_min_temp_mean_lin (mean (SD)) | 8.27 (3.19) | 10.72 (2.51) | 0.854 | |
| sd_min_temp_mean_lin (mean (SD)) | 3.01 (1.96) | 2.38 (1.14) | 0.396 | |
| mean_max_temp_mean_lin (mean (SD)) | 19.17 (3.62) | 21.04 (3.36) | 0.535 | |
| sd_max_temp_mean_lin (mean (SD)) | 3.86 (2.20) | 3.54 (1.39) | 0.170 |
acf(collisions_weather_corr_rect_synth_tr$rate_veh, lag.max = 50, main = "Treatments\nACF:high-speed tickets per 1MM vehicles")
pacf(collisions_weather_corr_rect_synth_tr$rate_veh, lag.max = 50, main = "Treatments\nPACF:high-speed tickets per 1MM vehicles")
We noticed a partial autocorrelation of 10 days in the treated.
acf(collisions_weather_corr_rect_synth_cntr$rate_veh*10, lag.max = 50, main = "Controls,\nACF:high-speed tickets per 1MM vehicles")
pacf(collisions_weather_corr_rect_synth_cntr$rate_veh*10, lag.max = 50, main = "Controls,\nPACF:high-speed tickets per 1MM vehicles")
We noticed a partial autocorrelation of 25 days in the controls.
x_vals <- time(forecast::msts(
collisions_weather_corr_rect_synth_cntr$rate_veh,
start = c(2019, 1),
seasonal.periods = c(7)
))
forecast::autoplot(forecast::mstl(forecast::msts(collisions_weather_corr_rect_synth_tr$rate_veh*10,start = c(2019, 1),seasonal.periods = c(7)), robust = TRUE))+ ggtitle("Treated (per10M)")+
scale_x_continuous(
breaks = seq(2019, 2025, by = 1/7)[seq(1, 43, by = 7)], # where the ticks sit
labels = seq(-35, 7)[seq(1, 43, by = 7)] # what the ticks say
) +
geom_vline(
xintercept = seq(2019, 2025, by = 1/7)[which.min(abs(seq(-abs(min_value)+1, 7)))], # rel_day 0
linetype = "dashed", colour = "black")
forecast::autoplot(forecast::mstl(forecast::msts(collisions_weather_corr_rect_synth_cntr$rate_veh*10,start = c(2019, 1),seasonal.periods = c(7)), robust = TRUE))+ ggtitle("Controls (per 10M)")+
scale_x_continuous(
breaks = seq(2019, 2025, by = 1/7)[seq(1, 43, by = 7)], # where the ticks sit
labels = seq(-35, 7)[seq(1, 43, by = 7)] # what the ticks say
) +
geom_vline(
xintercept = seq(2019, 2025, by = 1/7)[which.min(abs(seq(-abs(min_value)+1, 7)))], # rel_day 0
linetype = "dashed", colour = "black")
The decomposition shows that differences between treated and control groups emerge in the trend, not in weekly seasonality, with treated show a downward trend, while controls show an upward trend. Any modelling or graphical comparison should therefore adjust for the weekly effect and focus on the diverging long-term trajectory.
We visualize series differentiating by treated or control cities.
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = sum_velocidad)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = abs(min_value)+1, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Number of High‑speed Tickets (sum)") +
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = sum_velocidad)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = abs(min_value)+1, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Number of High‑speed Tickets (sum)") +
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = (sum_velocidad/exp_veh)*1e6)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Number of Daily high-speed collisions per 1MM vehicle counts (sum)") +
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = (sum_velocidad/exp_veh)*1e6)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Number of Daily high-speed collisions per 1MM vehicle counts (sum)") +
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_min_temp_mean_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Min temperature") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_min_temp_mean_lin,
ymin = p25_min_temp_mean_lin,
ymax = p75_min_temp_mean_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = mean_min_temp_mean_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Min temperature") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_min_temp_mean_lin,
ymin = p25_min_temp_mean_lin,
ymax = p75_min_temp_mean_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_max_temp_mean_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Max temperature") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_max_temp_mean_lin,
ymin = p25_max_temp_mean_lin,
ymax = p75_max_temp_mean_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = mean_max_temp_mean_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Max temperature") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_max_temp_mean_lin,
ymin = p25_max_temp_mean_lin,
ymax = p75_max_temp_mean_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_median_total_precip_median_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Precipitations \n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_median_total_precip_median_lin,
ymin = p25_total_precip_median_lin,
ymax = p75_total_precip_median_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = mean_median_total_precip_median_lin)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Precipitations \n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_median_total_precip_median_lin,
ymin = p25_total_precip_median_lin,
ymax = p75_total_precip_median_lin),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
ggplot(collisions_weather_corr_rect_synth_cntr_tr, aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 30, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-30 to 5 days post-treatment)", y = "Precipitations (2-day lagged)\n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +
geom_ribbon(
data = collisions_weather_corr_rect_synth_cntr,
aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp,
ymin = p25_lag_2_prec_median_imp,
ymax = p75_lag_2_prec_median_imp),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom")
ggplot(collisions_weather_corr_rect_quebec_synth_cntr_tr, aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp)) +
geom_line(aes(group = id, colour = treated_cntrl)) +
scale_colour_brewer(palette = "Dark2", name = "MRC") +
geom_vline(xintercept = 38, linetype = "dashed", color = "black", linewidth=1) +
labs(x = "Days (-38 to 5 days post-treatment)", y = "Precipitations (2-day lagged)\n(annual mean of city-level medians, where each city median\nis computed from the medians of its weather stations)") +
geom_ribbon(
data = collisions_weather_corr_rect_quebec_synth_cntr,
aes(x = yday_corr, y = mean_median_lag_2_prec_median_imp,
ymin = p25_lag_2_prec_median_imp,
ymax = p75_lag_2_prec_median_imp),
inherit.aes = FALSE,
fill = "grey70", # o un color que combine con tu paleta
alpha = 0.3
)+
theme_minimal() +
theme(panel.grid.minor = element_blank(), legend.position = "bottom") 
library(DiagrammeR)
gr <-
grViz("
digraph study_design {
graph [layout = dot, rankdir = TB]
node [shape = rectangle, style = filled, color = LightSkyBlue, fontname = Helvetica]
start [label = 'Analysis Start']
adelanto [label = 'Advance two days (D)']
no_adelanto [label = 'No advance\n(from race day) (D_off)']
quebec_only [label = 'Quebec only as control (D)']
quebec_sherbrooke [label = 'Quebec and Sherbrooke\nas control (D_sens)']
exp_3days [label = 'Exposure 3 days\nafter the race (D)']
exp_7days [label = 'Exposure 7 days\nafter the race (D7)']
# edges
start -> {adelanto no_adelanto}
adelanto -> {quebec_only quebec_sherbrooke}
quebec_only -> {exp_3days exp_7days}
}",
width = 1000,
height = 1400)
grunlink(paste0(getwd(),"/_figs/analysisplan_files"), recursive = TRUE)
htmlwidgets::saveWidget(gr, paste0(getwd(),"/_figs/analysisplan.html"))
webshot::webshot(paste0(getwd(),"/_figs/analysisplan.html"),
paste0(getwd(),"/_figs/analysisplan.png"),
vwidth = 300, vheight = 300*1.2, zoom=10, expand=100) # Test w/ diff coords:top, left, width, and height
cat(paste0("R library: ", Sys.getenv("R_LIBS_USER")))
cat(paste0("Date: ",withr::with_locale(new = c('LC_TIME' = 'C'), code =Sys.time())))
cat(paste0("Editor context: ", getwd()))
cat("quarto version: "); system("quarto --version")
quarto::quarto_version()
save.image("_data/step12.RData")
saveRDS(collisions_weather_corr, file = "_data/collisions_weather.rds", ascii = FALSE, version = NULL, compress = TRUE, refhook = NULL)R library: H:/My Drive/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1/renv/library/windows/R-4.4/x86_64-w64-mingw32Date: 2025-07-29 17:22:23.34514Editor context: H:/My Drive/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1quarto version: [1] 0
[1] '1.7.29'sesion_info <- devtools::session_info()Warning in system2("quarto", "-V", stdout = TRUE, env = paste0("TMPDIR=", : el
comando ejecutado '"quarto"
TMPDIR=C:/Users/andre/AppData/Local/Temp/RtmpmMdxFK/file16a4469b1488d -V' tiene
el estatus 1dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
) |>
knitr::kable(caption = "R packages", format = "html",
col.names = c("Row number", "Package", "Version"),
row.names = FALSE,
align = c("c", "l", "r")) |>
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|>
kableExtra::scroll_box(width = "100%", height = "375px") | Row number | Package | Version |
|---|---|---|
| abind | 1.4-8 | RSPM |
| assertthat | 0.2.1 | RSPM |
| backports | 1.5.0 | RSPM |
| bayesplot | 1.12.0 | RSPM |
| bdsmatrix | 1.3-7 | RSPM |
| bit | 4.6.0 | CRAN (R 4.4.3) |
| bit64 | 4.6.0-1 | CRAN (R 4.4.3) |
| Boom | 0.9.15 | RSPM |
| BoomSpikeSlab | 1.2.6 | RSPM |
| boot | 1.3-30 | CRAN (R 4.4.1) |
| bpmn | 0.1.0 | Github (bergant/bpmn@628d3efaa27544c221b2fd7f1895301a63b70c49) |
| bridgesampling | 1.1-2 | RSPM |
| brms | 2.22.0 | RSPM |
| Brobdingnag | 1.2-9 | RSPM |
| broom | 1.0.8 | RSPM |
| bsts | 0.9.10 | RSPM |
| cachem | 1.1.0 | CRAN (R 4.4.3) |
| callr | 3.7.6 | RSPM |
| car | 3.1-3 | RSPM |
| carData | 3.0-5 | RSPM |
| CausalImpact | 1.3.0 | RSPM |
| CBPS | 0.23 | RSPM |
| checkmate | 2.3.2 | RSPM |
| class | 7.3-22 | CRAN (R 4.4.1) |
| cli | 3.6.5 | RSPM |
| coda | 0.19-4.1 | RSPM |
| codetools | 0.2-20 | CRAN (R 4.4.1) |
| collapse | 2.1.2 | RSPM |
| colorspace | 2.1-1 | RSPM |
| conquer | 1.3.3 | RSPM |
| cubature | 2.1.4 | RSPM |
| curl | 6.2.3 | CRAN (R 4.4.1) |
| CVXR | 1.0-15 | RSPM |
| dagitty | 0.3-4 | RSPM |
| data.table | 1.17.4 | RSPM |
| DBI | 1.2.3 | RSPM |
| devtools | 2.4.5 | RSPM |
| DHARMa | 0.4.7 | RSPM |
| DiagrammeR | 1.0.11 | RSPM |
| digest | 0.6.37 | RSPM |
| distributional | 0.5.0 | RSPM |
| doParallel | 1.0.17 | RSPM |
| doRNG | 1.8.6.2 | RSPM |
| doSNOW | 1.0.20 | RSPM |
| dplyr | 1.1.4 | RSPM |
| dreamerr | 1.5.0 | RSPM |
| e1071 | 1.7-16 | RSPM |
| ECOSolveR | 0.5.5 | RSPM |
| ellipsis | 0.3.2 | RSPM |
| emmeans | 1.11.1 | RSPM |
| estimability | 1.5.1 | RSPM |
| evaluate | 1.0.3 | RSPM |
| farver | 2.1.2 | RSPM |
| fastDummies | 1.7.5 | RSPM |
| fastmap | 1.2.0 | CRAN (R 4.4.3) |
| fect | 1.0.0 | RSPM |
| fixest | 0.12.1 | RSPM |
| forcats | 1.0.0 | RSPM |
| foreach | 1.5.2 | RSPM |
| forecast | 8.24.0 | RSPM |
| Formula | 1.2-5 | RSPM |
| fracdiff | 1.5-3 | RSPM |
| fs | 1.6.6 | RSPM |
| future | 1.49.0 | RSPM |
| future.apply | 1.11.3 | RSPM |
| geeM | 0.10.1 | RSPM |
| geepack | 1.3.12 | RSPM |
| generics | 0.1.4 | RSPM |
| geosphere | 1.5-20 | RSPM |
| GGally | 2.2.1 | RSPM |
| ggdag | 0.2.13 | RSPM |
| ggplot2 | 3.5.2 | RSPM |
| ggstats | 0.9.0 | RSPM |
| glmmTMB | 1.1.11 | RSPM |
| glmnet | 4.1-9 | RSPM |
| glmx | 0.2-1 | RSPM |
| globals | 0.18.0 | RSPM |
| glue | 1.8.0 | RSPM |
| gmp | 0.7-5 | RSPM |
| gnm | 1.1-5 | RSPM |
| gridExtra | 2.3 | RSPM |
| gtable | 0.3.6 | RSPM |
| gtools | 3.9.5 | RSPM |
| haven | 2.5.4 | RSPM |
| hms | 1.1.3 | CRAN (R 4.4.3) |
| htmltools | 0.5.8.1 | RSPM |
| htmlwidgets | 1.6.4 | RSPM |
| httpuv | 1.6.16 | RSPM |
| httr2 | 1.1.2 | RSPM |
| igraph | 2.1.4 | RSPM |
| iterators | 1.0.14 | RSPM |
| jsonlite | 2.0.0 | CRAN (R 4.4.3) |
| kableExtra | 1.4.0 | RSPM |
| kernlab | 0.9-33 | RSPM |
| knitr | 1.50 | RSPM |
| labeling | 0.4.3 | RSPM |
| labelled | 2.14.1 | RSPM |
| later | 1.4.2 | RSPM |
| lattice | 0.22-6 | CRAN (R 4.4.1) |
| lfe | 3.1.1 | RSPM |
| lifecycle | 1.0.4 | RSPM |
| listenv | 0.9.1 | RSPM |
| lme4 | 1.1-37 | RSPM |
| lmtest | 0.9-40 | RSPM |
| loo | 2.8.0 | RSPM |
| lubridate | 1.9.4 | RSPM |
| magrittr | 2.0.3 | RSPM |
| MASS | 7.3-60.2 | CRAN (R 4.4.1) |
| MatchIt | 4.7.2 | RSPM |
| Matrix | 1.7-0 | CRAN (R 4.4.1) |
| MatrixModels | 0.5-4 | RSPM |
| matrixStats | 1.5.0 | RSPM |
| maxLik | 1.5-2.1 | RSPM |
| memoise | 2.0.1 | CRAN (R 4.4.3) |
| mgcv | 1.9-1 | CRAN (R 4.4.1) |
| mime | 0.13 | CRAN (R 4.4.3) |
| miniUI | 0.1.2 | RSPM |
| minqa | 1.2.8 | RSPM |
| miscTools | 0.6-28 | RSPM |
| mitools | 2.4 | RSPM |
| mvtnorm | 1.3-3 | RSPM |
| nixtlar | 0.6.2 | RSPM |
| nlme | 3.1-164 | CRAN (R 4.4.1) |
| nloptr | 2.2.1 | RSPM |
| nnet | 7.3-19 | CRAN (R 4.4.1) |
| np | 0.60-18 | RSPM |
| numDeriv | 2016.8-1.1 | RSPM |
| openxlsx | 4.2.8 | RSPM |
| optimx | 2025-4.9 | RSPM |
| PanelMatch | 3.1.1 | RSPM |
| parallelly | 1.44.0 | RSPM |
| pillar | 1.10.2 | RSPM |
| pkgbuild | 1.4.8 | RSPM |
| pkgconfig | 2.0.3 | RSPM |
| pkgload | 1.4.0 | RSPM |
| plm | 2.6-6 | RSPM |
| plyr | 1.8.9 | RSPM |
| posterior | 1.6.1 | RSPM |
| pracma | 2.4.4 | CRAN (R 4.4.1) |
| processx | 3.8.6 | RSPM |
| profvis | 0.4.0 | RSPM |
| promises | 1.3.2 | RSPM |
| proxy | 0.4-27 | RSPM |
| ps | 1.9.1 | RSPM |
| purrr | 1.0.4 | RSPM |
| Qtools | 1.5.9 | RSPM |
| quadprog | 1.5-8 | RSPM |
| quantdr | 1.2.2 | RSPM |
| quantmod | 0.4.28 | RSPM |
| quantreg | 6.1 | RSPM |
| quarto | 1.4.4 | RSPM |
| qvcalc | 1.0.4 | RSPM |
| R6 | 2.6.1 | RSPM |
| rappdirs | 0.3.3 | CRAN (R 4.4.3) |
| rbibutils | 2.3 | RSPM |
| RColorBrewer | 1.1-3 | RSPM |
| Rcpp | 1.0.14 | RSPM |
| RcppParallel | 5.1.10 | RSPM |
| Rdpack | 2.6.4 | RSPM |
| readr | 2.1.5 | CRAN (R 4.4.3) |
| reformulas | 0.4.1 | RSPM |
| relimp | 1.0-5 | RSPM |
| remotes | 2.5.0 | RSPM |
| renv | 1.1.2 | CRAN (R 4.4.1) |
| reshape2 | 1.4.4 | RSPM |
| rgenoud | 5.9-0.11 | RSPM |
| rio | 1.2.3 | RSPM |
| rlang | 1.1.6 | RSPM |
| rmarkdown | 2.29 | RSPM |
| Rmpfr | 1.1-0 | RSPM |
| rngtools | 1.5.2 | RSPM |
| rstantools | 2.4.0 | RSPM |
| rstudioapi | 0.17.1 | RSPM |
| sandwich | 3.1-1 | RSPM |
| scales | 1.4.0 | RSPM |
| scpi | 3.0.0 | RSPM |
| sessioninfo | 1.2.3 | RSPM |
| shape | 1.4.6.1 | RSPM |
| shiny | 1.10.0 | RSPM |
| snow | 0.4-4 | RSPM |
| sp | 2.2-0 | RSPM |
| SparseM | 1.84-2 | RSPM |
| stringi | 1.8.7 | RSPM |
| stringmagic | 1.2.0 | RSPM |
| stringr | 1.5.1 | RSPM |
| survey | 4.4-2 | RSPM |
| survival | 3.6-4 | CRAN (R 4.4.1) |
| svglite | 2.2.1 | RSPM |
| Synth | 1.1-8 | RSPM |
| systemfonts | 1.2.3 | RSPM |
| tableone | 0.13.2 | RSPM |
| tensorA | 0.36.2.1 | RSPM |
| textshaping | 1.0.1 | RSPM |
| tibble | 3.2.1 | RSPM |
| tidygraph | 1.3.1 | RSPM |
| tidyr | 1.3.1 | RSPM |
| tidyselect | 1.2.1 | RSPM |
| tidyverse | 2.0.0 | RSPM |
| timechange | 0.3.0 | RSPM |
| timeDate | 4041.110 | RSPM |
| TMB | 1.9.17 | RSPM |
| tseries | 0.10-58 | RSPM |
| TTR | 0.24.4 | RSPM |
| tzdb | 0.5.0 | CRAN (R 4.4.3) |
| urca | 1.3-4 | RSPM |
| urlchecker | 1.0.1 | RSPM |
| usethis | 3.1.0 | RSPM |
| V8 | 6.0.3 | RSPM |
| vctrs | 0.6.5 | RSPM |
| viridisLite | 0.4.2 | RSPM |
| visNetwork | 2.1.2 | RSPM |
| weathercan | 0.7.3.9000 | https://ropensci.r-universe.dev (R 4.4.3) |
| webshot | 0.5.5 | RSPM |
| withr | 3.0.2 | RSPM |
| xfun | 0.52 | RSPM |
| xml2 | 1.3.8 | CRAN (R 4.4.3) |
| xtable | 1.8-4 | RSPM |
| xts | 0.14.1 | RSPM |
| yaml | 2.3.10 | RSPM |
| zip | 2.3.3 | RSPM |
| zoo | 1.8-14 | RSPM |
reticulate::py_list_packages()%>%
knitr::kable(caption = "Python packages", format = "html",
col.names = c("Package", "Version", "Requirement"),
row.names = FALSE,
align = c("c", "l", "r", "r"))%>%
kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 12)|>
kableExtra::scroll_box(width = "100%", height = "375px") Error in path.expand(path): argumento 'path' inválido