Step 3. Check analytical strategies (part 2)
# remove objects and memory
rm(list=ls());gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 840104 44.9 1664252 88.9 1127548 60.3
Vcells 1707842 13.1 8388608 64.0 3505321 26.8#remove images
while(!dev.cur())dev.off()
cat("\014")load(paste0(getwd(),"/_data/step2.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(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)
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 number
#nixtlar::nixtla_client_setup(api_key = readLines(paste0(gsub("f1/","", getwd() ),"/key.txt"))[[1]])
nixtlar::nixtla_client_setup(api_key = readLines(paste0(gsub("f1/","", getwd() ),"/key_alt.txt"))[[1]])Warning in readLines(paste0(gsub("f1/", "", getwd()), "/key_alt.txt")):
incomplete final line found on 'H:/My Drive/PERSONAL
ANDRES/UCH_salud_publica/pasantia/f1/key_alt.txt'API key has been set for the current session.cat("Treated (w/o Sherbrooke), effects after 2 days\n")Treated (w/o Sherbrooke), effects after 2 daysb_cmp <- brm(
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 + (1 | cluster) + offset(off_lic),#median_lag_2_prec_median_imp+
family = negbinomial(),#poisson(),
data = subset(collisions_weather_corr_rect_pre_post, grepl("^Pre|^Exp",exposure_window) & treated_sens==TRUE),
prior = c(
prior(normal(0, 5), class = "b"),#,prior(exponential(1), class = "nu") # ν > 0
prior(gamma(0.01, 0.01), class = "shape") # θ controla la sobre-disp.
),
chains = 4, cores = 16, iter = 4000*10, warmup = 4000*2,
control = list(adapt_delta = .99)
)Compiling Stan program...Start samplinga_cmp <- brm(
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 +(1 | cluster) + offset(off_lic),# median_lag_2_prec_median_imp+
family = negbinomial(),#poisson(),
data = subset(collisions_weather_corr_rect_pre_post, grepl("^Pre|^Exp",exposure_window) & treated_sens==FALSE),
prior = c(
prior(normal(0, 5), class = "b"),#,#,prior(exponential(1), class = "nu") # ν > 0
prior(gamma(0.01, 0.01), class = "shape") # θ controla la sobre-disp.
),
chains = 4, cores = 16, iter = 4000*10, warmup = 4000*2,
control = list(adapt_delta = .99)
)Compiling Stan program...
Start sampling# # 2.1 Histograma/overlay de densidades
# pp_check(b_cmp, type = "dens_overlay", ndraws = 200)
# pp_check(a_cmp, type = "dens_overlay", ndraws = 200)
#
# # 2.2 Rootogram (bueno para conteos)
# pp_check(b_cmp, type = "rootogram")
# pp_check(a_cmp, type = "rootogram")
#
# # 2.3 Distribución del número de ceros
# pp_check(b_cmp, type = "bars", statistic = "zero")
#
# # 2.4 Estadísticos resumidos
# pp_check(
# b_cmp,
# type = "stat",
# statistic = c("mean", "sd", "var", "max") # elige los que te interesen
# )#relevelfactorexposure_windowrefEQPreMexposureExposure
#
# y_rep <- posterior_predict(b_cmp, draws = 500) # 500 réplicas
# summ_mean <- apply(y_rep, 1, function(x) c(mean = mean(x)))
# summ_var <- apply(y_rep, 1, function(x) c(var = var(x)))
# df_predict_mean<- sapply(summ_mean, FUN = identity)
# df_predict_var<- sapply(summ_var, FUN = identity)
# tibble(mean= as.numeric(df_predict_mean), var= as.numeric(df_predict_var))|>
# ggplot(aes(x = mean, y = var)) +
# geom_point(alpha = .3) +
# geom_abline(slope = 1, intercept = 0, colour = "red") +
# labs(title = "Media vs Varianza en réplicas -> línea roja = Poisson ideal")
# #MUCHA SOBREDISPERSIÓN
#
# Extrae draws del coeficiente, p. ej. el de `exposure_windowExp`
betaA <- posterior_samples(a_cmp, "b_relevelfactorexposure_windowrefEQPreMexposureExposure")$b_relevelfactorexposure_windowrefEQPreMexposureExposureWarning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.betaB <- posterior_samples(b_cmp, "b_relevelfactorexposure_windowrefEQPreMexposureExposure")$b_relevelfactorexposure_windowrefEQPreMexposureExposureWarning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
recommended alternatives.cat("Cat B= Host-MRC-years\n")Cat B= Host-MRC-yearsmean(exp(betaB)); quantile(exp(betaB), c(.025, .975))[1] 1.354806 2.5% 97.5%
1.075355 1.691926 cat("Cat A= Non-Host-MRC-years\n")Cat A= Non-Host-MRC-yearsmean(exp(betaA)); quantile(exp(betaA), c(.025, .975))[1] 1.128858 2.5% 97.5%
0.7456823 1.6510213 cat("Delta\n")Deltadelta <- betaA - betaB # diferencia A − B
# Resumen con intervalo de credibilidad
quantile(delta, c(.025, .5, .975)) 2.5% 50% 97.5%
-0.6492694 -0.1970737 0.2623611 #> 2.5% 50% 97.5%
#> -0.12 | 0.08 | 0.28 (ejemplo)
#>
#Si el 95 % CrI de delta incluye 0 → el efecto puede ser igual en ambas poblaciones.
library(posterior)This is posterior version 1.6.1
Adjuntando el paquete: 'posterior'The following object is masked from 'package:bayesplot':
rhatThe following objects are masked from 'package:stats':
mad, sd, varThe following objects are masked from 'package:base':
%in%, matchdelta_draws <- as_draws_matrix(matrix(delta, ncol = 1))
colnames(delta_draws) <- "delta" # nombre de la columna → etiqueta en la leyenda
# 2. Grafica el área de credibilidad del 95 %
mcmc_areas(delta_draws, pars = "delta", prob = 0.95) +
geom_vline(xintercept = 0, linetype = "dashed") +
labs(title = "Δβ (A − B)",
x = "Coefficient difference",
y = NULL)
summ <- summarise_draws(
delta_draws,
prob_gt0 = ~ mean(.x > 0),
prob_lt0 = ~ mean(.x < 0),
prob_rope = ~ mean(abs(.x) < 0.05) # ROPE de ±0.05 “en línea”
)
summ# A tibble: 1 × 4
variable prob_gt0 prob_lt0 prob_rope
<chr> <dbl> <dbl> <dbl>
1 delta 0.199 0.801 0.119#ROPE: El 11.9 % de las simulaciones cae en la zona de irrelevancia práctica que fijaste.
# Usando reglas heurísticas típicas:
# • < 0.75 → débil, 0.75–0.90 → moderada, 0.90–0.97 → fuerte, > 0.97 → muy fuerte.
# Con 0.82 estás en la franja moderada: hay más apoyo por un efecto positivo que negativo, pero la incertidumbre todavía es apreciable.df <- subset(
collisions_weather_corr_rect_synth,
yday_corr <= (max_time - 7) # mismo recorte temporal que en countSynth
)
df$treat_unit <- df$tr_contrl_corr
df$post_treat <- ifelse(df$yday_corr >= 35, 1, 0)
df$D <- df$tr_contrl_corr * df$post_treat # interacción exposición x tiempo
df$log_rate <- log(df$sum_velocidad_lead2) - log(df$exp_lic) # log(count) – log(exposure)
#table(collisions_weather_corr_rect_synth$id_lab,collisions_weather_corr_rect_synth$tr_contrl_corr)
library(fect)
set.seed(2125)
road.main.gsynth <- fect::fect(
formula = rate_lead2_lic ~ D + mean_min_temp_mean_lin + mean_max_temp_mean_lin +#rate_lead2_lic
median_total_precip_median_lin + median_lag_2_prec_median_imp,
data = df,
index = c("id", "yday_corr"), # unit.variable y time.variable
force = "two-way", # fijamos efectos fijos unitarios y temporales
CV = TRUE, # selección automática de r con cross‑validation
r = 1:4, # ≈ K = 4 en tu modelo original
se = TRUE, # errores estándar por bootstrap
nboots = 1000, # ajusta según tu tolerancia de tiempo/precisión
parallel= F, # usa varios hilos si tu máquina lo permite
method= "mc"
)Cross-validating ...Criterion: Mean Squared Prediction ErrorMatrix completion method...
lambda.norm = 1.00000; MSPE = 1.44867; MSPTATT = 0.02709; MSE = 0.59267*
lambda.norm = 0.42170; MSPE = 1.61655; MSPTATT = 0.01223; MSE = 0.36869
lambda.norm = 0.17783; MSPE = 1.80265; MSPTATT = 0.00358; MSE = 0.11313
lambda.norm = 0.07499; MSPE = 1.84662; MSPTATT = 0.00072; MSE = 0.02422
lambda.norm* = 1Bootstrapping for uncertainties ... .
.
.
.
.
.
.
.
.
.1000 runsCannot use full pre-treatment periods in the F test. The first period is removed.# ATT:
# ATT S.E. CI.lower CI.upper p.value
# Tr obs equally weighted 0.223 0.1645 -0.09942 0.5455 0.1752
# Tr units equally weighted 0.223 0.1645 -0.09942 0.5455 0.1752
road.main.gsynth
plot(road.main.gsynth, main="-3 pre to +0, Quebec+Sherbrooke")
df2 <- subset(
collisions_weather_corr_rect_quebec_synth,
yday_corr <= (max_time - 7) # mismo recorte temporal que en countSynth
)
df2$treat_unit <- df2$tr_contrl_corr
df2$post_treat <- ifelse(df2$yday_corr >= 35, 1, 0)
df2$D <- df2$tr_contrl_corr * df2$post_treat # interacción exposición x tiempo
road.main.gsynth2 <- fect::fect(
formula = rate_lead2_lic ~ D + mean_min_temp_mean_lin + mean_max_temp_mean_lin +
median_total_precip_median_lin + median_lag_2_prec_median_imp,
data = df2,
index = c("id", "yday_corr"), # unit.variable y time.variable
force = "two-way", # fijamos efectos fijos unitarios y temporales
CV = TRUE, # selección automática de r con cross‑validation
r = 1:4, # ≈ K = 4 en tu modelo original
se = TRUE, # errores estándar por bootstrap
nboots = 1000, # ajusta según tu tolerancia de tiempo/precisión
parallel= F, # usa varios hilos si tu máquina lo permite
method= "mc"
)Cross-validating ...Criterion: Mean Squared Prediction ErrorMatrix completion method...
lambda.norm = 1.00000; MSPE = 2.00240; MSPTATT = 0.02521; MSE = 0.51710*
lambda.norm = 0.42170; MSPE = 2.32997; MSPTATT = 0.01583; MSE = 0.36842
lambda.norm = 0.17783; MSPE = 2.60047; MSPTATT = 0.00497; MSE = 0.13719
lambda.norm = 0.07499; MSPE = 2.65169; MSPTATT = 0.00100; MSE = 0.03253
lambda.norm* = 1Bootstrapping for uncertainties ... .
.
.
.
.
.
.
.
.
.1000 runsCannot use full pre-treatment periods in the F test. The first period is removed.cat("Quebec only as controls\n")
road.main.gsynth2
# ATT:
# ATT S.E. CI.lower CI.upper p.value
# Tr obs equally weighted 0.2647 0.1803 -0.08866 0.6181 0.142
# Tr units equally weighted 0.2647 0.1803 -0.08866 0.6181 0.142
plot(road.main.gsynth2, main="-3 pre to +0, Quebec only")
Call:
fect.formula(formula = rate_lead2_lic ~ D + mean_min_temp_mean_lin +
mean_max_temp_mean_lin + median_total_precip_median_lin +
median_lag_2_prec_median_imp, data = df, index = c("id",
"yday_corr"), force = "two-way", r = 1:4, CV = TRUE, method = "mc",
se = TRUE, nboots = 1000, parallel = F)
ATT:
ATT S.E. CI.lower CI.upper p.value
Tr obs equally weighted 0.2667 0.1757 -0.07765 0.611 0.129
Tr units equally weighted 0.2667 0.1757 -0.07765 0.611 0.129
Covariates:
Coef S.E. CI.lower CI.upper p.value
mean_min_temp_mean_lin 0.0215623 0.009682 0.002586 0.040538 0.02594
mean_max_temp_mean_lin -0.0022321 0.008198 -0.018300 0.013836 0.78542
median_total_precip_median_lin -0.0050004 0.005643 -0.016060 0.006059 0.37552
median_lag_2_prec_median_imp -0.0004109 0.003537 -0.007343 0.006521 0.90751
Quebec only as controls
Call:
fect.formula(formula = rate_lead2_lic ~ D + mean_min_temp_mean_lin +
mean_max_temp_mean_lin + median_total_precip_median_lin +
median_lag_2_prec_median_imp, data = df2, index = c("id",
"yday_corr"), force = "two-way", r = 1:4, CV = TRUE, method = "mc",
se = TRUE, nboots = 1000, parallel = F)
ATT:
ATT S.E. CI.lower CI.upper p.value
Tr obs equally weighted 0.2853 0.1797 -0.06692 0.6375 0.1124
Tr units equally weighted 0.2853 0.1797 -0.06692 0.6375 0.1124
Covariates:
Coef S.E. CI.lower CI.upper p.value
mean_min_temp_mean_lin 0.0207474 0.010834 -0.0004867 0.04198 0.05549
mean_max_temp_mean_lin -0.0032537 0.008305 -0.0195315 0.01302 0.69523
median_total_precip_median_lin 0.0033263 0.007128 -0.0106444 0.01730 0.64075
median_lag_2_prec_median_imp -0.0001074 0.004392 -0.0087147 0.00850 0.98048We defined databases for (t_int_man) = 38 days (df35) as pre-intervention period and a more restrictive (t_int_man-3) = 35 (df33).
#https://cran.r-project.org/web/packages/nixtlar/vignettes/prediction-intervals.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/exogenous-variables.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/cross-validation.html
#https://cran.r-project.org/web/packages/nixtlar/vignettes/anomaly-detection.html
collisions_weather_corr_rect_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_synth_tr|>
mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>
mutate(Y_lag10 = dplyr::lag(rate_veh, 10, order_by = yday_corr))|>
ungroup()|>
pull(Y_lag10)
collisions_weather_corr_rect_synth_tr$Y_lag10 <-
collisions_weather_corr_rect_synth_tr|>
mutate(
third_value = nth(na.omit(Y_lag10), 1),
Y_lag10 = ifelse(is.na(Y_lag10), third_value, Y_lag10)
)|>
dplyr::select(Y_lag10)|>
ungroup()|>
pull(Y_lag10)
df33 <- collisions_weather_corr_rect_synth_tr|>
arrange(year.x, yday_corr)|> # 1,2,3… days relative to race
filter(yday_corr<(t_int_man-3)) |>
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
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,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
rate_lead1_lic= rate_lead1_lic,
rate_lead2_lic= rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic,
Y_lag10= Y_lag10
)
df33c <- collisions_weather_corr_rect_synth_cntr|>
arrange(year.x, yday_corr)|> # 1,2,3… days relative to race
filter(yday_corr<(t_int_man-3)) |>
transmute(
unique_id = "control", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
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,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
rate_lead1_lic= rate_lead1_lic,
rate_lead2_lic= rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
df33_post <- collisions_weather_corr_rect_synth_tr|>
arrange(year.x, yday_corr)|> # 1,2,3… days relative to race
filter(yday_corr>=(t_int_man-3), yday_corr<=max_time) |> # take the first 30
transmute(
unique_id = "treated", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
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,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
rate_lead1_lic= rate_lead1_lic,
rate_lead2_lic= rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic,
Y_lag10= Y_lag10
)
df33c_post <- collisions_weather_corr_rect_synth_cntr|>
arrange(year.x, yday_corr)|> # 1,2,3… days relative to race
filter(yday_corr>=(t_int_man-3), yday_corr<=max_time) |> # take the first 30
transmute(
unique_id = "control", # exactly this name
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
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,
mean_min_temp_mean_lin= mean_min_temp_mean_lin,
mean_max_temp_mean_lin= mean_max_temp_mean_lin,
rate_lead1_lic= rate_lead1_lic,
rate_lead2_lic= rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
fc33 <- nixtla_client_forecast(
rbind(mutate(df33[1:1,], ds=ds-1), df33),
h = 7,
level = c(95),
freq = "D", # <— this fixes the infer_frequency error
finetune_steps = 10, #Number of steps used to finetune 'TimeGPT' in the new data.
finetune_loss ="rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, rate_lead1_lic, rate_lead2_lic, sum_velocidad, sum_velocidad_lead1, sum_velocidad_lead2, exp_lic, Y_lag10p33 <- nixtla_client_plot(rbind(mutate(df33[1:1,], ds=ds-1), df33, df33_post), fc33,
max_insample_length = 100)
p33$layers[[1]]$aes_params$fill <- "#eb8f8f" # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5
# 1) Banda de CI (layer 1)
p33$layers[[1]]$geom_params$fill <- "#FFC10740" # amarillo semitransparente
p33$layers[[1]]$geom_params$alpha <- 0.4
# 2) Línea histórica y pronóstico (layer 2)
# p$layers[[2]]$aes_params$colour <- "steelblue"
p33$layers[[2]]$aes_params$linewidth <- 1.0
p33$layers[[1]]$aes_params$fill <- "#eb8f8f" # verde semitransparente
p33$layers[[1]]$aes_params$alpha <- 0.5
t0_33 <- min(p33$data$ds)
t0_33_posix <- min(p33$data$ds)
xint35 <- t0_33_posix + days((t_int_man-3))
t0_33 <- as.Date(min(p33$data$ds))
p33+
scale_color_manual(
values = c(
# y = "#FFD600", # amarillo intenso
TimeGPT = "red" # violeta (DarkOrchid / BlueViolet)
),
labels= c( TimeGPT = "Prediction (95% IC in gray area)" )
) +
scale_x_datetime(
name = "Days of follow-up (33 days before race,\n7 days after the race)",
limits = c(min(p33$data$ds), max(p33$data$ds)),
breaks = seq(min(p33$data$ds), max(p33$data$ds), by = "7 days"),
labels = function(x) {
# convierto POSIXct a Date y resto t0
as.integer(as.Date(x) - t0_33 + 1)
}
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(), # quita el label de faceta
strip.background = element_blank(),
legend.position= "bottom"
)+
geom_vline(
xintercept = xint35-days(2),
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y="Highspeed collisions per 1M license holders")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
## ── Conjuntos 33 días pre-intervención y 7 post ───────────────────────────────
df33_qc <- collisions_weather_corr_rect_quebec_synth_tr|>
arrange(year.x, yday_corr)|>
filter(yday_corr < (t_int_man-3))|>
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
rate_lead1_lic,
rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
df33c_qc <- collisions_weather_corr_rect_quebec_synth_cntr|>
arrange(year.x, yday_corr)|>
filter(yday_corr < (t_int_man-3))|>
transmute(
unique_id = "control",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
rate_lead1_lic,
rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
df33_post_qc <- collisions_weather_corr_rect_quebec_synth_tr|>
arrange(year.x, yday_corr)|>
filter(yday_corr >= (t_int_man-3), yday_corr <= max_time)|>
transmute(
unique_id = "treated",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
rate_lead1_lic,
rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
df33c_post_qc <- collisions_weather_corr_rect_quebec_synth_cntr|>
arrange(year.x, yday_corr)|>
filter(yday_corr >= (t_int_man-3), yday_corr <= max_time)|>
transmute(
unique_id = "control",
ds = as.POSIXct(as.Date("2019-05-19") + (yday_corr - 1)),
y = rate_veh,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin,
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
rate_lead1_lic,
rate_lead2_lic,
sum_velocidad= sum_velocidad,
sum_velocidad_lead1= sum_velocidad_lead1,
sum_velocidad_lead2= sum_velocidad_lead2,
exp_lic= exp_lic
)
## ── TimeGPT ───────────────────────────────────────────────────────────────────
fc33_qc <- nixtla_client_forecast(
rbind(mutate(df33_qc[1:1,], ds = ds - 1), df33_qc),
h = 7,
level = c(95),
freq = "D",
finetune_steps = 10,
finetune_loss = "rmse"
)Using historical exogenous features: mean_median_lag_2_prec_median_imp, mean_median_total_precip_median_lin, mean_min_temp_mean_lin, mean_max_temp_mean_lin, rate_lead1_lic, rate_lead2_lic, sum_velocidad, sum_velocidad_lead1, sum_velocidad_lead2, exp_licp33_qc <- nixtla_client_plot(
rbind(mutate(df33_qc[1:1,], ds = ds - 1),
df33_qc,
df33_post_qc),
fc33_qc,
max_insample_length = 100
)
## ── Estética idéntica a tu original ───────────────────────────────────────────
p33_qc$layers[[1]]$geom_params$fill <- "#FFC10740"
p33_qc$layers[[1]]$geom_params$alpha <- 0.4
p33_qc$layers[[2]]$aes_params$linewidth <- 1.0
t0_33_qc <- as.Date(min(p33_qc$data$ds))
xint33_qc <- min(p33_qc$data$ds) + lubridate::days((t_int_man-3))
p33_qc +
scale_color_manual(
values = c(TimeGPT = "red"),
labels = c(TimeGPT = "Prediction (95% IC in gray area)")
) +
scale_x_datetime(
name = "Days of follow-up (33 before, 7 after)",
limits = c(min(p33_qc$data$ds), max(p33_qc$data$ds)),
breaks = seq(min(p33_qc$data$ds), max(p33_qc$data$ds), by = "7 days"),
labels = function(x) as.integer(as.Date(x) - t0_33_qc + 1)
) +
theme_classic(base_size = 13) +
theme(
strip.text = element_blank(),
strip.background = element_blank(),
legend.position = "bottom"
) +
geom_vline(
xintercept = xint33_qc,
linetype = "dashed",
color = "black",
linewidth = 1
) +
labs(y = "High-speed collisions per 1 M license holders")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
We used the CausalImpact package, useful for time series data, inferring a counterfactual post race period based on pre-intervention components1. Given that the offset was constant, we used the sum of the high-speed collisions instead of rates.
library(bsts)
# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates33 <- cbind(
# y = as.numeric(df35$y),
mean_min_temp_mean_lin = df33$mean_min_temp_mean_lin,
mean_max_temp_mean_lin = df33$mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp = df33$mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= df33$mean_median_total_precip_median_lin#,
#Y_lag10 = df33$Y_lag10
)
# Ensure names explicitly
colnames(covariates33) <- c(#"y",
"mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,#"Y_lag10")
)
covariates33<- as.data.frame(covariates33)
covariates_future33 <- as.matrix(df33_post[1:7, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,
#"Y_lag10"
)])
covariates_full33 <- rbind(covariates33, covariates_future33)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33$sum_velocidad_lead2), lambda_opt)
y33b <- c(as.numeric(df33$sum_velocidad_lead2), rep(NA, 7)) # respuesta escalada por 100 millones
y33b <- log(y33b)
y33b <- c(y_bc, rep(NA, 7)) # respuesta escalada por 100 millones
ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
#ss2b <- AddAutoAr(ss2b, y33b, lags = 1)
ss2b <- AddDynamicRegression(ss2b, covariates_full33)
set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
niter = 4e4, family = "gaussian")
#plot(model, "components", burn = 2000)
impact3d_ratio_resp2 <- CausalImpact(bsts.model = model33b,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = as.numeric(BoxCox(as.numeric(df33_post$sum_velocidad_lead2[1:7]), lambda_opt)))
#log(df33_post$sum_velocidad_lead2[1:7]))
#plot(impact3d, "original")
summary(impact3d_ratio_resp2)
plot(impact3d_ratio_resp2)
err_ci_00 <- impact3d_ratio_resp2$series$response - impact3d_ratio_resp2$series$point.pred
rmse_ci_00 <- sqrt(mean(err_ci_00[1:34]^2))
mean_pre_00 <- mean(df33$sum_velocidad_lead2[1:34])
sd_pre_00 <- sd(df33$sum_velocidad_lead2[1:34])
nrmse_00 <- rmse_ci_00 / mean_pre_00 # RMSE relativo;
sd_ratio_00 <- rmse_ci_00 / sd_pre_00 # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_00))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_00, mean_pre_00))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_00, sd_pre_00))
# Benchmark naïve (random walk de un paso):
naive_rmse_00 <- sqrt(mean((df33$sum_velocidad_lead2[2:34] - df33$sum_velocidad_lead2[1:33])^2))
improvement_00 <- 1 - rmse_ci_00 / naive_rmse_00 # >0 = mejora sobre el naïve=-=-=-=-= Iteration 0 Tue Aug 19 16:08:26 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:08:27 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:08:29 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:08:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:08:31 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:08:33 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:08:34 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:08:35 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:08:37 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:08:38 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 25 172
Prediction (s.d.) 21 (2.5) 150 (17.6)
95% CI [17, 27] [118, 188]
Absolute effect (s.d.) 3.2 (2.5) 22.6 (17.6)
95% CI [-2.2, 7.8] [-15.3, 54.4]
Relative effect (s.d.) 17% (18%) 17% (18%)
95% CI [-8.2%, 46%] [-8.2%, 46%]
Posterior tail-area probability p: 0.09249
Posterior prob. of a causal effect: 91%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 1.56 (average absolute forecast error)
NRMSE: 0.03 (RMSE relative to pre-period mean, 54.6)
SD-ratio: 0.16 (RMSE vs pre-period variability, SD=9.6)# Add external regressors explicitly (e.g., temperatures, precipitation)
covariates33c <- cbind(
# y = as.numeric(df35$y),
mean_min_temp_mean_lin = df33c$mean_min_temp_mean_lin,
mean_max_temp_mean_lin = df33c$mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp = df33c$mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin= df33c$mean_median_total_precip_median_lin#,
#Y_lag10 = df33$Y_lag10
)
# Ensure names explicitly
colnames(covariates33c) <- c(#"y",
"mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,#"Y_lag10")
)
covariates33c<- as.data.frame(covariates33c)
covariates_future33c <- as.matrix(df33c_post[1:7, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,
#"Y_lag10"
)])
covariates_full33c <- rbind(covariates33c, covariates_future33c)
# y33cb <- c(as.numeric(df33c$sum_velocidad_lead2), rep(NA, 7)) # respuesta escalada por 100 millones
# y33cb <- log(y33cb)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33c$sum_velocidad_lead2), lambda_opt)
y33cb <- c(y_bc, rep(NA, 7)) # respuesta escalada por 100 millones
ss2cb <- list()
ss2cb <- AddLocalLinearTrend(ss2cb, y33cb)
ss2cb <- AddSeasonal(ss2cb, y33cb, nseasons = 7, season.duration = 1)
#ss2b <- AddAutoAr(ss2b, y33cb, lags = 1)
ss2cb <- AddDynamicRegression(ss2cb, covariates_full33c)
set.seed(2125)
model33cb <- bsts(y33cb, state.specification = ss2cb, seed=2125,
niter = 4e4, family = "gaussian")
#plot(model, "components", burn = 2000)
impact3d_ratio_resp2_c <- CausalImpact(bsts.model = model33cb,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = as.numeric(BoxCox(as.numeric(df33c_post$sum_velocidad_lead2[1:7]), lambda_opt)))
#log(df33c_post$sum_velocidad_lead2[1:7]))
#plot(impact3d, "original")
summary(impact3d_ratio_resp2_c)
plot(impact3d_ratio_resp2_c)
err_ci_0c <- impact3d_ratio_resp2_c$series$response - impact3d_ratio_resp2_c$series$point.pred
rmse_ci_0c <- sqrt(mean(err_ci_0c[1:34]^2))
mean_pre_0c <- mean(df33c$sum_velocidad_lead2[1:34])
sd_pre_0c <- sd(df33c$sum_velocidad_lead2[1:34])
nrmse_0c <- rmse_ci_0c / mean_pre_0c # RMSE relativo;
sd_ratio_0c <- rmse_ci_0c / sd_pre_0c # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_0c))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_0c, mean_pre_0c))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_0c, sd_pre_0c))
# Benchmark naïve (random walk de un paso):
naive_rmse_0c <- sqrt(mean((df33c$sum_velocidad_lead2[2:34] - df33c$sum_velocidad_lead2[1:33])^2))
improvement_0c <- 1 - rmse_ci_0c / naive_rmse_0c # >0 = mejora sobre el naïve=-=-=-=-= Iteration 0 Tue Aug 19 16:08:42 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:08:44 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:08:45 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:08:46 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:08:48 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:08:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:08:50 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:08:52 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:08:53 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:08:54 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 0.74 5.20
Prediction (s.d.) 0.75 (0.13) 5.22 (0.90)
95% CI [0.49, 1] [3.45, 7]
Absolute effect (s.d.) -0.0021 (0.13) -0.0144 (0.90)
95% CI [-0.26, 0.25] [-1.81, 1.75]
Relative effect (s.d.) 3.1% (20%) 3.1% (20%)
95% CI [-26%, 51%] [-26%, 51%]
Posterior tail-area probability p: 0.49319
Posterior prob. of a causal effect: 51%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.00 (average absolute forecast error)
NRMSE: 0.00 (RMSE relative to pre-period mean, 21.2)
SD-ratio: 0.00 (RMSE vs pre-period variability, SD=5.9)# yss2ab<- c(as.numeric(df33$sum_velocidad_lead2), rep(NA,3))
# yss2ab<- log(yss2ab)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33$sum_velocidad_lead2), lambda_opt)
yss2ab <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
# Model 2
ss2ab <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2ab <- AddLocalLevel(ss2ab, yss2ab) #
# Add weekly seasonal
ss2ab <- AddSeasonal(ss2ab, yss2ab,
nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
#ss2ab <- AddAutoAr(ss2ab, y = c(as.numeric(df33$rate_lead2_lic), rep(NA,3)), lags = 10) #NO PUEDO AREGAR AR1 CON gaussian
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
covariates_future33b <- as.matrix(df33_post[1:4, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,
#"Y_lag10"
)])
covariates_full33b <- rbind(covariates33, covariates_future33b)
# Include regression terms
ss2ab <- AddDynamicRegression(ss2ab, covariates_full33b)
model2ab_ratio <- bsts(yss2ab,
state.specification = ss2ab, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
family ="gaussian", #A Bayesian Analysis of Time-Series Event Count Data
niter = 4e4,
#burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
seed= 2125)
## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().
impact3ab_ratio_resp <- CausalImpact(bsts.model = model2ab_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
post.period.response = as.numeric(BoxCox(as.numeric(df33_post$sum_velocidad_lead2[1:4]), lambda_opt)))
#post.period.response = log(df33_post$sum_velocidad_lead2[1:3]))
summary(impact3ab_ratio_resp)
plot(impact3ab_ratio_resp)
err_ci <- impact3ab_ratio_resp$series$response - impact3ab_ratio_resp$series$point.pred
rmse_ci <- sqrt(mean(err_ci[1:34]^2))
mean_pre <- mean(df33$sum_velocidad_lead2[1:34])
sd_pre <- sd(df33$sum_velocidad_lead2[1:34])
nrmse <- rmse_ci / mean_pre # RMSE relativo;
sd_ratio <- rmse_ci / sd_pre # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse, mean_pre))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio, sd_pre))
# Benchmark naïve (random walk de un paso):
naive_rmse <- sqrt(mean((df33$sum_velocidad_lead2[2:34] - df33$sum_velocidad_lead2[1:33])^2))
improvement <- 1 - rmse_ci / naive_rmse # >0 = mejora sobre el naïve=-=-=-=-= Iteration 0 Tue Aug 19 16:08:59 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:09:00 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:09:01 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:09:02 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:09:03 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:09:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:09:05 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:09:06 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:09:07 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:09:08 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 26 103
Prediction (s.d.) 21 (1.9) 82 (7.4)
95% CI [17, 24] [68, 98]
Absolute effect (s.d.) 5.3 (1.9) 21.1 (7.4)
95% CI [1.4, 8.8] [5.8, 35.1]
Relative effect (s.d.) 27% (12%) 27% (12%)
95% CI [5.9%, 51%] [5.9%, 51%]
Posterior tail-area probability p: 0.00585
Posterior prob. of a causal effect: 99.41528%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 1.55 (average absolute forecast error)
NRMSE: 0.03 (RMSE relative to pre-period mean, 54.6)
SD-ratio: 0.16 (RMSE vs pre-period variability, SD=9.6)# yss2cab <- c(as.numeric(df33c$sum_velocidad_lead2), rep(NA,3))
# yss2cab <- log(yss2cab)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33c$sum_velocidad_lead2), lambda_opt)
yss2cab <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
# Model 2
ss2cab <- list()
# Local trend, weekly-seasonal #https://qastack.mx/stats/209426/predictions-from-bsts-model-in-r-are-failing-completely - PUSE UN GENERALIZED LOCAL TREND
ss2cab <- AddLocalLevel(ss2cab, yss2cab) #
# Add weekly seasonal
ss2cab <- AddSeasonal(ss2cab, yss2cab,
nseasons=7, season.duration = 1) #weeks OJO, ESTOS NO SON WEEKS VERDADEROS. PORQUE TENGO MAS DE EUN AÑO
#ss2ab <- AddAutoAr(ss2ab, y = c(as.numeric(df33$rate_lead2_lic), rep(NA,3)), lags = 10) #NO PUEDO AREGAR AR1 CON gaussian
# For example, to add a day-of-week component to data with daily granularity, use model.args = list(nseasons = 7, season.duration = 1). To add a day-of-week component to data with hourly granularity, set model.args = list(nseasons = 7, season.duration = 24).
covariates_future33cb <- as.matrix(df33c_post[1:4, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,
#"Y_lag10"
)])
covariates_full33cb <- rbind(covariates33c, covariates_future33cb)
# Include regression terms
ss2cab <- AddDynamicRegression(ss2cab, covariates_full33cb)
model2cab_ratio <- bsts(yss2cab,
state.specification = ss2cab, #A list with elements created by AddLocalLinearTrend, AddSeasonal, and similar functions for adding components of state. See the help page for state.specification.
family ="gaussian", #A Bayesian Analysis of Time-Series Event Count Data
niter = 4e4,
#burn = 2000, #http://finzi.psych.upenn.edu/library/bsts/html/SuggestBurn.html Suggest the size of an MCMC burn in sample as a proportion of the total run.
seed= 2125)
## Since we have predictor variables in the model, we need to explicitly
# make their coefficients time-varying using AddDynamicRegression().
impact3cab_ratio_resp <- CausalImpact(bsts.model = model2cab_ratio,model.args = list(prior.level.sd=.1, dynamic.regression=T),
#post.period.response = log(df33c_post$sum_velocidad_lead2[1:3]))
post.period.response = as.numeric(BoxCox(as.numeric(df33c_post$sum_velocidad_lead2[1:4]), lambda_opt)))
summary(impact3cab_ratio_resp)
plot(impact3cab_ratio_resp)
err_ci_c <- impact3cab_ratio_resp$series$response - impact3cab_ratio_resp$series$point.pred
rmse_ci_c <- sqrt(mean(err_ci_c[1:34]^2))
mean_pre_c <- mean(df33c$sum_velocidad_lead2[1:34])
sd_pre_c <- sd(df33c$sum_velocidad_lead2[1:34])
nrmse_c <- rmse_ci_c / mean_pre_c # RMSE relativo;
sd_ratio_c <- rmse_ci_c / sd_pre_c # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_c))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_c, mean_pre_c))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_c, sd_pre_c))
# Benchmark naïve (random walk de un paso):
naive_rmse_c <- sqrt(mean((df33c$sum_velocidad_lead2[2:34] - df33c$sum_velocidad_lead2[1:33])^2))
improvement_c <- 1 - rmse_ci_c / naive_rmse_c # >0 = mejora sobre el naïve=-=-=-=-= Iteration 0 Tue Aug 19 16:09:12 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:09:13 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:09:14 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:09:15 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:09:16 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:09:18 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:09:19 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:09:20 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:09:21 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:09:22 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 0.74 2.97
Prediction (s.d.) 0.75 (0.098) 2.98 (0.393)
95% CI [0.55, 0.94] [2.20, 3.77]
Absolute effect (s.d.) -0.0029 (0.098) -0.0115 (0.393)
95% CI [-0.2, 0.19] [-0.8, 0.77]
Relative effect (s.d.) 1.5% (14%) 1.5% (14%)
95% CI [-21%, 35%] [-21%, 35%]
Posterior tail-area probability p: 0.48879
Posterior prob. of a causal effect: 51%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.00 (average absolute forecast error)
NRMSE: 0.00 (RMSE relative to pre-period mean, 21.2)
SD-ratio: 0.00 (RMSE vs pre-period variability, SD=5.9)# y33b <- c(as.numeric(df33$sum_velocidad_lead2), rep(NA, 3)) # respuesta escalada por 100 millones
# y33b <- log(y33b)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33$sum_velocidad_lead2), lambda_opt)
y33b <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
ss2b <- list()
ss2b <- AddLocalLinearTrend(ss2b, y33b)
ss2b <- AddSeasonal(ss2b, y33b, nseasons = 7, season.duration = 1)
#ss2b <- AddAutoAr(ss2b, y33b, lags = 1)
ss2b <- AddDynamicRegression(ss2b, covariates_full33b)
set.seed(2125)
model33b <- bsts(y33b, state.specification = ss2b, seed=2125,
niter = 4e4, family = "gaussian")
#plot(model, "components", burn = 2000)
impact3db_ratio_resp2 <- CausalImpact(bsts.model = model33b, model.args = list(prior.level.sd=.1, dynamic.regression=T),
#post.period.response = log(df33_post$sum_velocidad_lead2[1:3]))
post.period.response = as.numeric(BoxCox(as.numeric(df33_post$sum_velocidad_lead2[1:4]), lambda_opt)))
summary(impact3db_ratio_resp2)
plot(impact3db_ratio_resp2)
err_ci_corr <- impact3db_ratio_resp2$series$response - impact3db_ratio_resp2$series$point.pred
rmse_ci_corr <- sqrt(mean(err_ci_corr[1:34]^2))
mean_pre_corr <- mean(df33$sum_velocidad_lead2[1:34])
sd_pre_corr <- sd(df33$sum_velocidad_lead2[1:34])
nrmse_corr <- rmse_ci_corr / mean_pre_corr # RMSE relativo;
sd_ratio_corr <- rmse_ci_corr / sd_pre_corr # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_corr))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_corr, mean_pre_corr))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_corr, sd_pre_corr))=-=-=-=-= Iteration 0 Tue Aug 19 16:09:26 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:09:27 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:09:28 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:09:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:09:31 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:09:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:09:33 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:09:34 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:09:36 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:09:37 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 26 103
Prediction (s.d.) 21 (2.3) 85 (9.3)
95% CI [17, 26] [68, 105]
Absolute effect (s.d.) 4.5 (2.3) 18.1 (9.3)
95% CI [-0.28, 8.8] [-1.11, 35.4]
Relative effect (s.d.) 23% (14%) 23% (14%)
95% CI [-1.1%, 52%] [-1.1%, 52%]
Posterior tail-area probability p: 0.03138
Posterior prob. of a causal effect: 96.862%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 1.57 (average absolute forecast error)
NRMSE: 0.03 (RMSE relative to pre-period mean, 54.6)
SD-ratio: 0.16 (RMSE vs pre-period variability, SD=9.6)# y33b_c <- c(as.numeric(df33c$sum_velocidad_lead2), rep(NA, 3)) # respuesta escalada por 100 millones
# y33b_c <- log(y33b_c)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33c$sum_velocidad_lead2), lambda_opt)
y33b_c <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
ss2b_c <- list()
ss2b_c <- AddLocalLinearTrend(ss2b_c, y33b_c)
ss2b_c <- AddSeasonal(ss2b_c, y33b_c, nseasons = 7, season.duration = 1)
#ss2b <- AddAutoAr(ss2b, y33b, lags = 1)
covariates_future33c2 <- as.matrix(df33c_post[1:4, c("mean_min_temp_mean_lin",
"mean_max_temp_mean_lin",
"mean_median_lag_2_prec_median_imp",
"mean_median_total_precip_median_lin"#,
#"Y_lag10"
)])
covariates_full33c2 <- rbind(covariates33c, covariates_future33c2)
ss2b_c <- AddDynamicRegression(ss2b_c, covariates_full33c2)
set.seed(2125)
model33b_c <- bsts(y33b_c, state.specification = ss2b_c, seed=2125,
niter = 4e4, family = "gaussian")
#plot(model, "components", burn = 2000)
impact3db_ratio_resp2_c <- CausalImpact(bsts.model = model33b_c, model.args = list(prior.level.sd=.1, dynamic.regression=T),
#post.period.response = log(df33c_post$sum_velocidad_lead2[1:3]))
post.period.response = as.numeric(BoxCox(as.numeric(df33c_post$sum_velocidad_lead2[1:4]), lambda_opt)))
summary(impact3db_ratio_resp2_c)
plot(impact3db_ratio_resp2_c)
err_ci_c_corr <- impact3db_ratio_resp2_c$series$response - impact3db_ratio_resp2_c$series$point.pred
rmse_ci_c_corr <- sqrt(mean(err_ci_c_corr[1:34]^2))
mean_pre_c_corr <- mean(df33c$sum_velocidad_lead2[1:34])
sd_pre_c_corr <- sd(df33c$sum_velocidad_lead2[1:34])
nrmse_c_corr <- rmse_ci_c_corr / mean_pre_c_corr # RMSE relativo;
sd_ratio_c_corr <- rmse_ci_c_corr / sd_pre_c_corr # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_c_corr))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_c_corr, mean_pre_c_corr))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_c_corr, sd_pre_c_corr))=-=-=-=-= Iteration 0 Tue Aug 19 16:09:41 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:09:42 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:09:43 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:09:45 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:09:46 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:09:47 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:09:48 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:09:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:09:51 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:09:52 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 0.74 2.97
Prediction (s.d.) 0.74 (0.099) 2.98 (0.398)
95% CI [0.55, 0.94] [2.20, 3.76]
Absolute effect (s.d.) -0.0022 (0.099) -0.0090 (0.398)
95% CI [-0.20, 0.19] [-0.79, 0.77]
Relative effect (s.d.) 1.6% (15%) 1.6% (15%)
95% CI [-21%, 35%] [-21%, 35%]
Posterior tail-area probability p: 0.49158
Posterior prob. of a causal effect: 51%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.00 (average absolute forecast error)
NRMSE: 0.00 (RMSE relative to pre-period mean, 21.2)
SD-ratio: 0.00 (RMSE vs pre-period variability, SD=5.9)## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33_qc <- c(as.numeric(df33_qc$sum_velocidad_lead2), rep(NA, 3)) # 7 post-días + 2 de ajuste
# y33_qc <- log(y33_qc)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33_qc$sum_velocidad_lead2), lambda_opt)
y33_qc <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates33_qc <- df33_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin#,
#Y_lag10
)
covariates_future33_qc <- df33_post_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin#,
#Y_lag10
)
covariates_full33_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:7,])
covariates_full33_qc2 <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:4,])## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33_qc_c <- c(as.numeric(df33c_qc$sum_velocidad_lead2), rep(NA, 3)) # 7 post-días + 2 de ajuste
# y33_qc_c <- log(y33_qc_c)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- BoxCox(as.numeric(df33c_qc$sum_velocidad_lead2), lambda_opt)
y33_qc_c <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates33_qc_c <- df33c_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin#,
#Y_lag10
)
covariates_future33_qc_c <- df33c_post_qc |> dplyr::select(
mean_min_temp_mean_lin,
mean_max_temp_mean_lin,
mean_median_lag_2_prec_median_imp,
mean_median_total_precip_median_lin#,
#Y_lag10
)
covariates_full33_qc_c <- dplyr::bind_rows(covariates33_qc_c, covariates_future33_qc_c[1:7,])
covariates_full33_qc_c2 <- dplyr::bind_rows(covariates33_qc_c, covariates_future33_qc_c[1:4,])## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33b_qc <- c(as.numeric(df33_qc$sum_velocidad_lead2), rep(NA, 3)) # 7 post-días + 2 de ajuste
# y33b_qc <- log(y33b_qc)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y_bc <- forecast::BoxCox(as.numeric(df33_qc$sum_velocidad_lead2), lambda_opt)
y33b_qc <- c(y_bc, rep(NA, 4)) # respuesta escalada por 100 millones
covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc[1:4,])
## ── Especificación del estado: nivel local ───────────────────────────────────
ss2ab_qc <- list()
ss2ab_qc <- AddLocalLevel(ss2ab_qc, y33b_qc)
ss2ab_qc <- AddSeasonal(ss2ab_qc, y33b_qc, nseasons = 7, season.duration = 1)
#ss2ab_qc <- AddAutoAr (ss2ab_qc, y33b_qc, lags = 10)
ss2ab_qc <- AddDynamicRegression(ss2ab_qc, covariates_full33b_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2ab_qc <- bsts(
y33b_qc,
state.specification = ss2ab_qc,
family = "gaussian",
niter = 4e4,
seed = 2125
)
impact33b_qc <- CausalImpact(
bsts.model = model2ab_qc,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
#post.period.response = log(df33_post_qc$sum_velocidad_lead2[1:3])
post.period.response = as.numeric(forecast::BoxCox(as.numeric(df33_post_qc$sum_velocidad_lead2[1:4]), lambda_opt))
)
summary(impact33b_qc)
err_ci_qb_corr <- impact33b_qc$series$response - impact33b_qc$series$point.pred
rmse_ci_qb_corr <- sqrt(mean(err_ci_qb_corr[1:34]^2))
mean_pre_qb_corr <- mean(df33_qc$sum_velocidad_lead2[1:34], na.rm=T)
sd_pre_qb_corr <- sd(df33_qc$sum_velocidad_lead2[1:34], na.rm=T)
nrmse_qb_corr <- rmse_ci_qb_corr / mean_pre_qb_corr # RMSE relativo;
sd_ratio_qb_corr <- rmse_ci_qb_corr / sd_pre_qb_corr # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_qb_corr))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_qb_corr, mean_pre_qb_corr))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_qb_corr, sd_pre_qb_corr))
plot(impact33b_qc)
=-=-=-=-= Iteration 0 Tue Aug 19 16:09:57 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:09:58 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:09:59 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:10:00 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:10:02 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:10:03 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:10:04 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:10:05 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:10:06 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:10:07 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 15 61
Prediction (s.d.) 13 (0.95) 52 (3.82)
95% CI [11, 15] [44, 59]
Absolute effect (s.d.) 2.3 (0.95) 9.1 (3.82)
95% CI [0.36, 4.1] [1.44, 16.5]
Relative effect (s.d.) 18% (8.9%) 18% (8.9%)
95% CI [2.4%, 37%] [2.4%, 37%]
Posterior tail-area probability p: 0.01214
Posterior prob. of a causal effect: 98.786%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.60 (average absolute forecast error)
NRMSE: 0.01 (RMSE relative to pre-period mean, 56.7)
SD-ratio: 0.06 (RMSE vs pre-period variability, SD=9.5)## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33b_qc_c <- c(as.numeric(df33c_qc$sum_velocidad_lead2), rep(NA, 3)) # 7 post-días + 2 de ajuste
# y33b_qc_c <- log(y33b_qc_c)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y33b_qc_c <- BoxCox(as.numeric(df33c_qc$sum_velocidad_lead2), lambda_opt)
y33b_qc_c <- c(as.numeric(y33b_qc_c), rep(NA, 4))
## ── Especificación del estado: nivel local ───────────────────────────────────
ss2ab_qc_c <- list()
ss2ab_qc_c <- AddLocalLevel(ss2ab_qc_c, y33b_qc_c)
ss2ab_qc_c <- AddSeasonal(ss2ab_qc_c, y33b_qc_c, nseasons = 7, season.duration = 1)
#ss2ab_qc <- AddAutoAr (ss2ab_qc, y33b_qc, lags = 10)
ss2ab_qc_c <- AddDynamicRegression(ss2ab_qc_c, covariates_full33_qc_c2)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model2ab_qc_c <- bsts(
y33b_qc_c,
state.specification = ss2ab_qc_c,
family = "gaussian",
niter = 4e4,
seed = 2125
)
impact33ab_qc_c <- CausalImpact(
bsts.model = model2ab_qc_c,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
#post.period.response = log(df33c_post_qc$sum_velocidad_lead2[1:3])
post.period.response = as.numeric(BoxCox(as.numeric(df33c_post_qc$sum_velocidad_lead2[1:4]), lambda_opt))
)
summary(impact33ab_qc_c)
err_ci_c_qb_corr <- impact33ab_qc_c$series$response - impact33ab_qc_c$series$point.pred
rmse_ci_c_qb_corr <- sqrt(mean(err_ci_c_qb_corr[1:34]^2))
mean_pre_c_qb_corr <- mean(df33c_qc$sum_velocidad_lead2[1:34], na.rm=T)
sd_pre_c_qb_corr <- sd(df33c_qc$sum_velocidad_lead2[1:34], na.rm=T)
nrmse_c_qb_corr <- rmse_ci_c_qb_corr / mean_pre_c_qb_corr # RMSE relativo;
sd_ratio_c_qb_corr <- rmse_ci_c_qb_corr / sd_pre_c_qb_corr # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_c_qb_corr))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_c_qb_corr, mean_pre_c_qb_corr))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_c_qb_corr, sd_pre_c_qb_corr))
plot(impact33ab_qc_c)
=-=-=-=-= Iteration 0 Tue Aug 19 16:10:11 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:10:12 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:10:13 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:10:14 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:10:15 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:10:16 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:10:17 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:10:18 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:10:19 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:10:20 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 0.58 2.31
Prediction (s.d.) 0.58 (0.091) 2.32 (0.362)
95% CI [0.4, 0.76] [1.6, 3.03]
Absolute effect (s.d.) -0.0017 (0.091) -0.0070 (0.362)
95% CI [-0.18, 0.18] [-0.72, 0.71]
Relative effect (s.d.) 2.4% (18%) 2.4% (18%)
95% CI [-24%, 44%] [-24%, 44%]
Posterior tail-area probability p: 0.49217
Posterior prob. of a causal effect: 51%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.00 (average absolute forecast error)
NRMSE: 0.00 (RMSE relative to pre-period mean, 19.0)
SD-ratio: 0.00 (RMSE vs pre-period variability, SD=5.2)## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33b_lin_qc <- c(as.numeric(df33_qc$sum_velocidad_lead2), rep(NA, 3))
# y33b_lin_qc <- log(y33b_lin_qc)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y33b_lin_qc <- BoxCox(as.numeric(df33_qc$sum_velocidad_lead2), lambda_opt)
y33b_lin_qc <- c(as.numeric(y33b_lin_qc), rep(NA, 4))
## ── Covariables (data.frame) ─────────────────────────────────────────────────
covariates_full33b_qc <- dplyr::bind_rows(covariates33_qc, covariates_future33_qc)
## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2bb_qc <- list()
ss2bb_qc <- AddLocalLinearTrend(ss2bb_qc, y33b_lin_qc)
ss2bb_qc <- AddSeasonal (ss2bb_qc, y33b_lin_qc, nseasons = 7, season.duration = 1)
#ss2bb_qc <- AddAutoAr (ss2bb_qc, y33b_lin_qc, lags = 1)
ss2bb_qc <- AddDynamicRegression(ss2bb_qc, covariates_full33b_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc2 <- bsts(
y33b_lin_qc,
state.specification = ss2bb_qc,
family = "gaussian",
niter = 4e4,
seed = 2125
)
impact2_33b_qc2 <- CausalImpact(
bsts.model = model33b_qc2,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
#post.period.response = log(df33_post_qc$sum_velocidad_lead2[1:3]),
post.period.response = as.numeric(forecast::BoxCox(as.numeric(df33_post_qc$sum_velocidad_lead2[1:4]), lambda_opt))
)
summary(impact2_33b_qc2)
err_ci_qb_corr2 <- impact2_33b_qc2$series$response - impact2_33b_qc2$series$point.pred
rmse_ci_qb_corr2 <- sqrt(mean(err_ci_qb_corr2[1:34]^2))
mean_pre_qb_corr2 <- mean(df33_qc$sum_velocidad_lead2[1:34], na.rm=T)
sd_pre_qb_corr2 <- sd(df33_qc$sum_velocidad_lead2[1:34], na.rm=T)
nrmse_qb_corr2 <- rmse_ci_qb_corr2 / mean_pre_qb_corr2 # RMSE relativo;
sd_ratio_qb_corr2 <- rmse_ci_qb_corr2 / sd_pre_qb_corr2 # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_qb_corr2))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_qb_corr2, mean_pre_qb_corr2))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_qb_corr2, sd_pre_qb_corr2))
plot(impact2_33b_qc2)
=-=-=-=-= Iteration 0 Tue Aug 19 16:10:25 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:10:26 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:10:27 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:10:28 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:10:30 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:10:31 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:10:32 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:10:33 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:10:35 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:10:36 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 15 61
Prediction (s.d.) 13 (1.2) 52 (4.6)
95% CI [11, 15] [43, 62]
Absolute effect (s.d.) 2.1 (1.2) 8.5 (4.6)
95% CI [-0.21, 4.4] [-0.86, 17.4]
Relative effect (s.d.) 17% (11%) 17% (11%)
95% CI [-1.4%, 40%] [-1.4%, 40%]
Posterior tail-area probability p: 0.03546
Posterior prob. of a causal effect: 96.454%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.62 (average absolute forecast error)
NRMSE: 0.01 (RMSE relative to pre-period mean, 56.7)
SD-ratio: 0.07 (RMSE vs pre-period variability, SD=9.5)## ── Vector de respuesta ───────────────────────────────────────────────────────
# y33b_lin_qc_c <- c(as.numeric(df33c_qc$sum_velocidad_lead2), rep(NA, 3))
# y33b_lin_qc_c <- log(y33b_lin_qc_c)
lambda_opt <- forecast::BoxCox.lambda(as.numeric(df33c_qc$sum_velocidad_lead2), lower = -2, upper = 2)
y33b_lin_qc_c <- BoxCox(as.numeric(df33c_qc$sum_velocidad_lead2), lambda_opt)
y33b_lin_qc_c <- c(as.numeric(y33b_lin_qc_c), rep(NA, 4))
## ── Estado: tendencia local lineal ───────────────────────────────────────────
ss2bb_qc_c <- list()
ss2bb_qc_c <- AddLocalLinearTrend(ss2bb_qc_c, y33b_lin_qc_c)
ss2bb_qc_c <- AddSeasonal (ss2bb_qc_c, y33b_lin_qc_c, nseasons = 7, season.duration = 1)
#ss2bb_qc <- AddAutoAr (ss2bb_qc, y33b_lin_qc, lags = 1)
ss2bb_qc_c <- AddDynamicRegression(ss2bb_qc_c, covariates_full33b_qc)
## ── Ajuste BSTS y CausalImpact ───────────────────────────────────────────────
model33b_qc_c <- bsts(
y33b_lin_qc_c,
state.specification = ss2bb_qc_c,
family = "gaussian",
niter = 4e4,
seed = 2125
)
impact2_33b_qc_c <- CausalImpact(
bsts.model = model33b_qc_c,
model.args = list(prior.level.sd = .1, dynamic.regression = TRUE),
#post.period.response = log(df33c_post_qc$sum_velocidad_lead2[1:3])
post.period.response = as.numeric(BoxCox(as.numeric(df33c_post_qc$sum_velocidad_lead2[1:4]), lambda_opt))
)
summary(impact2_33b_qc_c)
err_ci_c_qb_corr2 <- impact2_33b_qc_c$series$response - impact2_33b_qc_c$series$point.pred
rmse_ci_c_qb_corr2 <- sqrt(mean(err_ci_c_qb_corr2[1:34]^2))
mean_pre_c_qb_corr2 <- mean(df33c_qc$sum_velocidad_lead2[1:34], na.rm=T)
sd_pre_c_qb_corr2 <- sd(df33c_qc$sum_velocidad_lead2[1:34], na.rm=T)
nrmse_c_qb_corr2 <- rmse_ci_c_qb_corr2 / mean_pre_c_qb_corr2 # RMSE relativo;
sd_ratio_c_qb_corr2 <- rmse_ci_c_qb_corr2 / sd_pre_c_qb_corr2 # RMSE relativo a la variabilidad natural
cat(sprintf("RMSE on 35 post-period points: %.2f (average absolute forecast error)\n", rmse_ci_c_qb_corr2))
cat(sprintf("NRMSE: %.2f (RMSE relative to pre-period mean, %.1f)\n", nrmse_c_qb_corr2, mean_pre_c_qb_corr2))
cat(sprintf("SD-ratio: %.2f (RMSE vs pre-period variability, SD=%.1f)\n", sd_ratio_c_qb_corr2, sd_pre_c_qb_corr2))
plot(impact2_33b_qc_c)
=-=-=-=-= Iteration 0 Tue Aug 19 16:10:40 2025
=-=-=-=-=
=-=-=-=-= Iteration 4000 Tue Aug 19 16:10:41 2025
=-=-=-=-=
=-=-=-=-= Iteration 8000 Tue Aug 19 16:10:43 2025
=-=-=-=-=
=-=-=-=-= Iteration 12000 Tue Aug 19 16:10:44 2025
=-=-=-=-=
=-=-=-=-= Iteration 16000 Tue Aug 19 16:10:45 2025
=-=-=-=-=
=-=-=-=-= Iteration 20000 Tue Aug 19 16:10:46 2025
=-=-=-=-=
=-=-=-=-= Iteration 24000 Tue Aug 19 16:10:47 2025
=-=-=-=-=
=-=-=-=-= Iteration 28000 Tue Aug 19 16:10:49 2025
=-=-=-=-=
=-=-=-=-= Iteration 32000 Tue Aug 19 16:10:50 2025
=-=-=-=-=
=-=-=-=-= Iteration 36000 Tue Aug 19 16:10:51 2025
=-=-=-=-=
Posterior inference {CausalImpact}
Average Cumulative
Actual 0.58 2.31
Prediction (s.d.) 0.58 (0.097) 2.31 (0.390)
95% CI [0.39, 0.77] [1.55, 3.09]
Absolute effect (s.d.) -0.00089 (0.097) -0.00357 (0.390)
95% CI [-0.19, 0.19] [-0.78, 0.76]
Relative effect (s.d.) 3% (20%) 3% (20%)
95% CI [-25%, 49%] [-25%, 49%]
Posterior tail-area probability p: 0.49647
Posterior prob. of a causal effect: 50%
For more details, type: summary(impact, "report")
RMSE on 35 post-period points: 0.00 (average absolute forecast error)
NRMSE: 0.00 (RMSE relative to pre-period mean, 19.0)
SD-ratio: 0.00 (RMSE vs pre-period variability, SD=5.2)# Inversa de Box–Cox (con salvaguarda numérica)
inv_boxcox <- function(z, lambda) {
if (abs(lambda) < 1e-12) return(exp(z)) # caso λ≈0 (log)
base <- lambda * z + 1
base[base <= 0] <- .Machine$double.eps # evita NaN por redondeo/colas
base^(1 / lambda)
}
# RR a partir de un objeto CausalImpact usando IC para obtener sd (post = últimos Tpost)
rr_from_impact <- function(impact, lambda, Tpost = 4, n_sims = 20000, level = 0.95, seed = 2125) {
stopifnot(!is.null(impact$series))
ser <- as.data.frame(impact$series)
post <- tail(ser, Tpost)
req <- c("response", "point.pred", "point.pred.lower", "point.pred.upper")
if (!all(req %in% names(post))) {
stop("Faltan columnas en impact$series: se requieren 'response', 'point.pred', 'point.pred.lower', 'point.pred.upper'.")
}
z <- qnorm(0.5 + level/2)
y_t <- as.numeric(post$response) # observado (escala Box–Cox)
mu_t <- as.numeric(post$point.pred) # media contrafactual (escala Box–Cox)
sd_t <- as.numeric(post$point.pred.upper - post$point.pred.lower) / (2*z)
# Observado a escala original
y <- inv_boxcox(y_t, lambda)
# Simulación del contrafactual en escala Box–Cox y back-transform
set.seed(seed)
Tpost <- length(mu_t)
cf_t <- matrix(
rnorm(Tpost * n_sims, mean = rep(mu_t, n_sims), sd = rep(sd_t, n_sims)),
nrow = Tpost, ncol = n_sims, byrow = FALSE
)
cf <- inv_boxcox(cf_t, lambda)
# RR puntual (usando la media contrafactual) y RR por simulación
rr_hat <- sum(y) / sum(inv_boxcox(mu_t, lambda))
rr_draws <- sum(y) / pmax(colSums(cf), .Machine$double.eps)
list(rr_hat = rr_hat, rr_draws = rr_draws)
}
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:#:
## --- Asigna los λ que realmente usaste en cada impacto ---
lambda_imp1 <- 0.6779003
lambda_imp2 <- -1.325838
# Calcula RR de cada impacto (post = 4 tiempos)
res1 <- rr_from_impact(impact3db_ratio_resp2, lambda = lambda_imp1, Tpost = 4, n_sims = 2e5, seed = 2125)
sprintf("Back-transformed RRs: %.2f (%.2f,%.2f)", quantile(res1$rr_draws, c(0.5), na.rm = TRUE),quantile(res1$rr_draws, c(0.025), na.rm = TRUE),quantile(res1$rr_draws, c(0.975), na.rm = TRUE))[1] "Back-transformed RRs: 1.30 (1.07,1.63)"res2 <- rr_from_impact(impact3db_ratio_resp2_c, lambda = lambda_imp2, Tpost = 4, n_sims = 2e5, seed = 2125)
sprintf("Back-transformed RRs: %.2f (%.2f,%.2f)", quantile(res2$rr_draws, c(0.5), na.rm = TRUE),quantile(res2$rr_draws, c(0.025), na.rm = TRUE),quantile(res2$rr_draws, c(0.975), na.rm = TRUE))[1] "Back-transformed RRs: 0.00 (0.00,3.89)"# RRR e intervalo (cuantílico)
rrr_draws <- res1$rr_draws / res2$rr_draws
RRR_hat <- res1$rr_hat / res2$rr_hat
ci <- quantile(rrr_draws, c(0.025, 0.5, 0.975), na.rm = TRUE)
cat("Québec only, -3 días pre, local linear trend (escala original via back-transform)\n")Québec only, -3 días pre, local linear trend (escala original via back-transform)cat(sprintf("RRR = %.3f (95%% CI %.3f–%.3f)\n", RRR_hat, ci[1], ci[3]))RRR = 1.535 (95% CI 0.404–29841960211.971)#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_#_
## --- Asigna los λ que realmente usaste en cada impacto ---
lambda_imp1q <- 0.4903961 # para impact33b_qc (df33_qc)
lambda_imp2q <- -1.720255 # para impact2_33b_qc_c (df33c_qc)
# Calcula RR de cada impacto (post = 4 tiempos)
res1q <- rr_from_impact(impact33b_qc, lambda = lambda_imp1q, Tpost = 4, n_sims = 2e5, seed = 2125)
sprintf("Back-transformed RRs: %.2f (%.2f,%.2f)", quantile(res1q$rr_draws, c(0.5), na.rm = TRUE),quantile(res1q$rr_draws, c(0.025), na.rm = TRUE),quantile(res1q$rr_draws, c(0.975), na.rm = TRUE))[1] "Back-transformed RRs: 1.32 (1.10,1.61)"res2q <- rr_from_impact(impact2_33b_qc_c, lambda = lambda_imp2q, Tpost = 4, n_sims = 2e5, seed = 2125)
sprintf("Back-transformed RRs: %.2f (%.2f,%.2f)", quantile(res2q$rr_draws, c(0.5), na.rm = TRUE),quantile(res2q$rr_draws, c(0.025), na.rm = TRUE),quantile(res2q$rr_draws, c(0.975), na.rm = TRUE))[1] "Back-transformed RRs: 0.00 (0.00,5.50)"# RRR e intervalo (cuantílico)
rrr_drawsq <- res1q$rr_draws / res2q$rr_draws
RRR_hatq <- res1q$rr_hat / res2q$rr_hat
ciq <- quantile(rrr_drawsq, c(0.025, 0.5, 0.975), na.rm = TRUE)
cat("Québec only, -3 días pre, local linear trend (escala original via back-transform)\n")Québec only, -3 días pre, local linear trend (escala original via back-transform)cat(sprintf("RRR = %.3f (95%% CI %.3f–%.3f)\n", RRR_hatq, ciq[1], ciq[3]))RRR = 1.568 (95% CI 0.285–73843935.211)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: 2025-08-19 16:10:56.225483Editor 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/Rtmp2ZIvTs/file4600278e614f -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 |
| 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 |
| cowplot | 1.1.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 |
| devtools | 2.4.5 | RSPM |
| DHARMa | 0.4.7 | 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 |
| 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 |
| ggridges | 0.5.6 | 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 |
| 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 |
| inline | 0.3.21 | 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 |
| 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 |
| 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 |
| 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 |
| 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 |
| 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 |
| 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 |
| 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 |
| 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 |
| Synth | 1.1-8 | RSPM |
| systemfonts | 1.2.3 | 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 |
| utf8 | 1.2.5 | RSPM |
| V8 | 6.0.3 | RSPM |
| vctrs | 0.6.5 | RSPM |
| viridisLite | 0.4.2 | RSPM |
| weathercan | 0.7.3.9000 | https://ropensci.r-universe.dev (R 4.4.3) |
| 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