Step 4. Post review (1)
# remove objects and memory
rm(list=ls());gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 839044 44.9 1661868 88.8 1126260 60.2
Vcells 1697676 13.0 8388608 64.0 3395449 26.0#remove images
while(!dev.cur())dev.off()
cat("\014")load(paste0(getwd(),"/_data/step3.RData"))#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 gaussian 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(rio); library(purrr); library(kableExtra); library(geosphere); library(geepack); library(lme4); library(glmmTMB); library(DHARMa); library(car); library(sandwich); library(emmeans); library(gnm); library(splines); library(geeM); library(plm); library(PanelMatch)
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, union
Adjuntando el paquete: 'kableExtra'The following object is masked from 'package:dplyr':
group_rowsCargando paquete requerido: Matrix
Adjuntando el paquete: 'lme4'The following object is masked from 'package:rio':
factorizeThis 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':
someThe following object is masked from 'package:dplyr':
recodeWelcome 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:stats':
weights#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
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")
)
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
compare_gnm_effects <- function(model_ref, model_cmp, comparison_name = "Ref vs Cmp",
coef_index = 1, digits = 3) {
# extract coefficients and robustly handle named/unnamed coef vectors
coef_ref <- stats::coef(model_ref)
coef_cmp <- stats::coef(model_cmp)
if (is.null(coef_ref) || is.null(coef_cmp)) stop("Both inputs must be model-like objects with coef().")
b1 <- coef_ref[[coef_index]]
b2 <- coef_cmp[[coef_index]]
# extract standard errors from vcov; allow either matrix or numeric
vc1 <- tryCatch(stats::vcov(model_ref), error = function(e) NULL)
vc2 <- tryCatch(stats::vcov(model_cmp), error = function(e) NULL)
if (is.null(vc1) || is.null(vc2)) stop("Both inputs must support vcov().")
se1 <- sqrt(diag(vc1))[[coef_index]]
se2 <- sqrt(diag(vc2))[[coef_index]]
# difference test
z <- (b2 - b1) / sqrt(se1^2 + se2^2)
p <- 2 * stats::pnorm(-abs(z))
# relative reduction ratio (RRR) and 95% CI
log_rrr <- b2 - b1
se_rrr <- sqrt(se1^2 + se2^2)
RRR <- exp(log_rrr)
lo95 <- exp(log_rrr - 1.96 * se_rrr)
hi95 <- exp(log_rrr + 1.96 * se_rrr)
# formatted output
cat(sprintf("%s — coefficient index %d\n", comparison_name, coef_index))
cat(sprintf(" z = %.*f; p = %.*f\n", digits, z, digits, p))
cat(sprintf(" log(RRR) = %.*f; se = %.*f\n", digits, log_rrr, digits, se_rrr))
cat(sprintf(" RRR = %.*f (95%% CI %.*f–%.*f)\n", digits, RRR, digits, lo95, digits, hi95))
# invisibly return numeric results for further use
invisible(list(
coef_ref = b1,
se_ref = se1,
coef_cmp = b2,
se_cmp = se2,
z = z,
p = p,
log_rrr = log_rrr,
se_rrr = se_rrr,
RRR = RRR,
lo95 = lo95,
hi95 = hi95
))
}
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#CONFIG #######################################################################
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
options(scipen=2) #display numbers rather scientific numbercollisions_weather_corr_rect_pre_post_placebo <-
collisions_weather_corr_rect |>
#4 days , 2 weeks before
mutate(
pre_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - lubridate::days((14+3)),
race_date_rec - lubridate::days((14-0)))), 0L),
exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - lubridate::days(3),
race_date_rec + lubridate::days(0))), 0L),
#4 days , 2 weeks after
neg_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec + lubridate::days(28-3),
race_date_rec + lubridate::days(28+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)
)|>
mutate(
treated_placebo = !grepl("\\.(23)$", cluster) &
grepl(paste0("^(", paste(as.character(c(2009, 2020, 2021)), collapse="|"), ")"), cluster)
) |>
mutate(
treated_sens_placebo = !grepl("\\.(23|43)$", cluster) &
grepl(paste0("^(", paste(as.character(c(2009, 2020, 2021)), collapse="|"), ")"), cluster),
D= paste0(exposure_window, ".", treated),
D_sens= paste0(exposure_window, ".", treated_sens),
D_placebo= paste0(exposure_window, ".", treated_placebo),
D_placebo_sens= paste0(exposure_window, ".", treated_sens_placebo)
)# subset(collisions_weather_corr_rect_pre_post_placebo, grepl("^Pre|^Exp",exposure_window) &
# treated_placebo==TRUE)|> mutate(cluster=as.character(cluster))|> pull(cluster)|> unique()|> dput()
cluster_placebo_treated<-
c("2009.43", "2020.43", "2021.43", "2009.58", "2020.58", "2021.58",
"2009.65", "2020.65", "2021.65", "2009.66", "2020.66", "2021.66")
# subset(collisions_weather_corr_rect_pre_post_placebo, grepl("^Pre|^Exp",exposure_window) &
# treated_placebo==FALSE) |> mutate(cluster=as.character(cluster)) |> pull(cluster) |> unique()
cluster_placebo_controls<-
c("2000.23", "2001.23", "2002.23", "2003.23", "2004.23", "2005.23",
"2006.23", "2007.23", "2008.23", "2009.23", "2010.23", "2011.23",
"2012.23", "2013.23", "2014.23", "2015.23", "2016.23", "2017.23",
"2018.23", "2019.23", "2020.23", "2021.23", "2022.23", "2000.43",
"2001.43", "2002.43", "2003.43", "2004.43", "2005.43", "2006.43",
"2007.43", "2008.43", "2010.43", "2011.43", "2012.43", "2013.43",
"2014.43", "2015.43", "2016.43", "2017.43", "2018.43", "2019.43",
"2022.43", "2000.58", "2001.58", "2002.58", "2003.58", "2004.58",
"2005.58", "2006.58", "2007.58", "2008.58", "2010.58", "2011.58",
"2012.58", "2013.58", "2014.58", "2015.58", "2016.58", "2017.58",
"2018.58", "2019.58", "2022.58", "2000.65", "2001.65", "2002.65",
"2003.65", "2004.65", "2005.65", "2006.65", "2007.65", "2008.65",
"2010.65", "2011.65", "2012.65", "2013.65", "2014.65", "2015.65",
"2016.65", "2017.65", "2018.65", "2019.65", "2022.65", "2000.66",
"2001.66", "2002.66", "2003.66", "2004.66", "2005.66", "2006.66",
"2007.66", "2008.66", "2010.66", "2011.66", "2012.66", "2013.66",
"2014.66", "2015.66", "2016.66", "2017.66", "2018.66", "2019.66",
"2022.66")
cat("*_*_*_*_*_*_*_*_*_")
cat("*_*_*_*_PLACEBO!*_*_*_*_*_")
cat("*_*_*_*_*_*_*_*_*_")
cat("Treated (w/ Sherbrooke), effects after 2 days\n")
model_gnm_tratados_lead2_PL <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) & treated_placebo==TRUE)
)
summary(model_gnm_tratados_lead2_PL)
#exp(1.12608)
#[1] 3.083545
#NS
cat("Controls (Quebec only), effects after 2 days\n")
model_gnm_controles_lead2_PL <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) & treated_placebo==FALSE)
)
summary(model_gnm_controles_lead2_PL)
#exp(0.245724)
#[1] 1.278547
#SIG
lead2_PL_main <- compare_gnm_effects(model_gnm_controles_lead2_PL,
model_gnm_tratados_lead2_PL,
comparison_name = "Treated vs. Controls (QC only; w/o Sherbrooke), placebo",
coef_index = 1, digits = 3)
# Treated vs. Controls (QC only; w/o Sherbrooke), placebo — coefficient index 1
# z = 1.212; p = 0.226
# log(RRR) = 0.880; se = 0.726
# RRR = 2.412 (95% CI 0.581–10.015)
cat("Treated (w/ Sherbrooke), effects after 2 days\n")
model_gnm_tratados_lead2_sens_PL <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) & treated_sens_placebo==TRUE)
)
summary(model_gnm_tratados_lead2_sens_PL)
# exp(0.75969)
#[1] 2.137613
#NS
cat("Controls (w/ Sherbrooke), effects after 2 days\n")
model_gnm_controles_lead2_sens_PL <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) & treated_sens_placebo==FALSE)
)
summary(model_gnm_controles_lead2_sens_PL)
#[1] 1.275913
#SIG
lead2_PL_sens <- compare_gnm_effects(model_gnm_controles_lead2_sens_PL,
model_gnm_tratados_lead2_sens_PL,
comparison_name = "Treated vs. Controls (QC + Sherbrooke), placebo",
coef_index = 1, digits = 3)
# Treated vs. Controls (QC + Sherbrooke), placebo — coefficient index 1
# z = 0.562; p = 0.574
# log(RRR) = 0.516; se = 0.918
# RRR = 1.675 (95% CI 0.277–10.128)
cat(paste0("Control years in Host-MRCs: ", paste(setdiff(paste0("20",sprintf("%02.0f",0:22)), c("2009", "2020", "2021")), collapse=", ")))
cbind.data.frame(type=c("Placebo-post (controls=QC+Sherbrooke)","","Placebo-post (controls=QC only)",""),
rbind.data.frame(
tidy_exp_gnm_corr(model_gnm_tratados_lead2_sens_PL),
tidy_exp_gnm_corr(model_gnm_controles_lead2_sens_PL),
tidy_exp_gnm_corr(model_gnm_tratados_lead2_PL),
tidy_exp_gnm_corr(model_gnm_controles_lead2_PL)))|>
mutate(RR= sprintf("%.2f (%.2f–%.2f)", `exp(Est.)`, `2.5%`, `97.5%`))|>
mutate(sig= sprintf("%.3f", P))|>
rename("p.value"="P")|>
mutate(term="Exposure")|>
dplyr::select(type, term, RR, sig)|>
knitr::kable("markdown",
caption="High-speed collisions, 2-day lagged effect", digits=3, na= "-")
cat("PLACEBO= Main analysis (sherbrooke+ quebec as controls):\n")
sprintf("RRR = %.2f (95%% CI: %.2f–%.2f), p = %s",
lead2_PL_sens$RRR,
lead2_PL_sens$lo95,
lead2_PL_sens$hi95,
ifelse(lead2_PL_sens$p < 0.001,
"<.001",
sprintf("%.2f", lead2_PL_sens$p)))
cat("PLACEBO= Secondary analysis (sherbrooke+ quebec as controls):\n")
sprintf("RRR = %.2f (95%% CI: %.2f–%.2f), p = %s",
lead2_PL_main$RRR,
lead2_PL_main$lo95,
lead2_PL_main$hi95,
ifelse(lead2_PL_main$p < 0.001,
"<.001",
sprintf("%.2f", lead2_PL_main$p)))*_*_*_*_*_*_*_*_*_*_*_*_*_PLACEBO!*_*_*_*_*_*_*_*_*_*_*_*_*_*_Treated (w/ Sherbrooke), effects after 2 days
Call:
gnm(formula = sum_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp", exposure_window) & treated_placebo ==
TRUE))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.9769824 -0.0734078 -0.0002638 0.0633610 0.5854904
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.12608
mean_min_temp_mean_lin 0.71766
mean_max_temp_mean_lin -0.86878
median_total_precip_median_lin -0.74645
median_lag_2_prec_median_imp -0.08658
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.72094
mean_min_temp_mean_lin 0.38491
mean_max_temp_mean_lin 0.43539
median_total_precip_median_lin 0.38229
median_lag_2_prec_median_imp 0.07383
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.562 0.1623
mean_min_temp_mean_lin 1.864 0.1045
mean_max_temp_mean_lin -1.995 0.0862
median_total_precip_median_lin -1.953 0.0918
median_lag_2_prec_median_imp -1.173 0.2793
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure
mean_min_temp_mean_lin
mean_max_temp_mean_lin .
median_total_precip_median_lin .
median_lag_2_prec_median_imp
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 0.1994859)
Residual deviance: 1.9391 on 7 degrees of freedom
AIC: NA
Number of iterations: 11
Controls (Quebec only), effects after 2 days
Call:
gnm(formula = sum_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp", exposure_window) & treated_placebo ==
FALSE))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.254476 -0.396152 -0.001368 0.385185 1.339121
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.245724
mean_min_temp_mean_lin 0.023783
mean_max_temp_mean_lin 0.024355
median_total_precip_median_lin -0.003734
median_lag_2_prec_median_imp 0.011276
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.088996
mean_min_temp_mean_lin 0.034686
mean_max_temp_mean_lin 0.023919
median_total_precip_median_lin 0.034355
median_lag_2_prec_median_imp 0.030340
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 2.761 0.00688
mean_min_temp_mean_lin 0.686 0.49453
mean_max_temp_mean_lin 1.018 0.31108
median_total_precip_median_lin -0.109 0.91366
median_lag_2_prec_median_imp 0.372 0.71095
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure **
mean_min_temp_mean_lin
mean_max_temp_mean_lin
median_total_precip_median_lin
median_lag_2_prec_median_imp
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 0.7742434)
Residual deviance: 92.283 on 98 degrees of freedom
AIC: NA
Number of iterations: 8
Treated vs. Controls (QC only; w/o Sherbrooke), placebo — coefficient index 1
z = 1.212; p = 0.226
log(RRR) = 0.880; se = 0.726
RRR = 2.412 (95% CI 0.581–10.015)
Treated (w/ Sherbrooke), effects after 2 days
Call:
gnm(formula = sum_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp", exposure_window) & treated_sens_placebo ==
TRUE))
Deviance Residuals:
Min 1Q Median 3Q Max
-0.9829379 -0.0420524 -0.0003016 0.0331756 0.5946266
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.75969
mean_min_temp_mean_lin 0.49547
mean_max_temp_mean_lin -0.61431
median_total_precip_median_lin -0.51702
median_lag_2_prec_median_imp -0.06850
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.91370
mean_min_temp_mean_lin 0.49590
mean_max_temp_mean_lin 0.56110
median_total_precip_median_lin 0.49357
median_lag_2_prec_median_imp 0.08949
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.831 0.452
mean_min_temp_mean_lin 0.999 0.374
mean_max_temp_mean_lin -1.095 0.335
median_total_precip_median_lin -1.048 0.354
median_lag_2_prec_median_imp -0.765 0.487
(Dispersion parameter for quasipoisson family taken to be 0.2834946)
Residual deviance: 1.5192 on 4 degrees of freedom
AIC: NA
Number of iterations: 11
Controls (w/ Sherbrooke), effects after 2 days
Call:
gnm(formula = sum_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp", exposure_window) & treated_sens_placebo ==
FALSE))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.256575 -0.388010 -0.001409 0.375711 1.332871
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.243662
mean_min_temp_mean_lin 0.023055
mean_max_temp_mean_lin 0.023860
median_total_precip_median_lin -0.004686
median_lag_2_prec_median_imp 0.011430
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.088601
mean_min_temp_mean_lin 0.034516
mean_max_temp_mean_lin 0.023802
median_total_precip_median_lin 0.034198
median_lag_2_prec_median_imp 0.030208
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 2.750 0.00706
mean_min_temp_mean_lin 0.668 0.50569
mean_max_temp_mean_lin 1.002 0.31853
median_total_precip_median_lin -0.137 0.89128
median_lag_2_prec_median_imp 0.378 0.70594
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure **
mean_min_temp_mean_lin
mean_max_temp_mean_lin
median_total_precip_median_lin
median_lag_2_prec_median_imp
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 0.7679163)
Residual deviance: 94.308 on 101 degrees of freedom
AIC: NA
Number of iterations: 8
Treated vs. Controls (QC + Sherbrooke), placebo — coefficient index 1
z = 0.562; p = 0.574
log(RRR) = 0.516; se = 0.918
RRR = 1.675 (95% CI 0.277–10.128)
Control years in Host-MRCs: 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2022| type | term | RR | sig |
|---|---|---|---|
| Placebo-post (controls=QC+Sherbrooke) | Exposure | 2.14 (0.36–12.81) | 0.406 |
| Exposure | 1.28 (1.07–1.52) | 0.006 | |
| Placebo-post (controls=QC only) | Exposure | 3.08 (0.75–12.67) | 0.118 |
| Exposure | 1.28 (1.07–1.52) | 0.006 |
PLACEBO= Main analysis (sherbrooke+ quebec as controls):
[1] "RRR = 1.68 (95% CI: 0.28–10.13), p = 0.57"
PLACEBO= Secondary analysis (sherbrooke+ quebec as controls):
[1] "RRR = 2.41 (95% CI: 0.58–10.02), p = 0.23"francos_yrs<- c("2006", "2010", "2011", "2012", "2016", "2017", "2018", "2022")
cat("Treated (w/ Sherbrooke), effects after 2 days, 4-day summary\n")
model_gnm_tratados_lead2_sinfrancos <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) &
treated==TRUE &
!grepl(paste(francos_yrs, collapse="|"),cluster))
)
summary(model_gnm_tratados_lead2_sinfrancos)
tidy_exp_gnm_corr(model_gnm_tratados_lead2_sinfrancos) |>
knitr::kable("markdown", caption="Host MRCs (without Sherbrooke) pre vs. race weekend. High-speed collisions, 2-day lagged effect")Treated (w/ Sherbrooke), effects after 2 days, 4-day summary
Call:
gnm(formula = sum_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp", exposure_window) & treated == TRUE &
!grepl(paste(francos_yrs, collapse = "|"), cluster)))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.090517 -0.387309 -0.001522 0.374840 1.379435
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.463274
mean_min_temp_mean_lin 0.048339
mean_max_temp_mean_lin -0.009404
median_total_precip_median_lin 0.014995
median_lag_2_prec_median_imp -0.024983
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.163644
mean_min_temp_mean_lin 0.049781
mean_max_temp_mean_lin 0.041369
median_total_precip_median_lin 0.064013
median_lag_2_prec_median_imp 0.044714
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 2.831 0.00703
mean_min_temp_mean_lin 0.971 0.33696
mean_max_temp_mean_lin -0.227 0.82124
median_total_precip_median_lin 0.234 0.81591
median_lag_2_prec_median_imp -0.559 0.57924
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure **
mean_min_temp_mean_lin
mean_max_temp_mean_lin
median_total_precip_median_lin
median_lag_2_prec_median_imp
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasipoisson family taken to be 0.840592)
Residual deviance: 45.886 on 43 degrees of freedom
AIC: NA
Number of iterations: 8| term | exp(Est.) | 2.5% | 97.5% | P |
|---|---|---|---|---|
| relevel(factor(exposure_window), ref = “Pre-exposure”)Exposure | 1.589269 | 1.153196 | 2.19024 | 0.0046405 |
# 1. Calculate the Ratio based on your 2025 figures
# Projected Host License Holders 2025 (Linear extrapolation of your data): ~2,538,839
# https://docs.google.com/spreadsheets/d/1znDbWVXO2A5uSKiVazQwFJKZcrrH-VZXnQsG-INBQ-w/edit?gid=819465529#gid=819465529
# influx of more than 100,000 F1 fans.
# https://gpdestinations.com/2025-canadian-grand-prix-weekend-attendance/
# Haldenby, N. (2025, June 16). **2025 Canadian Grand Prix weekend attendance exceeds 350,000**. GPDestinations. https://gpdestinations.com/2025-canadian-grand-prix-weekend-attendance/
# Attendance numbers have grown in Montreal in recent years. The last race before the coronavirus pandemic, in 2019,
# had a weekend total of 307,000 spectators. When F1 returned in 2022, the race weekend had 338,000 attendees.
# License hodlers in 2019: 2449509 (Google docs)
# Expected Visitors 2019: 350,000
#https://montreal.eater.com/2022/6/13/23166358/montreal-grand-prix-party-restaurants-bars-clubs-f1-events
visitor_ratio <- 307000 / 2449509 # Returns ~0.125
# 2. Prepare the dataset
# We create a specific dataframe for the model to ensure the offset is calculated correctly
df_model_4dsumm <- subset(collisions_weather_corr_rect_pre_post,
grepl("^Pre|^Exp", exposure_window) & treated_sens == TRUE)
# 3. Create the Adjusted Offset
# Logic: Total_Population = Locals + Visitors
# Where Visitors = Locals * visitor_ratio
# Therefore: Total_Population = Locals * (1 + visitor_ratio)
# Initialize with original license holders
df_model_4dsumm$population_adjusted <- df_model_4dsumm$exp_lic
# Identify exposure rows (Race periods)
is_exposure_4dsumm <- grepl("^Exp", df_model_4dsumm$exposure_window)
# Apply the increase ONLY to the exposure period
df_model_4dsumm$population_adjusted[is_exposure_4dsumm] <- df_model_4dsumm$population_adjusted[is_exposure_4dsumm] * (1 + visitor_ratio)cat("Treated (w/o Sherbrooke), effects after 2 days, 4-day summary\n")
model_gnm_tratados_lead2_sens_adj_offs <- gnm(
sum_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(log(population_adjusted)),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = df_model_4dsumm
)
lead2_sens_adj_offset <- compare_gnm_effects(model_gnm_controles_lead2_sens,
model_gnm_tratados_lead2_sens_adj_offs,
comparison_name = "Treated vs. Controls (QC + Sherbrooke), offset adj",
coef_index = 1, digits = 3)
# Treated vs. Controls (QC + Sherbrooke), offset adj — coefficient index 1
# z = 0.939; p = 0.348
# log(RRR) = 0.192; se = 0.204
# RRR = 1.212 (95% CI 0.812–1.809)
rbind.data.frame(
tidy_exp_gnm_corr(model_gnm_tratados_lead2_sens_adj_offs)
)|>
mutate(RR= sprintf("%.2f (%.2f–%.2f)", `exp(Est.)`, `2.5%`, `97.5%`))|>
mutate(sig= sprintf("%.3f", P))|>
rename("p.value"="P")|>
mutate(term=c("Exposure, 4-day summ"))|>
dplyr::select(term, RR, sig)|>
knitr::kable("markdown",
caption="High-speed collisions, 2-day lagged effect, offset * 1.25", digits=3, na= "-")
# |term |RR |sig |
# |:--------------------|:----------------|:-----|
# |Exposure, 4-day summ |1.22 (0.99–1.50) |0.058 |
sprintf("RRR = %.2f (95%% CI: %.2f–%.2f), p = %s",
lead2_sens_adj_offset$RRR,
lead2_sens_adj_offset$lo95,
lead2_sens_adj_offset$hi95,
ifelse(lead2_sens_adj_offset$p < 0.001,
"<.001",
sprintf("%.2f", lead2_sens_adj_offset$p)))Treated (w/o Sherbrooke), effects after 2 days, 4-day summary
Treated vs. Controls (QC + Sherbrooke), offset adj — coefficient index 1
z = 0.939; p = 0.348
log(RRR) = 0.192; se = 0.204
RRR = 1.212 (95% CI 0.812–1.809)| term | RR | sig |
|---|---|---|
| Exposure, 4-day summ | 1.22 (0.99–1.50) | 0.058 |
[1] "RRR = 1.21 (95% CI: 0.81–1.81), p = 0.35"# asegúrate de que este sea el nombre correcto en tu entorno
data_df <- collisions_weather_corr_rect_pre_post_full
mrc_names <- c("montreal", "laval", "longueuil", "quebec", "shrebrooke")
yday_pre <- c(21, 22, 23, 24)
yday_exp <- c(35, 36, 37, 38)
summarise_by_mrc <- function(df, window_pattern, treated_flag, group_name) {
df |>
dplyr::filter(grepl(window_pattern, exposure_window), treated_sens == treated_flag) |>
dplyr::mutate(
mrc = dplyr::case_when(
mrc == 66 ~ "montreal",
mrc == 65 ~ "laval",
mrc == 58 ~ "longueuil",
mrc == 23 ~ "quebec",
mrc == 43 ~ "shrebrooke",
TRUE ~ as.character(mrc)
)
) |>
dplyr::group_by(mrc) |>
dplyr::summarise(
n = dplyr::n(),
sum_vel = base::sum(velocidad, na.rm = TRUE),
mean = mean(velocidad, na.rm = TRUE),
median = stats::quantile(velocidad, .5, na.rm = TRUE),
p25 = stats::quantile(velocidad, .25, na.rm = TRUE),
p75 = stats::quantile(velocidad, .75, na.rm = TRUE),
sd = stats::sd(velocidad, na.rm = TRUE),
.groups = "drop"
) |>
tidyr::complete(mrc = mrc_names) |>
dplyr::mutate(
group = group_name,
sum_nonnorm = ifelse(is.na(median), NA_character_, sprintf("%.2f [%.2f, %.2f]", median, p25, p75)),
sum_nonnorm_n_sum = ifelse(is.na(n), NA_character_, sprintf("%s; n=%d; sum=%.2f", sum_nonnorm, n, sum_vel))
) |>
dplyr::select(mrc, group, n, sum_vel, mean, sd, median, p25, p75, sum_nonnorm, sum_nonnorm_n_sum)
}
summarise_by_yday <- function(df, window_pattern, treated_flag, group_name) {
days <- if (grepl("^Pre", window_pattern)) yday_pre else yday_exp
df |>
dplyr::filter(grepl(window_pattern, exposure_window), treated_sens == treated_flag) |>
dplyr::group_by(yday_corr) |>
dplyr::summarise(
n = dplyr::n(),
sum_vel = base::sum(velocidad, na.rm = TRUE),
mean = mean(velocidad, na.rm = TRUE),
median = stats::quantile(velocidad, .5, na.rm = TRUE),
p25 = stats::quantile(velocidad, .25, na.rm = TRUE),
p75 = stats::quantile(velocidad, .75, na.rm = TRUE),
sd = stats::sd(velocidad, na.rm = TRUE),
.groups = "drop"
) |>
tidyr::complete(yday_corr = days) |>
dplyr::mutate(
group = group_name,
sum_nonnorm = ifelse(is.na(median), NA_character_, sprintf("%.2f [%.2f, %.2f]", median, p25, p75)),
sum_nonnorm_n_sum = ifelse(is.na(n), NA_character_, sprintf("%s; n=%d; sum=%.2f", sum_nonnorm, n, sum_vel))
) |>
dplyr::select(yday_corr, group, n, sum_vel, mean, sd, median, p25, p75, sum_nonnorm, sum_nonnorm_n_sum)
}
params <- tibble::tibble(
window_pattern = c("^Pre", "^Exp", "^Pre", "^Exp"),
treated_flag = c(FALSE, FALSE, TRUE, TRUE),
group_name = c("Control_Nonhost", "GP_Nonhost", "Control_Host", "GP_Host")
)
# pmap pasando el data frame explícitamente para evitar errores de nombre
summaries_mrc_list <- purrr::pmap(
list(params$window_pattern, params$treated_flag, params$group_name),
function(window_pattern, treated_flag, group_name) {
summarise_by_mrc(data_df, window_pattern, treated_flag, group_name)
}
)
all_mrc <- dplyr::bind_rows(summaries_mrc_list)
table_mrc_wide <- all_mrc |>
tidyr::pivot_wider(
id_cols = mrc,
names_from = group,
values_from = sum_nonnorm_n_sum
) |>
dplyr::select(mrc, Control_Nonhost, GP_Nonhost, Control_Host, GP_Host)
# análogo para yday
summaries_yday_list <- purrr::pmap(
list(params$window_pattern, params$treated_flag, params$group_name),
function(window_pattern, treated_flag, group_name) {
summarise_by_yday(data_df, window_pattern, treated_flag, group_name)
}
)
all_yday <- dplyr::bind_rows(summaries_yday_list)
table_yday_wide <- all_yday |>
tidyr::pivot_wider(
id_cols = yday_corr,
names_from = group,
values_from = sum_nonnorm_n_sum
) |>
dplyr::arrange(yday_corr) |>
dplyr::select(yday_corr, Control_Nonhost, GP_Nonhost, Control_Host, GP_Host)
dplyr::bind_rows(
table_yday_wide |>
dplyr::mutate(across(-yday_corr, ~as.integer(stringr::str_extract(.x, "(?<=sum=)\\d+\\.?\\d*")))) |>
dplyr::mutate(yday_corr= dplyr::case_when(
yday_corr == 35~ "Thursday",
yday_corr == 36 ~ "Friday",
yday_corr == 37 ~ "Saturday",
yday_corr == 38 ~ "Sunday",
yday_corr == 21~ "Thursday",
yday_corr == 22 ~ "Friday",
yday_corr == 23 ~ "Saturday",
yday_corr == 24 ~ "Sunday"
)) |>
dplyr::select(yday_corr, GP_Nonhost, Control_Nonhost, GP_Host, Control_Host)|>
dplyr::group_by(yday_corr) |>
summarise(GP_Nonhost= max(GP_Nonhost, na.rm=T),
Control_Nonhost= max(Control_Nonhost, na.rm=T),
GP_Host= max(GP_Host, na.rm=T),
Control_Host= max(Control_Host, na.rm=T)) |>
dplyr::mutate(yday_corr= factor(yday_corr,
levels= c("Thursday", "Friday", "Saturday", "Sunday"))) |>
dplyr::arrange(yday_corr),
#second table
dplyr::mutate(table_mrc_wide, across(-mrc, ~as.integer(stringr::str_extract(.x, "(?<=sum=)\\d+\\.?\\d*"))))
) |>
dplyr::mutate(category= ifelse(is.na(mrc), as.character(yday_corr), as.character(mrc))) |>
dplyr::select(category, everything()) |>
dplyr::select(-any_of(c("mrc", "yday_corr"))) |>
knitr::kable("markdown", caption= "Table 2. High-speed collisions by time windows (two weeks before and during the F1 Grand Prix race weekend) and Non-host/Host MRC-years")| category | GP_Nonhost | Control_Nonhost | GP_Host | Control_Host |
|---|---|---|---|---|
| Thursday | 18 | 16 | 76 | 52 |
| Friday | 25 | 26 | 65 | 74 |
| Saturday | 24 | 23 | 85 | 59 |
| Sunday | 26 | 23 | 66 | 56 |
| laval | 2 | 3 | 39 | 43 |
| longueuil | 5 | 6 | 52 | 44 |
| montreal | 14 | 14 | 201 | 154 |
| quebec | 61 | 58 | NA | NA |
| shrebrooke | 11 | 7 | NA | NA |
df_speed_per_lic_yr <- collisions_weather_corr_rect|>
#4 days , 2 weeks before
mutate(
pre_exposure_post_conv = coalesce(
as.integer(dplyr::between(
date,
race_date_rec - lubridate::days(14+3),
race_date_rec - lubridate::days(4) # miércoles previo al jueves de carrera
)),
0L
),
#2025-12-22: Teniendo en cuenta el desfase de seguimiento
exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - lubridate::days(1),
race_date_rec + lubridate::days(2))), 0L))|>
group_by(mrc_f, year, exposure_post_conv) %>%
reframe(high_speed_tickets_per_lic = (sum(velocidad, na.rm=T)/sum(license_holders_sex_age, na.rm=T))*1e6) |>
ungroup() |>
dplyr::mutate(
mrc = dplyr::case_when(
mrc_f == 66 ~ "montreal",
mrc_f == 65 ~ "laval",
mrc_f == 58 ~ "longueuil",
mrc_f == 23 ~ "quebec",
mrc_f == 43 ~ "shrebrooke",
TRUE ~ as.character(mrc_f)
)
)
df_speed_per_lic_yr |>
filter(exposure_post_conv==0) |>
group_by(mrc) |>
dplyr::summarise(
n = n(),
mean = mean(high_speed_tickets_per_lic, na.rm = TRUE),
median = quantile(high_speed_tickets_per_lic, .5, na.rm = TRUE),
p25 = quantile(high_speed_tickets_per_lic, .25, na.rm = TRUE),
p75 = quantile(high_speed_tickets_per_lic, .75, na.rm = TRUE),
sd = sd(high_speed_tickets_per_lic, na.rm = TRUE),
iqr = IQR(high_speed_tickets_per_lic, na.rm = TRUE),
skew= mean/sd #<2 serious skewness
) |> mutate(print= sprintf("(Mdn= %2.2f; IQR= %2.2f-%2.2f)", median, p25, p75)) |>
knitr::kable("markdown", caption="Yearly high‑speed collision rates by region (per 1,000,000 license holders; May–Sep)")| mrc | n | mean | median | p25 | p75 | sd | iqr | skew | |
|---|---|---|---|---|---|---|---|---|---|
| laval | 23 | 1.7129131 | 1.4906999 | 1.0987943 | 2.4327199 | 0.9135975 | 1.3339256 | 1.874910 | (Mdn= 1.49; IQR= 1.10-2.43) |
| longueuil | 23 | 0.5765156 | 0.5860928 | 0.3753684 | 0.7880717 | 0.2529888 | 0.4127033 | 2.278818 | (Mdn= 0.59; IQR= 0.38-0.79) |
| montreal | 23 | 1.8189576 | 1.4584196 | 1.2037574 | 2.5241584 | 0.8211711 | 1.3204010 | 2.215078 | (Mdn= 1.46; IQR= 1.20-2.52) |
| quebec | 23 | 1.2837297 | 1.3909654 | 0.8413894 | 1.6015473 | 0.5295038 | 0.7601578 | 2.424401 | (Mdn= 1.39; IQR= 0.84-1.60) |
| shrebrooke | 23 | 0.6201147 | 0.5217501 | 0.4135006 | 0.7869631 | 0.3343950 | 0.3734625 | 1.854438 | (Mdn= 0.52; IQR= 0.41-0.79) |
qic_gnm <- function(mod_quasi) {
## 1) φ from quasipoisson
phi <- get_phi(mod_quasi)
## 2) refit Poisson (mismo call) → logLik
mod_pois <- update(mod_quasi, family = poisson)
## 3) qAIC = -2·logLik/φ + 2·k
bbmle::qAIC(mod_pois, dispersion = phi)
}
cat("Treated (w/o Sherbrooke), instantaneous effects\n")
model_gnm_tratados_lead0_sens <- gnm(
sum_velocidad ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(collisions_weather_corr_rect_pre_post_placebo,
grepl("^Pre|^Exp",exposure_window) & treated_sens==TRUE)
)
cat("QIC, 2-day lagged high-speed collisions\n")
qic_gnm(model_gnm_tratados_lead2_sens)
#[1] 524.0127
cat("QIC, 1-day lagged high-speed collisions\n")
qic_gnm(model_gnm_tratados_lead1_sens)
#[1] 530.3881
cat("QIC, high-speed collisions, contemporaneous\n")
qic_gnm(model_gnm_tratados_lead0_sens)
#[1] 530.3753Treated (w/o Sherbrooke), instantaneous effects
QIC, 2-day lagged high-speed collisions
[1] 524.0127
QIC, 1-day lagged high-speed collisions
[1] 530.3881
QIC, high-speed collisions, contemporaneous
[1] 530.3753collisions_weather_corr_rect$non_high_speed<-
collisions_weather_corr_rect$nb_collisions- collisions_weather_corr_rect$velocidad
collisions_weather_corr_rect$prueba<-
collisions_weather_corr_rect$nb_collisions- (collisions_weather_corr_rect$velocidad+ collisions_weather_corr_rect$alcohol)
collisions_weather_corr_rect$nb_blese_any<-
collisions_weather_corr_rect$nb_blese_grave_c+collisions_weather_corr_rect$nb_blese_leger_c
stopifnot(all(collisions_weather_corr_rect$non_high_speed >= 0, na.rm = TRUE))
# 2) Is nb_blese_any consistent with severity components? (if you have them)
# nb_blese_any should often equal deaths + serious + minor (depending on definitions)
summary(collisions_weather_corr_rect$nb_blese_any -
(collisions_weather_corr_rect$nb_mort_c +
collisions_weather_corr_rect$nb_blese_grave_c +
collisions_weather_corr_rect$nb_blese_leger_c))
total_casualty <- with(collisions_weather_corr_rect,
nb_mort_c + nb_blese_grave_c + nb_blese_leger_c)
summary(total_casualty - collisions_weather_corr_rect$velocidad)
sum((total_casualty - collisions_weather_corr_rect$velocidad) < 0, na.rm = TRUE)
cat("“non–high-speed casualty collisions (fatal + injured)")
message("“As a specificity check, we repeated the analysis using
non–high-speed casualty collisions (fatal and non-fatal injuries) as a negative-control outcome")“As a specificity check, we repeated the analysis using
non–high-speed casualty collisions (fatal and non-fatal injuries) as a negative-control outcomecollisions_weather_corr_rect$non_high_speed_lead2_post_conv<- (total_casualty - collisions_weather_corr_rect$velocidad) Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.00000 0.00000 0.00000 -0.04412 0.00000 0.00000
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 2.000 5.000 8.476 10.000 70.000
[1] 0
“non–high-speed casualty collisions (fatal + injured)collisions_weather_corr_rect$non_high_speed_lead2_post_conv <-
collisions_weather_corr_rect|>
group_by(mrc, year)|>
mutate(non_high_speed_lead2_post_conv = dplyr::lead(non_high_speed, 2, order_by = iso_yday))|>
ungroup()|>
pull(non_high_speed_lead2_post_conv)
cat("Last value carried forward approach\n")
collisions_weather_corr_rect|>
group_by(mrc, year)|>
arrange(iso_yday, .by_group = TRUE) |>
# create *_imp copies
mutate(across(all_of(starts_with("non_high_speed_lead2_post_conv")), ~ .x, .names = "{.col}_imp")) |>
# LOCF
tidyr::fill(ends_with("_imp"), .direction = "down") |>
ungroup()|>
(\(df) {
pull(df, non_high_speed_lead2_post_conv_imp) ->> collisions_weather_corr_rect$non_high_speed_lead2_post_conv_imp
})()
nc_collisions_weather_corr_rect_pre_post_lag2 <-
collisions_weather_corr_rect |>
#4 days , 2 weeks before
mutate(
pre_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - lubridate::days((14+3)),
race_date_rec - lubridate::days((14-0)))), 0L),
exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec - lubridate::days(3),
race_date_rec + lubridate::days(0))), 0L),
#4 days , 2 weeks after
neg_exposure_post_conv = coalesce(
as.integer(dplyr::between(date, race_date_rec + lubridate::days(28-3),
race_date_rec + lubridate::days(28+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),
sum_no_velocidad_lead2= sum(non_high_speed_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
nc_rate_lead2_veh = (sum_no_velocidad_lead2 / exp_veh) * 1e6, # not to model
nc_rate_lead2_lic = (sum_no_velocidad_lead2 / 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)
)
# If you already have non_high_speed_lead2:
summary(nc_collisions_weather_corr_rect_pre_post_lag2$sum_no_velocidad_lead2)
sum(nc_collisions_weather_corr_rect_pre_post_lag2$sum_no_velocidad_lead2 < 0, na.rm = TRUE)
# Check complement identity (should be ~0 if definitions match):
# summary((nc_collisions_weather_corr_rect_pre_post_lag2$sum_velocidad_lead2 +
# nc_collisions_weather_corr_rect_pre_post_lag2$sum_no_velocidad_lead2) - nc_collisions_weather_corr_rect_pre_post_lag2$overall_lead2)
tmp <- subset(nc_collisions_weather_corr_rect_pre_post_lag2,
grepl("^Pre|^Exp", exposure_window))
table(tmp$exposure_window) # should look balanced-ish
range(with(tmp, tapply(exposure_window, cluster, function(x) length(unique(x)))))
# Stronger:
any(table(tmp$cluster) != 2)
with(nc_collisions_weather_corr_rect_pre_post_lag2,
table(treated, sub(".*\\.", "", cluster))) # shows MRC codes by treatedLast value carried forward approach
Min. 1st Qu. Median Mean 3rd Qu. Max.
7.00 33.00 46.00 80.28 97.00 320.00
[1] 0
Exposure Pre-exposure
115 115
[1] 2 2
[1] FALSE
treated 23 43 58 65 66
FALSE 69 9 9 9 9
TRUE 0 60 60 60 60cat("Treated (w/ Sherbrooke), effects after 2 days, 4-day summary\n")
nc_model_gnm_tratados_lead2 <- gnm(
sum_no_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2, grepl("^Pre|^Exp",exposure_window) & treated==TRUE)
)
summary(nc_model_gnm_tratados_lead2)
cat("Negative-Control= Controls (w/o Sherbrooke), effects after 2 days, 4-day summary\n")
nc_model_gnm_controles_lead2 <- gnm(
sum_no_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2, grepl("^Pre|^Exp",exposure_window) & treated==FALSE)
)
summary(nc_model_gnm_controles_lead2)
nc_lead2_main <- compare_gnm_effects(nc_model_gnm_controles_lead2,
nc_model_gnm_tratados_lead2,
comparison_name = "Treated vs. Controls (QC only; w/o Sherbrooke)",
coef_index = 1, digits = 3)
# Treated vs. Controls (QC only; w/o Sherbrooke) — coefficient index 1
# z = 1.791; p = 0.073
# log(RRR) = 0.391; se = 0.218
# RRR = 1.478 (95% CI 0.964–2.267)
cat("Negative-Control= Treated (w/o Sherbrooke), effects after 2 days, 4-day summary\n")
nc_model_gnm_tratados_lead2_sens <- gnm(
sum_no_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2 ,
grepl("^Pre|^Exp",exposure_window) & treated_sens==TRUE)
)
summary(nc_model_gnm_tratados_lead2_sens)
cat("Negative-Control= Controls (Quebec+Sherbrooke), effects after 2 days, 4-day summary\n")
nc_model_gnm_controles_lead2_sens <- gnm(
sum_no_velocidad_lead2 ~ relevel(factor(exposure_window), ref="Pre-exposure")+ mean_min_temp_mean_lin + mean_max_temp_mean_lin + median_total_precip_median_lin + median_lag_2_prec_median_imp+# spline cúbico (suaviza)
offset(off_lic),
eliminate = cluster, # <── “conditional” part
family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2 ,
grepl("^Pre|^Exp",exposure_window) & treated_sens==FALSE)
)
summary(nc_model_gnm_controles_lead2_sens)
nc_lead2_sens <- compare_gnm_effects(nc_model_gnm_controles_lead2_sens,
nc_model_gnm_tratados_lead2_sens,
comparison_name = "Treated vs. Controls (QC + Sherbrooke)",
coef_index = 1, digits = 3)
# Treated vs. Controls (QC + Sherbrooke) — coefficient index 1
# z = 1.516; p = 0.129
# log(RRR) = 0.310; se = 0.204
# RRR = 1.363 (95% CI 0.913–2.035)
# Term RR (95% CI) a RRR (95% CI)
# Host vs. Non-host
# Host MRC-years 1.37 (1.12–1.69)
# Non-host MRC-years & host MRCs in ‘09, ‘20, ‘21 1.01 (0.71–1.42)
# 1.36 (0.91–2.04); P= 0.13
# Host MRC-years b 1.36 (1.12–1.65)
# Non-host MRC-years (Quebec City only b) & host MRCs in ‘09, ‘20, ‘21 0.92 (0.63–1.35)
# 1.48 (0.96–2.27); P= 0.07
cbind.data.frame(type=c("controls=QC+Sherbrooke","","controls=QC only",""),
rbind.data.frame(
tidy_exp_gnm_corr(nc_model_gnm_tratados_lead2_sens),
tidy_exp_gnm_corr(nc_model_gnm_controles_lead2_sens),
tidy_exp_gnm_corr(nc_model_gnm_tratados_lead2),
tidy_exp_gnm_corr(nc_model_gnm_controles_lead2)))|>
mutate(RR= sprintf("%.2f (%.2f–%.2f)", `exp(Est.)`, `2.5%`, `97.5%`))|>
mutate(sig= sprintf("%.3f", P))|>
rename("p.value"="P")|>
mutate(term="Exposure")|>
dplyr::select(type, term, RR, sig)|>
knitr::kable("markdown",
caption="Non-High-speed collisions (NC), 4-day summ, 2-day lagged effect", digits=3, na= "-")
cat("Main analysis (sherbrooke+ quebec as controls):\n")
sprintf("RRR = %.2f (95%% CI: %.2f–%.2f), p = %s",
nc_lead2_sens$RRR,
nc_lead2_sens$lo95,
nc_lead2_sens$hi95,
ifelse(nc_lead2_sens$p < 0.001,
"<.001",
sprintf("%.2f", nc_lead2_sens$p)))
cat("Sensitivity analysis (only Quebec as controls):\n")
sprintf("RRR = %.2f (95%% CI: %.2f–%.2f), p = %s",
nc_lead2_main$RRR,
nc_lead2_main$lo95,
nc_lead2_main$hi95,
ifelse(nc_lead2_main$p < 0.001,
"<.001",
sprintf("%.2f", nc_lead2_main$p)))Treated (w/ Sherbrooke), effects after 2 days, 4-day summary
Call:
gnm(formula = sum_no_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2,
grepl("^Pre|^Exp", exposure_window) & treated == TRUE))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.3280749 -0.5732217 0.0001243 0.5655793 2.0242002
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.02539269
mean_min_temp_mean_lin 0.00697976
mean_max_temp_mean_lin 0.00002631
median_total_precip_median_lin 0.01007663
median_lag_2_prec_median_imp 0.00286638
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.02418364
mean_min_temp_mean_lin 0.00962899
mean_max_temp_mean_lin 0.00687030
median_total_precip_median_lin 0.00900557
median_lag_2_prec_median_imp 0.00897756
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.050 0.297
mean_min_temp_mean_lin 0.725 0.471
mean_max_temp_mean_lin 0.004 0.997
median_total_precip_median_lin 1.119 0.267
median_lag_2_prec_median_imp 0.319 0.750
(Dispersion parameter for quasipoisson family taken to be 1.38789)
Residual deviance: 104.85 on 75 degrees of freedom
AIC: NA
Number of iterations: 2
Negative-Control= Controls (w/o Sherbrooke), effects after 2 days, 4-day summary
Call:
gnm(formula = sum_no_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2,
grepl("^Pre|^Exp", exposure_window) & treated == FALSE))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.629e+00 -5.216e-01 7.349e-07 5.147e-01 1.528e+00
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.043356
mean_min_temp_mean_lin 0.003877
mean_max_temp_mean_lin 0.004112
median_total_precip_median_lin 0.011219
median_lag_2_prec_median_imp -0.008775
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.040313
mean_min_temp_mean_lin 0.011894
mean_max_temp_mean_lin 0.009360
median_total_precip_median_lin 0.012522
median_lag_2_prec_median_imp 0.009988
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.075 0.291
mean_min_temp_mean_lin 0.326 0.747
mean_max_temp_mean_lin 0.439 0.664
median_total_precip_median_lin 0.896 0.377
median_lag_2_prec_median_imp -0.879 0.387
(Dispersion parameter for quasipoisson family taken to be 0.9937378)
Residual deviance: 29.901 on 30 degrees of freedom
AIC: NA
Number of iterations: 2
Treated vs. Controls (QC only; w/o Sherbrooke) — coefficient index 1
z = -0.382; p = 0.702
log(RRR) = -0.018; se = 0.047
RRR = 0.982 (95% CI 0.896–1.077)
Negative-Control= Treated (w/o Sherbrooke), effects after 2 days, 4-day summary
Call:
gnm(formula = sum_no_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2,
grepl("^Pre|^Exp", exposure_window) & treated_sens == TRUE))
Deviance Residuals:
Min 1Q Median 3Q Max
-2.27747808 -0.56231958 0.00005884 0.56401268 1.97127623
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.026163
mean_min_temp_mean_lin 0.009887
mean_max_temp_mean_lin -0.002113
median_total_precip_median_lin 0.009771
median_lag_2_prec_median_imp -0.001097
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.026107
mean_min_temp_mean_lin 0.010723
mean_max_temp_mean_lin 0.007673
median_total_precip_median_lin 0.010967
median_lag_2_prec_median_imp 0.009755
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.002 0.321
mean_min_temp_mean_lin 0.922 0.361
mean_max_temp_mean_lin -0.275 0.784
median_total_precip_median_lin 0.891 0.377
median_lag_2_prec_median_imp -0.112 0.911
(Dispersion parameter for quasipoisson family taken to be 1.497436)
Residual deviance: 82.944 on 55 degrees of freedom
AIC: NA
Number of iterations: 2
Negative-Control= Controls (Quebec+Sherbrooke), effects after 2 days, 4-day summary
Call:
gnm(formula = sum_no_velocidad_lead2 ~ relevel(factor(exposure_window),
ref = "Pre-exposure") + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp +
offset(off_lic), eliminate = cluster, family = quasipoisson(link = "log"),
data = subset(nc_collisions_weather_corr_rect_pre_post_lag2,
grepl("^Pre|^Exp", exposure_window) & treated_sens ==
FALSE))
Deviance Residuals:
Min 1Q Median 3Q Max
-1.7413298 -0.5233244 0.0002809 0.4942871 1.6079154
Coefficients of interest:
Estimate
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.0395171
mean_min_temp_mean_lin -0.0006133
mean_max_temp_mean_lin 0.0067849
median_total_precip_median_lin 0.0121251
median_lag_2_prec_median_imp -0.0027687
Std. Error
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 0.0364237
mean_min_temp_mean_lin 0.0107272
mean_max_temp_mean_lin 0.0082671
median_total_precip_median_lin 0.0101113
median_lag_2_prec_median_imp 0.0096483
t value Pr(>|t|)
relevel(factor(exposure_window), ref = "Pre-exposure")Exposure 1.085 0.283
mean_min_temp_mean_lin -0.057 0.955
mean_max_temp_mean_lin 0.821 0.416
median_total_precip_median_lin 1.199 0.236
median_lag_2_prec_median_imp -0.287 0.775
(Dispersion parameter for quasipoisson family taken to be 1.048714)
Residual deviance: 52.694 on 50 degrees of freedom
AIC: NA
Number of iterations: 2
Treated vs. Controls (QC + Sherbrooke) — coefficient index 1
z = -0.298; p = 0.766
log(RRR) = -0.013; se = 0.045
RRR = 0.987 (95% CI 0.904–1.077)| type | term | RR | sig |
|---|---|---|---|
| controls=QC+Sherbrooke | Exposure | 1.03 (0.98–1.08) | 0.316 |
| Exposure | 1.04 (0.97–1.12) | 0.278 | |
| controls=QC only | Exposure | 1.03 (0.98–1.08) | 0.294 |
| Exposure | 1.04 (0.96–1.13) | 0.282 |
Main analysis (sherbrooke+ quebec as controls):
[1] "RRR = 0.99 (95% CI: 0.90–1.08), p = 0.77"
Sensitivity analysis (only Quebec as controls):
[1] "RRR = 0.98 (95% CI: 0.90–1.08), p = 0.70"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/step3.RData")R library: H:/My Drive/PERSONAL ANDRES/UCH_salud_publica/pasantia/f1/f1/renv/library/windows/R-4.4/x86_64-w64-mingw32Date: 2026-01-02 13:12:41.813599Editor 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/Rtmp2BwyLy/file11aa8485c5c31 -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 |
| backports | 1.5.0 | RSPM |
| bayesplot | 1.12.0 | RSPM |
| bbmle | 1.0.25.1 | 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) |
| boot | 1.3-30 | CRAN (R 4.4.1) |
| bridgesampling | 1.1-2 | RSPM |
| brms | 2.22.0 | RSPM |
| Brobdingnag | 1.2-9 | RSPM |
| broom | 1.0.8 | RSPM |
| cachem | 1.1.0 | CRAN (R 4.4.3) |
| car | 3.1-3 | RSPM |
| carData | 3.0-5 | RSPM |
| CBPS | 0.23 | RSPM |
| checkmate | 2.3.2 | RSPM |
| cli | 3.6.5 | RSPM |
| cmprsk | 2.2-12 | 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 |
| cowplot | 1.1.3 | RSPM |
| curl | 6.2.3 | CRAN (R 4.4.1) |
| CVXR | 1.0-15 | RSPM |
| data.table | 1.17.4 | RSPM |
| devtools | 2.4.5 | RSPM |
| DHARMa | 0.4.7 | RSPM |
| digest | 0.6.37 | RSPM |
| distributional | 0.5.0 | RSPM |
| dplyr | 1.1.4 | RSPM |
| dreamerr | 1.5.0 | RSPM |
| ellipsis | 0.3.2 | RSPM |
| emmeans | 1.11.1 | RSPM |
| Epi | 2.60 | RSPM |
| estimability | 1.5.1 | RSPM |
| etm | 1.1.2 | RSPM |
| evaluate | 1.0.3 | RSPM |
| farver | 2.1.2 | RSPM |
| fastmap | 1.2.0 | CRAN (R 4.4.3) |
| fixest | 0.12.1 | RSPM |
| foreach | 1.5.2 | RSPM |
| Formula | 1.2-5 | 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 |
| ggplot2 | 3.5.2 | RSPM |
| glmmTMB | 1.1.11 | RSPM |
| glmnet | 4.1-9 | 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 |
| htmltools | 0.5.8.1 | RSPM |
| htmlwidgets | 1.6.4 | RSPM |
| httpuv | 1.6.16 | RSPM |
| httr2 | 1.1.2 | RSPM |
| inline | 0.3.21 | RSPM |
| iterators | 1.0.14 | RSPM |
| jsonlite | 2.0.0 | CRAN (R 4.4.3) |
| kableExtra | 1.4.0 | RSPM |
| knitr | 1.50 | 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) |
| 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 |
| 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) |
| numDeriv | 2016.8-1.1 | 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 |
| processx | 3.8.6 | RSPM |
| profvis | 0.4.0 | RSPM |
| promises | 1.3.2 | RSPM |
| ps | 1.9.1 | RSPM |
| purrr | 1.0.4 | RSPM |
| quarto | 1.4.4 | RSPM |
| QuickJSR | 1.7.0 | 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 |
| reformulas | 0.4.1 | RSPM |
| relimp | 1.0-5 | RSPM |
| remotes | 2.5.0 | RSPM |
| renv | 1.1.2 | CRAN (R 4.4.1) |
| rio | 1.2.3 | RSPM |
| rlang | 1.1.6 | RSPM |
| rmarkdown | 2.29 | RSPM |
| Rmpfr | 1.1-0 | RSPM |
| rstan | 2.32.7 | RSPM |
| rstantools | 2.4.0 | RSPM |
| rstudioapi | 0.17.1 | RSPM |
| sandwich | 3.1-1 | RSPM |
| scales | 1.4.0 | RSPM |
| sessioninfo | 1.2.3 | RSPM |
| shape | 1.4.6.1 | RSPM |
| shiny | 1.10.0 | RSPM |
| sp | 2.2-0 | RSPM |
| SparseM | 1.84-2 | RSPM |
| StanHeaders | 2.32.10 | RSPM |
| stringi | 1.8.7 | RSPM |
| stringmagic | 1.2.0 | RSPM |
| stringr | 1.5.1 | RSPM |
| survival | 3.6-4 | CRAN (R 4.4.1) |
| svglite | 2.2.1 | RSPM |
| systemfonts | 1.2.3 | RSPM |
| tensorA | 0.36.2.1 | RSPM |
| textshaping | 1.0.1 | RSPM |
| tibble | 3.2.1 | RSPM |
| tidyr | 1.3.1 | RSPM |
| tidyselect | 1.2.1 | RSPM |
| timechange | 0.3.0 | RSPM |
| timeDate | 4041.110 | RSPM |
| TMB | 1.9.17 | 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 |
| withr | 3.0.2 | RSPM |
| xfun | 0.52 | RSPM |
| xml2 | 1.3.8 | CRAN (R 4.4.3) |
| xtable | 1.8-4 | RSPM |
| yaml | 2.3.10 | 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