Clima y resultados perinatales: reproducción autocontenida del análisis
R Markdown reproducible
23/08/2026
1 1. Objetivo y principio de reproducibilidad
Este documento contiene todo el código necesario para ejecutar el análisis dentro de una sola sesión de R.
No utiliza source().
No utiliza resultados previamente generados.
No lee archivos Results_*.csv ni figuras históricas.
Las tablas, Odds Ratios, intervalos de confianza y gráficas que aparecen
al final son calculados nuevamente durante el Knit.
El flujo es:
Datos_clima_final.rds
↓
limpieza de la cohorte
↓
desenlaces perinatales
↓
case-crossover tiempo-estratificado
↓
extracción espacial de clima
↓
TEMP / WBGT / UTCI / Heat Index
↓
lags 0–6
↓
DLNM: crossbasis()
↓
regresión logística condicional: clogit()
↓
crosspred()
↓
OR + IC95% + curvas exposición-respuesta2 2. Decisiones analíticas reproducidas
El análisis principal utiliza:
- periodo: 2012–2021;
- diseño: time-stratified case-crossover;
- días control: mismo año, mismo mes y mismo día de la semana;
- rezagos: 0–6 días;
- exposición principal: MAX UTCI;
- función exposición-respuesta: natural spline;
- knots de exposición: percentiles 10, 75 y 90;
- función de lag: natural spline con 2 knots logarítmicos;
- referencia: percentil 75;
- comparación alta: percentil 99;
- modelo: regresión logística condicional con
strata(PTID); - desenlaces: PTB, muerte fetal, Apgar bajo, cesárea y compuesto.
3 3. Paquetes
if (!requireNamespace("pacman", quietly = TRUE)) {
install.packages("pacman")
}
pacman::p_load(
here,
dplyr,
tidyr,
purrr,
tibble,
stringr,
lubridate,
data.table,
raster,
sp,
ncdf4,
survival,
dlnm,
mixmeta,
splines,
ggplot2,
scales,
DT,
htmltools
)4 3.1 Evitar conflictos
entre raster y dplyr
El paquete raster también define una función llamada
select(). Para evitar que R use accidentalmente
raster::select() sobre data frames, este documento utiliza
explícitamente dplyr::select(),
dplyr::filter(), dplyr::rename() y otros
verbos de manipulación de datos.
5 4. Verificación de la raíz del proyecto
El Rmd debe estar guardado en la carpeta principal
Clima, al mismo nivel que
Datos_clima_final.rds y Clean.
cat("Directorio de trabajo:\n")## Directorio de trabajo:print(getwd())## [1] "C:/Users/GERMAN/Dropbox/Clima (1)"cat("\nRaíz de here():\n")##
## Raíz de here():print(here::here())## [1] "C:/Users/GERMAN/Dropbox/Clima (1)"rutas_necesarias <- tibble::tibble(
Recurso = c(
"Base perinatal",
"Temperature",
"WBGT",
"UTCI",
"Heat Index"
),
Ruta = c(
here::here("Datos_clima_final.rds"),
here::here("Clean", "Temperature"),
here::here("Clean", "WBGT"),
here::here("Clean", "UTCI"),
here::here("Clean", "Heat_Index")
)
) %>%
dplyr::mutate(
Existe = c(
file.exists(Ruta[1]),
dir.exists(Ruta[-1])
)
)
rutas_necesariasif (!all(rutas_necesarias$Existe)) {
stop("Faltan archivos o carpetas necesarias. Revise la tabla anterior.")
}6 5. Auditoría del entorno computacional
El ejercicio original fue desarrollado con R 4.4.1 y
dlnm 2.4.7.
Aquí se registra la versión utilizada en esta reproducción, pero
no se detiene el análisis si es diferente.
auditoria_entorno <- tibble::tibble(
Componente = c("R", "dlnm", "survival"),
Version_referencia = c(
"4.4.1",
"2.4.7",
NA_character_
),
Version_actual = c(
paste0(R.version$major, ".", R.version$minor),
as.character(packageVersion("dlnm")),
as.character(packageVersion("survival"))
)
) %>%
dplyr::mutate(
Coincide = ifelse(
is.na(Version_referencia),
NA,
Version_referencia == Version_actual
)
)
auditoria_entorno7 6. Parámetros del análisis
PAISES <- c(
ARG = "Argentina",
BOL = "Bolivia",
GTM = "Guatemala",
HND = "Honduras",
DOM = "Dominican Republic",
URY = "Uruguay"
)
ANIOS <- 2012:2021
LAG_MAX <- 6
STAT_MAIN <- "MAX"
EXPOSURE_MAIN <- "UTCI"
ARGVAR_KNOTS <- c(10, 75, 90)
ARGLAG_NK <- 2
CENTER_PERCENTILE <- 0.75
COMPARATOR_PERCENTILE <- 0.99
RUN_SENSITIVITY <- TRUE
PAISES## ARG BOL GTM
## "Argentina" "Bolivia" "Guatemala"
## HND DOM URY
## "Honduras" "Dominican Republic" "Uruguay"8 7. Leer la base principal
df_raw <- readRDS(
here::here("Datos_clima_final.rds")
)
cat("Filas:", format(nrow(df_raw), big.mark = "."), "\n")## Warning in prettyNum(.Internal(format(x, trim, digits, nsmall, width, 3L, :
## 'big.mark' y 'decimal.mark' son ambos '.', lo cual puede ser confuso## Filas: 201.636cat("Columnas:", ncol(df_raw), "\n")## Columnas: 45dplyr::glimpse(df_raw)## Rows: 201,636
## Columns: 45
## $ Pais <chr> "Rep.Dom", "Rep.Dom", "Rep.Dom", "Rep.Dom", "Rep.Dom"…
## $ Institucion <chr> "AltaGraciaDO", "AltaGraciaDO", "AltaGraciaDO", "Alta…
## $ VAR_0284 <chr> "2021-06-13", "2022-02-28", "2022-01-13", "2021-02-13…
## $ VAR_0009 <chr> "27", "25", "35", "21", "32", "30", "31", "36", "28",…
## $ VAR_0010 <chr> NA, NA, NA, NA, NA, NA, NA, "X", NA, NA, NA, NA, NA, …
## $ Anio <dbl> 2021, 2022, 2022, 2021, 2021, 2020, 2021, 2022, 2022,…
## $ VAR_0006 <chr> "1994-06-18", "1996-03-09", "1986-02-05", "1999-05-13…
## $ Longuitud <dbl> -69.90681, -69.90681, -69.90681, -69.90681, -69.90681…
## $ Latitud <dbl> 18.47343, 18.47343, 18.47343, 18.47343, 18.47343, 18.…
## $ VAR_0282 <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A"…
## $ VAR_0310 <chr> "B", "B", "A", "A", "A", "B", "A", "A", "B", "A", "A"…
## $ VAR_0311 <chr> "2892", "3629", "2948", "3175", "3232", "3459", "4309…
## $ VAR_0312 <chr> NA, NA, NA, NA, NA, NA, "B", "A", NA, NA, NA, NA, NA,…
## $ VAR_0315 <chr> "38", "41", "37", "39", "39", "39", "40", "36", "38",…
## $ VAR_0316 <chr> NA, NA, NA, "0", NA, NA, "0", NA, "0", NA, NA, "0", "…
## $ VAR_0320 <chr> "A", "A", "A", "A", NA, NA, "A", "A", "A", "A", "A", …
## $ VAR_0321 <chr> "8", "8", "8", "8", "8", "8", "8", "7", NA, "8", "8",…
## $ VAR_0322 <chr> "9", "9", "9", "9", "9", "9", "9", "8", NA, "9", "9",…
## $ VAR_0371 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "A", NA, …
## $ VAR_0329 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, "A", "A", "A", NA…
## $ VAR_0702 <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ VAR_0055 <chr> "0", "0", "140", "0", "0", "0", "196", "90", NA, NA, …
## $ VAR_0056 <chr> "65", "66", "72", "62", "60", "56", "80", "62", NA, N…
## $ VAR_0011 <chr> "A", "C", "C", "C", "D", "D", "D", "C", NA, "C", "C",…
## $ VAR_0012 <chr> "B", "B", "B", "B", "B", "B", "B", "B", NA, "B", NA, …
## $ VAR_0013 <chr> "D", "D", "D", "B", "C", "B", "C", "C", NA, "C", NA, …
## $ VAR_0014 <chr> "4", "5", "0", "8", "2", "6", "2", "4", NA, NA, NA, N…
## $ VAR_0015 <chr> "A", "A", "A", "B", NA, "B", "B", "B", NA, "B", NA, N…
## $ VAR_0016 <chr> "A", "A", "A", "A", "A", NA, NA, "A", NA, NA, NA, NA,…
## $ VAR_0047 <chr> "0", "0", "0", "2", "0", "3", "0", "2", "1", "0", "1"…
## $ VAR_0287 <chr> "B", "A", "B", "B", "B", "B", "A", "B", "A", "A", "B"…
## $ VAR_0040 <chr> "1", "1", "4", "3", "1", "4", "2", "3", "4", "5", "2"…
## $ edad_materna <dbl> 26, 25, 35, 21, 33, 29, 31, 36, NA, 30, NA, 29, 18, 2…
## $ diff_edad <dbl> -1, 0, 0, 0, 1, -1, 0, 0, NA, 0, NA, 2, 0, 0, NA, 1, …
## $ edad_materna_cat <chr> "25-35", "25-35", "25-35", "<25", "25-35", "25-35", "…
## $ .default <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ VAR_0009_cat <chr> "25-35", "25-35", "25-35", "<25", "25-35", "25-35", "…
## $ edad_final_cat <chr> "25-35", "25-35", "25-35", "<25", "25-35", "25-35", "…
## $ M_AGE <chr> "25-35", "25-35", "25-35", "<25", "25-35", "25-35", "…
## $ Pretermino <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,…
## $ sb <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ apgar5 <dbl> 9, 9, 9, 9, 9, 9, 9, 8, NA, 9, 9, 9, 9, NA, 9, 9, 9, …
## $ sexo <chr> "Male", "Male", "Female", "Female", "Female", "Male",…
## $ LATITUDE <dbl> 18.47343, 18.47343, 18.47343, 18.47343, 18.47343, 18.…
## $ LONGITUDE <dbl> -69.90681, -69.90681, -69.90681, -69.90681, -69.90681…9 8. Caracterización descriptiva inicial por país
Esta sección no transforma la base usada posteriormente; describe la fuente original.
df_descriptivo <- df_raw %>%
dplyr::mutate(
Pais = ifelse(
Pais == "Rep.Dom",
"Dominican Republic",
as.character(Pais)
)
)
tabla_pais <- df_descriptivo %>%
dplyr::count(Pais, name = "N") %>%
dplyr::mutate(
Porcentaje = round(100 * N / sum(N), 2)
) %>%
dplyr::arrange(desc(N))
DT::datatable(
tabla_pais,
rownames = FALSE,
options = list(scrollX = TRUE)
)ggplot(
tabla_pais,
aes(
x = reorder(Pais, N),
y = N
)
) +
geom_col() +
coord_flip() +
scale_y_continuous(labels = comma) +
labs(
title = "Registros disponibles por país",
x = NULL,
y = "Número de registros"
) +
theme_minimal()
10 9. Funciones auxiliares para los datos climáticos
10.1 9.1 Normalización de coordenadas
normalize_coords <- function(coords) {
x <- as.data.frame(coords)
if (ncol(x) < 2) {
stop("coords debe tener al menos dos columnas.")
}
nms <- tolower(names(x))
lon_idx <- which(
nms %in% c("lon", "longitude", "x")
)
lat_idx <- which(
nms %in% c("lat", "latitude", "y")
)
if (
length(lon_idx) == 1 &&
length(lat_idx) == 1
) {
x <- x[, c(lon_idx, lat_idx), drop = FALSE]
} else {
x <- x[, 1:2, drop = FALSE]
}
names(x) <- c("lon", "lat")
x$lon <- as.numeric(
as.character(x$lon)
)
x$lat <- as.numeric(
as.character(x$lat)
)
x
}10.2 9.2 Ajuste mínimo para Uruguay
clamp_uruguay <- function(r, coords, country_code, eps = 1e-6) {
coords <- normalize_coords(coords)
if (country_code != "URY") {
return(coords)
}
ex <- raster::extent(r)
coords$lon <- pmin(
pmax(coords$lon, ex@xmin + eps),
ex@xmax - eps
)
coords$lat <- pmin(
pmax(coords$lat, ex@ymin + eps),
ex@ymax - eps
)
coords
}10.3 9.3 Fechas de las capas raster
Primero se intenta interpretar las fechas desde los nombres de las
capas.
Si no es posible, se obtiene el año del nombre del archivo y se
construye una secuencia diaria desde el 1 de enero.
dates_from_raster <- function(r, file) {
layer_names <- names(r)
parse_layer_date <- function(z) {
z <- sub("^X", "", z)
candidatos <- c(
"%Y.%m.%d",
"%Y-%m-%d",
"%Y_%m_%d"
)
for (fmt in candidatos) {
ans <- suppressWarnings(
as.Date(z, format = fmt)
)
if (!is.na(ans)) {
return(ans)
}
}
as.Date(NA)
}
fechas <- as.Date(
sapply(layer_names, parse_layer_date),
origin = "1970-01-01"
)
if (
length(fechas) == raster::nlayers(r) &&
all(!is.na(fechas))
) {
return(fechas)
}
year <- stringr::str_extract(
basename(file),
"\\d{4}"
)
year <- as.integer(year)
if (is.na(year)) {
stop("No fue posible identificar el año en: ", file)
}
as.Date(
paste0(year, "-01-01")
) + seq_len(raster::nlayers(r)) - 1
}10.4 9.4 Verificación de unidades
nc_units <- function(file) {
nc <- try(
ncdf4::nc_open(file),
silent = TRUE
)
if (inherits(nc, "try-error")) {
return(NA_character_)
}
on.exit(
try(ncdf4::nc_close(nc), silent = TRUE),
add = TRUE
)
vars <- names(nc$var)
if (length(vars) == 0) {
return(NA_character_)
}
as.character(
nc$var[[vars[1]]]$units
)
}
convert_to_celsius_if_needed <- function(values, file, exposure) {
if (
exposure %in% c(
"WBGT",
"HEAT_IND"
)
) {
return(values)
}
units <- tolower(
trimws(
nc_units(file)
)
)
kelvin_by_units <- !is.na(units) &&
(
units == "k" ||
grepl("kelvin", units)
)
med <- suppressWarnings(
median(
as.matrix(values),
na.rm = TRUE
)
)
kelvin_by_values <- is.finite(med) &&
med > 100
if (
kelvin_by_units ||
kelvin_by_values
) {
values[] <- lapply(
values,
function(x) x - 273.15
)
}
values
}10.5 9.5 Convertir extracción raster a formato largo
to_long_climate <- function(
values,
dates,
varname,
coords
) {
values <- as.data.frame(
values,
check.names = FALSE
)
if (nrow(values) != length(dates)) {
stop(
"El número de fechas no coincide con las capas raster."
)
}
colnames(values) <- paste0(
"V",
seq_len(ncol(values))
)
dt <- data.table::as.data.table(
values
)
dt[, DATE := as.Date(dates)]
long <- data.table::melt(
dt,
id.vars = "DATE",
variable.name = "coord_id",
value.name = varname
)
long[
,
coord_id := as.integer(
sub("^V", "", coord_id)
)
]
coords2 <- normalize_coords(
coords
)
coords2$coord_id <- seq_len(
nrow(coords2)
)
long <- merge(
long,
coords2,
by = "coord_id",
all.x = TRUE,
sort = FALSE
)
long[
order(DATE, coord_id)
]
}11 10. Función general de extracción climática
Esta función reemplaza las rutas absolutas por
here::here() y utiliza file.path() para que el
análisis funcione en cualquier computador.
La interpolación espacial continúa siendo bilineal.
extract_climate <- function(
coords,
country_name,
country_code,
exposure,
temp_types = c(
"max",
"mean",
"min"
)
) {
folder_name <- switch(
exposure,
TEMP = "Temperature",
WBGT = "WBGT",
UTCI = "UTCI",
HEAT_IND = "Heat_Index"
)
folder <- here::here(
"Clean",
folder_name
)
if (!dir.exists(folder)) {
stop(
"No existe la carpeta: ",
folder
)
}
files <- list.files(
folder,
full.names = FALSE
)
country_variants <- unique(
c(
country_name,
gsub(
" ",
"_",
country_name
)
)
)
country_match <- Reduce(
`|`,
lapply(
country_variants,
function(x) {
grepl(
x,
files,
fixed = TRUE
)
}
)
)
files <- files[
country_match
]
if (length(files) == 0) {
stop(
"No se encontraron archivos ",
exposure,
" para ",
country_name
)
}
final_data <- NULL
for (tt in temp_types) {
tag <- paste0(
"_",
tolower(tt),
"_"
)
tt_files <- files[
grepl(
tag,
tolower(files),
fixed = TRUE
)
]
if (length(tt_files) == 0) {
next
}
all_data <- NULL
for (f in tt_files) {
full_path <- file.path(
folder,
f
)
r <- tryCatch(
{
if (exposure == "HEAT_IND") {
tryCatch(
raster::brick(
full_path,
varname = "heatx",
level = 1,
stopIfNotEqualSpaced = FALSE
),
error = function(e) {
raster::brick(
full_path,
stopIfNotEqualSpaced = FALSE
)
}
)
} else {
raster::brick(
full_path,
stopIfNotEqualSpaced = FALSE
)
}
},
error = function(e) {
stop(
"Error abriendo ",
full_path,
": ",
conditionMessage(e)
)
}
)
coords_use <- clamp_uruguay(
r,
coords,
country_code
)
raw <- raster::extract(
x = r,
y = as.matrix(
coords_use[, c("lon", "lat")]
),
method = "bilinear",
df = TRUE
)
raw <- raw[, -1, drop = FALSE]
values <- as.data.frame(
t(raw)
)
values <- convert_to_celsius_if_needed(
values,
full_path,
exposure
)
dates <- dates_from_raster(
r,
full_path
)
varname <- paste0(
toupper(tt),
"_",
exposure
)
long <- to_long_climate(
values,
dates,
varname,
coords
)
all_data <- dplyr::bind_rows(
all_data,
long
)
}
if (
is.null(all_data) ||
nrow(all_data) == 0
) {
next
}
keys <- c(
"DATE",
"coord_id",
"lon",
"lat"
)
if (is.null(final_data)) {
final_data <- all_data
} else {
final_data <- dplyr::full_join(
final_data,
all_data,
by = keys
)
}
}
if (
is.null(final_data) ||
nrow(final_data) == 0
) {
stop(
"La extracción quedó vacía para ",
exposure,
" / ",
country_name
)
}
final_data %>%
dplyr::arrange(
DATE,
coord_id
)
}12 11. Preparación de la cohorte por país
prepare_country_cohort <- function(
data,
country_name,
country_code
) {
x <- data %>%
dplyr::mutate(
Pais = ifelse(
Pais == "Rep.Dom",
"Dominican Republic",
as.character(Pais)
),
LONGITUDE = ifelse(
Pais == "Guatemala",
-91.5160418,
LONGITUDE
),
LATITUDE = ifelse(
Pais == "Guatemala",
14.8438898,
LATITUDE
),
LONGITUDE = ifelse(
Institucion == "BolivianoJaponesBO",
-64.89618,
LONGITUDE
),
LATITUDE = ifelse(
Institucion == "BolivianoJaponesBO",
-14.82157,
LATITUDE
),
lon = LONGITUDE,
lat = LATITUDE
) %>%
dplyr::filter(
Pais == country_name,
Institucion != "CHPR_UY",
Anio %in% ANIOS
) %>%
dplyr::mutate(
date = as.Date(
VAR_0284,
format = "%Y-%m-%d"
),
year = year(date),
month = month(date),
dow = factor(
weekdays(date)
)
) %>%
dplyr::filter(
!is.na(date)
)
x$PTID <- factor(
paste0(
country_code,
"_ID_",
seq_len(nrow(x))
)
)
x$PTB <- factor(
x$Pretermino,
levels = c(0, 1)
)
x$ga <- as.numeric(
as.character(
x$VAR_0315
)
)
x$SB <- as.numeric(
as.character(
x$sb
)
)
x$SB[
x$ga < 28 |
is.na(x$ga)
] <- 0
x$SB <- factor(
x$SB,
levels = c(0, 1)
)
x$apgar5 <- as.numeric(
as.character(
x$apgar5
)
)
x$APGAR <- ifelse(
x$apgar5 < 7,
1,
0
)
x$APGAR <- factor(
x$APGAR,
levels = c(0, 1)
)
x$ANY_CS <- ifelse(
x$VAR_0287 == "B",
1,
0
)
x$ANY_CS <- factor(
x$ANY_CS,
levels = c(0, 1)
)
x$M_AGE <- factor(
x$M_AGE,
levels = c(
"<25",
"25-35",
">35"
)
)
x$I_SEX <- factor(
x$sexo,
levels = c(
"Female",
"Male"
)
)
x
}13 12. Construcción del diseño case-crossover
Se conservan como casos iniciales los registros con PTB, muerte fetal o Apgar bajo, tal como en el flujo de referencia.
Cada caso se compara con días control del mismo año, mismo mes y mismo día de la semana.
make_case_crossover <- function(x) {
cc <- x %>%
dplyr::select(
date,
year,
month,
dow,
PTID,
PTB,
SB,
APGAR,
M_AGE,
I_SEX,
ANY_CS,
LONGITUDE,
LATITUDE,
lon,
lat
) %>%
dplyr::filter(
PTB == 1 |
SB == 1 |
APGAR == 1
) %>%
dplyr::rename(
DATE = date
)
cc$APGAR[
cc$SB == 1
] <- 0
cc <- cc %>%
dplyr::mutate(
PTB = as.integer(
as.character(PTB)
),
SB = as.integer(
as.character(SB)
),
APGAR = as.integer(
as.character(APGAR)
),
ANY_CS = as.integer(
as.character(ANY_CS)
)
) %>%
dplyr::arrange(DATE)
all_dates <- tibble(
DATE = seq(
min(cc$DATE),
max(cc$DATE),
by = "days"
)
) %>%
dplyr::mutate(
year = year(DATE),
month = month(DATE),
dow = factor(
weekdays(DATE)
)
)
all_dates <- data.table::setDT(
all_dates
)
cc_dt <- data.table::setDT(
cc
)
out <- all_dates[
cc_dt[, .(
PTID,
PTB = 0L,
SB = 0L,
APGAR = 0L,
ANY_CS = 0L,
M_AGE,
I_SEX,
year,
month,
dow,
LATITUDE,
LONGITUDE,
lat,
lon
)],
on = .(
year,
month,
dow
),
allow.cartesian = TRUE
][
cc_dt,
on = .(
PTID,
DATE = DATE
),
`:=`(
PTB = i.PTB,
SB = i.SB,
APGAR = i.APGAR,
ANY_CS = i.ANY_CS
)
][]
as.data.frame(
out
) %>%
dplyr::arrange(
PTID,
DATE
)
}14 13. Asignación de exposiciones y creación de lags
deduplicate_climate <- function(x) {
x %>%
dplyr::mutate(
lon_r = round(lon, 5),
lat_r = round(lat, 5)
) %>%
dplyr::group_by(
DATE,
lon_r,
lat_r
) %>%
dplyr::summarise(
dplyr::across(
everything(),
\(z) {
if (is.numeric(z)) {
mean(z, na.rm = TRUE)
} else {
dplyr::first(z)
}
}
),
.groups = "drop"
) %>%
dplyr::select(
-dplyr::any_of(
c(
"lon",
"lat",
"coord_id"
)
)
)
}
join_lag <- function(
data,
climate,
pattern,
offset
) {
before <- names(data)
out <- data %>%
dplyr::mutate(
date_minus_offset =
DATE + 1 - days(offset)
) %>%
dplyr::left_join(
climate,
by = c(
"date_minus_offset" = "DATE",
"lon_r",
"lat_r"
)
)
new_cols <- setdiff(
names(out),
c(
before,
"date_minus_offset"
)
)
selected <- new_cols[
grepl(
pattern,
new_cols
)
]
names(out)[
match(
selected,
names(out)
)
] <- paste0(
selected,
"_",
offset - 1
)
out
}assign_exposures <- function(
cc,
country_name,
country_code
) {
coords <- cc %>%
dplyr::distinct(
lon,
lat
) %>%
as.data.frame()
temp <- extract_climate(
coords,
country_name,
country_code,
"TEMP"
)
wbgt <- extract_climate(
coords,
country_name,
country_code,
"WBGT"
)
utci <- extract_climate(
coords,
country_name,
country_code,
"UTCI"
)
heat <- extract_climate(
coords,
country_name,
country_code,
"HEAT_IND"
)
temp <- deduplicate_climate(temp)
wbgt <- deduplicate_climate(wbgt)
utci <- deduplicate_climate(utci)
heat <- deduplicate_climate(heat)
temp <- temp %>%
dplyr::select(
DATE,
lon_r,
lat_r,
dplyr::ends_with("_TEMP")
)
wbgt <- wbgt %>%
dplyr::select(
DATE,
lon_r,
lat_r,
dplyr::ends_with("_WBGT")
)
utci <- utci %>%
dplyr::select(
DATE,
lon_r,
lat_r,
dplyr::ends_with("_UTCI")
)
heat <- heat %>%
dplyr::select(
DATE,
lon_r,
lat_r,
dplyr::ends_with("_HEAT_IND")
)
out <- cc %>%
dplyr::mutate(
lon_r = round(lon, 5),
lat_r = round(lat, 5)
)
for (offset in 1:7) {
out <- join_lag(
out,
temp,
"_TEMP$",
offset
)
}
for (offset in 1:7) {
out <- join_lag(
out,
wbgt,
"_WBGT$",
offset
)
}
for (offset in 1:7) {
out <- join_lag(
out,
utci,
"_UTCI$",
offset
)
}
for (offset in 1:7) {
out <- join_lag(
out,
heat,
"_HEAT_IND$",
offset
)
}
out %>%
dplyr::arrange(DATE) %>%
dplyr::select(
-date_minus_offset
)
}15 14. Rango térmico diurno y desenlace compuesto
add_analysis_variables <- function(cc) {
cc <- cc %>%
dplyr::mutate(
COMPOSITE = as.integer(
rowSums(
dplyr::across(
c(
PTB,
SB,
APGAR,
ANY_CS
)
),
na.rm = TRUE
) > 0
)
)
exposures <- c(
"TEMP",
"WBGT",
"UTCI",
"HEAT_IND"
)
for (lag in 0:6) {
for (exposure in exposures) {
max_name <- paste0(
"MAX_",
exposure,
"_",
lag
)
min_name <- paste0(
"MIN_",
exposure,
"_",
lag
)
diurnal_name <- paste0(
"DIURNAL_",
exposure,
"_",
lag
)
if (
all(
c(
max_name,
min_name
) %in% names(cc)
)
) {
cc[[diurnal_name]] <-
cc[[max_name]] -
cc[[min_name]]
}
}
}
cc
}16 15. Procesar todos los países
Esta es la etapa computacional más pesada porque lee archivos NetCDF y realiza interpolación bilineal.
cc_by_country <- list()
for (code in names(PAISES)) {
country <- unname(
PAISES[code]
)
cat(
"\nProcesando:",
country,
"(", code, ")\n"
)
cohort <- prepare_country_cohort(
df_raw,
country,
code
)
if (nrow(cohort) == 0) {
warning(
"No hay registros para ",
country
)
next
}
cc <- make_case_crossover(
cohort
)
cc <- assign_exposures(
cc,
country,
code
)
cc <- add_analysis_variables(
cc
)
cc$Pais <- country
cc$country_code <- code
cc_by_country[[code]] <- cc
cat(
" Estratos:",
dplyr::n_distinct(cc$PTID),
"\n"
)
cat(
" Filas:",
nrow(cc),
"\n"
)
}##
## Procesando: Argentina ( ARG )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 9493 of `x` matches multiple rows in `y`.
## ℹ Row 9493 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 2003
## Filas: 8806
##
## Procesando: Bolivia ( BOL )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 9493 of `x` matches multiple rows in `y`.
## ℹ Row 9493 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 988
## Filas: 4350
##
## Procesando: Guatemala ( GTM )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4747 of `x` matches multiple rows in `y`.
## ℹ Row 4747 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4747 of `x` matches multiple rows in `y`.
## ℹ Row 4747 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 1271
## Filas: 5559
##
## Procesando: Honduras ( HND )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 14239 of `x` matches multiple rows in `y`.
## ℹ Row 14239 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 14239 of `x` matches multiple rows in `y`.
## ℹ Row 14239 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 964
## Filas: 4222
##
## Procesando: Dominican Republic ( DOM )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 9493 of `x` matches multiple rows in `y`.
## ℹ Row 9493 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 9493 of `x` matches multiple rows in `y`.
## ℹ Row 9493 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 2208
## Filas: 9711
##
## Procesando: Uruguay ( URY )## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4747 of `x` matches multiple rows in `y`.
## ℹ Row 4747 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Warning in dplyr::full_join(final_data, all_data, by = keys): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 4747 of `x` matches multiple rows in `y`.
## ℹ Row 4747 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Estratos: 897
## Filas: 394617 16. Control de calidad de la base analítica
quality_check_country <- function(cc, code) {
required_outcomes <- c(
"PTID",
"PTB",
"SB",
"APGAR",
"ANY_CS",
"COMPOSITE"
)
exposure_names <- unlist(
lapply(
c(
"TEMP",
"UTCI",
"WBGT",
"HEAT_IND"
),
function(e) {
unlist(
lapply(
c(
"MIN",
"MEAN",
"MAX",
"DIURNAL"
),
function(s) {
paste0(
s,
"_",
e,
"_",
0:6
)
}
)
)
}
)
)
missing_cols <- setdiff(
c(
required_outcomes,
exposure_names
),
names(cc)
)
estratos <- cc %>%
dplyr::count(
PTID,
name = "N"
)
tibble(
country_code = code,
n_rows = nrow(cc),
n_strata = dplyr::n_distinct(cc$PTID),
min_rows_per_stratum = min(estratos$N),
max_rows_per_stratum = max(estratos$N),
missing_required_columns = length(missing_cols),
SB_APGAR_overlap = sum(
cc$SB == 1 &
cc$APGAR == 1,
na.rm = TRUE
)
)
}quality_summary <- purrr::imap_dfr(
cc_by_country,
quality_check_country
)
quality_summary18 17. Descriptivos de desenlaces de la base analítica
outcome_summary <- purrr::imap_dfr(
cc_by_country,
function(cc, code) {
tibble(
country_code = code,
Pais = unique(cc$Pais),
Estratos = dplyr::n_distinct(cc$PTID),
PTB = sum(cc$PTB == 1, na.rm = TRUE),
SB = sum(cc$SB == 1, na.rm = TRUE),
APGAR = sum(cc$APGAR == 1, na.rm = TRUE),
ANY_CS = sum(cc$ANY_CS == 1, na.rm = TRUE),
COMPOSITE = sum(cc$COMPOSITE == 1, na.rm = TRUE)
)
}
)
DT::datatable(
outcome_summary,
rownames = FALSE,
options = list(scrollX = TRUE)
)19 18. Función central DLNM + regresión logística condicional
run_dlnm_model <- function(
model_data,
outcome,
outcome_label,
stat = STAT_MAIN,
exposure = EXPOSURE_MAIN,
lag_max = LAG_MAX,
fun = "ns",
degree = 2,
argvar_knots = ARGVAR_KNOTS,
arglag_nk = ARGLAG_NK,
center_percentile = CENTER_PERCENTILE,
comparator_percentile = COMPARATOR_PERCENTILE
) {
data <- model_data
data$OUTCOME <- data[[outcome]]
case_ids <- unique(
data$PTID[
data$OUTCOME == 1
]
)
data <- data %>%
dplyr::filter(
PTID %in% case_ids
)
exposure_vars <- paste0(
stat,
"_",
exposure,
"_",
0:lag_max
)
missing_exposures <- setdiff(
exposure_vars,
names(data)
)
if (length(missing_exposures) > 0) {
stop(
"Faltan variables: ",
paste(
missing_exposures,
collapse = ", "
)
)
}
Tavs <- as.matrix(
data[, exposure_vars]
)
valid_rows <- complete.cases(
Tavs,
data$OUTCOME,
data$PTID
)
data <- data[
valid_rows,
,
drop = FALSE
]
Tavs <- Tavs[
valid_rows,
,
drop = FALSE
]
if (nrow(data) == 0) {
stop("No quedan datos completos para el modelo.")
}
if (fun == "poly") {
if (lag_max < 3) {
cb <- dlnm::crossbasis(
Tavs,
lag = c(
0,
lag_max
),
argvar = list(
fun = "poly",
degree = degree
),
arglag = list(
fun = "poly"
)
)
} else {
cb <- dlnm::crossbasis(
Tavs,
lag = c(
0,
lag_max
),
argvar = list(
fun = "poly",
degree = degree
),
arglag = list(
fun = "poly",
degree = degree,
intercept = FALSE
)
)
}
} else {
exposure_knots <- quantile(
Tavs,
probs = argvar_knots / 100,
na.rm = TRUE
)
if (lag_max < 3) {
cb <- dlnm::crossbasis(
Tavs,
lag = c(
0,
lag_max
),
argvar = list(
fun = "ns",
knots = exposure_knots
),
arglag = list(
fun = "ns"
)
)
} else {
cb <- dlnm::crossbasis(
Tavs,
lag = c(
0,
lag_max
),
argvar = list(
fun = "ns",
knots = exposure_knots
),
arglag = list(
fun = "ns",
knots = dlnm::logknots(
lag_max,
arglag_nk
)
)
)
}
}
fit <- survival::clogit(
OUTCOME ~ cb + strata(PTID),
data = data
)
exposure0 <- data[[paste0(
stat,
"_",
exposure,
"_0"
)]]
percentiles <- quantile(
exposure0,
probs = c(
.01,
.05,
.50,
.75,
.95,
.99
),
na.rm = TRUE
)
center <- unname(
quantile(
exposure0,
center_percentile,
na.rm = TRUE
)
)
comparator <- unname(
quantile(
exposure0,
comparator_percentile,
na.rm = TRUE
)
)
center <- round(
center,
1
)
pred <- dlnm::crosspred(
cb,
fit,
cen = center,
by = 0.1
)
curve <- tibble(
Exposure = as.numeric(
pred$predvar
),
OR = as.numeric(
pred$allRRfit
),
Low = as.numeric(
pred$allRRlow
),
High = as.numeric(
pred$allRRhigh
)
)
comparator_est <- tibble(
OR = approx(
curve$Exposure,
curve$OR,
xout = comparator,
ties = "ordered"
)$y,
Low = approx(
curve$Exposure,
curve$Low,
xout = comparator,
ties = "ordered"
)$y,
High = approx(
curve$Exposure,
curve$High,
xout = comparator,
ties = "ordered"
)$y
)
summary <- tibble(
OUTCOME = outcome_label,
N = nrow(data),
N_STRATA = dplyr::n_distinct(
data$PTID
),
Exposure = paste0(
stat,
"_",
exposure
),
LAG = lag_max,
FUN = fun,
DEGREE = ifelse(
fun == "poly",
degree,
NA
),
ARGVAR_KNOTS = ifelse(
fun == "ns",
paste(
argvar_knots,
collapse = "_"
),
NA_character_
),
ARGLAG_LAGNK = ifelse(
fun == "ns",
arglag_nk,
NA
),
AIC = AIC(fit),
CENTERING_PERCENT = paste0(
center_percentile * 100,
"th"
),
CENTER_EXPOSURE = center,
HIGH_PERCENT = paste0(
comparator_percentile * 100,
"th"
),
HIGH_EXPOSURE = comparator,
OR = comparator_est$OR,
CILow = comparator_est$Low,
CIHigh = comparator_est$High,
Significant = (
comparator_est$Low > 1 |
comparator_est$High < 1
)
)
list(
data = data,
Tavs = Tavs,
crossbasis = cb,
model = fit,
prediction = pred,
curve = curve,
summary = summary,
percentiles = percentiles,
center = center,
comparator = comparator
)
}20 19. Ejecutar los modelos principales
outcomes <- c(
PTB = "Preterm Birth",
SB = "Still Birth",
APGAR = "APGAR Score",
ANY_CS = "Any CS",
COMPOSITE = "Composite"
)
main_models <- list()
main_results <- list()
main_curves <- list()
for (code in names(cc_by_country)) {
cc <- cc_by_country[[code]]
country <- unique(
cc$Pais
)
main_models[[code]] <- list()
for (outcome in names(outcomes)) {
cat(
"\nModelo principal:",
country,
"-",
outcomes[[outcome]],
"\n"
)
fit <- tryCatch(
run_dlnm_model(
model_data = cc,
outcome = outcome,
outcome_label = outcomes[[outcome]]
),
error = function(e) {
warning(
country,
" / ",
outcome,
": ",
conditionMessage(e)
)
NULL
}
)
if (is.null(fit)) {
next
}
main_models[[code]][[outcome]] <- fit
main_results[[paste(
code,
outcome,
sep = "_"
)]] <- fit$summary %>%
dplyr::mutate(
country_code = code,
Pais = country
)
main_curves[[paste(
code,
outcome,
sep = "_"
)]] <- fit$curve %>%
dplyr::mutate(
country_code = code,
Pais = country,
OUTCOME = outcomes[[outcome]],
Center = fit$center
)
}
}##
## Modelo principal: Argentina - Preterm Birth
##
## Modelo principal: Argentina - Still Birth
##
## Modelo principal: Argentina - APGAR Score
##
## Modelo principal: Argentina - Any CS
##
## Modelo principal: Argentina - Composite
##
## Modelo principal: Bolivia - Preterm Birth
##
## Modelo principal: Bolivia - Still Birth
##
## Modelo principal: Bolivia - APGAR Score
##
## Modelo principal: Bolivia - Any CS
##
## Modelo principal: Bolivia - Composite
##
## Modelo principal: Guatemala - Preterm Birth
##
## Modelo principal: Guatemala - Still Birth
##
## Modelo principal: Guatemala - APGAR Score
##
## Modelo principal: Guatemala - Any CS
##
## Modelo principal: Guatemala - Composite
##
## Modelo principal: Honduras - Preterm Birth
##
## Modelo principal: Honduras - Still Birth
##
## Modelo principal: Honduras - APGAR Score
##
## Modelo principal: Honduras - Any CS
##
## Modelo principal: Honduras - Composite
##
## Modelo principal: Dominican Republic - Preterm Birth
##
## Modelo principal: Dominican Republic - Still Birth
##
## Modelo principal: Dominican Republic - APGAR Score
##
## Modelo principal: Dominican Republic - Any CS
##
## Modelo principal: Dominican Republic - Composite
##
## Modelo principal: Uruguay - Preterm Birth
##
## Modelo principal: Uruguay - Still Birth
##
## Modelo principal: Uruguay - APGAR Score
##
## Modelo principal: Uruguay - Any CS
##
## Modelo principal: Uruguay - Compositemain_results_df <- dplyr::bind_rows(
main_results
)
main_curves_df <- dplyr::bind_rows(
main_curves
)21 20. Tabla de resultados principales
main_results_table <- main_results_df %>%
dplyr::mutate(
OR = round(OR, 2),
CILow = round(CILow, 2),
CIHigh = round(CIHigh, 2),
AIC = round(AIC, 2),
`OR (IC95%)` = paste0(
OR,
" (",
CILow,
"–",
CIHigh,
")"
)
) %>%
dplyr::select(
Pais,
OUTCOME,
N,
N_STRATA,
Exposure,
LAG,
FUN,
ARGVAR_KNOTS,
CENTER_EXPOSURE,
HIGH_EXPOSURE,
`OR (IC95%)`,
Significant,
AIC
)
DT::datatable(
main_results_table,
rownames = FALSE,
filter = "top",
extensions = "Buttons",
options = list(
pageLength = 25,
scrollX = TRUE,
dom = "Bfrtip",
buttons = c(
"copy",
"csv",
"excel"
)
),
caption = "Resultados principales recalculados durante esta sesión"
)22 21. Curvas exposición–respuesta por desenlace y país
Cada panel corresponde a un modelo independiente para ese país.
for (outcome_name in unique(main_curves_df$OUTCOME)) {
cat(
"\n\n## ",
outcome_name,
"\n\n",
sep = ""
)
p <- main_curves_df %>%
dplyr::filter(
OUTCOME == outcome_name
) %>%
ggplot(
aes(
x = Exposure,
y = OR
)
) +
geom_ribbon(
aes(
ymin = Low,
ymax = High
),
alpha = 0.2
) +
geom_line(
linewidth = 0.8
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
facet_wrap(
~ Pais,
scales = "free_x"
) +
labs(
title = paste0(
outcome_name,
": asociación acumulada MAX UTCI, lag 0–6"
),
subtitle = "Referencia: percentil 75. IC95% sombreado.",
x = "MAX UTCI (°C)",
y = "Odds Ratio"
) +
theme_minimal()
print(p)
cat("\n\n")
}22.1 Preterm Birth
22.2 Still Birth
22.3 APGAR Score
22.4 Any CS
22.5 Composite
23 22. Curvas individuales por país
for (code in names(cc_by_country)) {
country <- unique(
cc_by_country[[code]]$Pais
)
cat(
"\n\n## ",
country,
"\n\n",
sep = ""
)
for (outcome_name in unique(main_curves_df$OUTCOME)) {
datos <- main_curves_df %>%
dplyr::filter(
country_code == code,
OUTCOME == outcome_name
)
if (nrow(datos) == 0) {
next
}
p <- ggplot(
datos,
aes(
x = Exposure,
y = OR
)
) +
geom_ribbon(
aes(
ymin = Low,
ymax = High
),
alpha = 0.2
) +
geom_line(
linewidth = 0.8
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
labs(
title = paste(
country,
"—",
outcome_name
),
subtitle = "MAX UTCI; efecto acumulado lag 0–6; referencia P75",
x = "MAX UTCI (°C)",
y = "Odds Ratio"
) +
theme_minimal()
print(p)
}
}23.1 Argentina
23.2 Bolivia
23.3 Guatemala
23.4 Honduras
23.5 Dominican Republic
23.6 Uruguay
24 23. Percentiles de exposición utilizados
percentile_table <- purrr::imap_dfr(
main_models,
function(country_models, code) {
purrr::imap_dfr(
country_models,
function(fit, outcome) {
tibble(
country_code = code,
Pais = unique(
fit$data$Pais
),
OUTCOME = outcomes[[outcome]],
Percentile = names(
fit$percentiles
),
Exposure = as.numeric(
fit$percentiles
)
)
}
)
}
)
DT::datatable(
percentile_table,
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
)
)25 24. Análisis de sensibilidad
Esta sección utiliza los mismos datos recién calculados.
No lee resultados históricos.
Se prueban:
- polinomios de grado 2, 3 y 4;
- diferentes ventanas de lag;
- alternativas en los knots de exposición;
- otras métricas climáticas.
run_sensitivity_set <- function(
cc,
country,
code,
outcome,
outcome_label
) {
ans <- list()
# A. Polinomios
for (deg in c(2, 3, 4)) {
fit <- tryCatch(
run_dlnm_model(
cc,
outcome,
outcome_label,
fun = "poly",
degree = deg
),
error = function(e) NULL
)
if (!is.null(fit)) {
ans[[paste0(
"poly_",
deg
)]] <- fit$summary %>%
dplyr::mutate(
Sensitivity = "Parameterization",
Scenario = paste0(
"Polynomial degree ",
deg
)
)
}
}
# B. Lags
for (lag in 0:5) {
fit <- tryCatch(
run_dlnm_model(
cc,
outcome,
outcome_label,
lag_max = lag
),
error = function(e) NULL
)
if (!is.null(fit)) {
ans[[paste0(
"lag_",
lag
)]] <- fit$summary %>%
dplyr::mutate(
Sensitivity = "Lag",
Scenario = paste0(
"Lag 0–",
lag
)
)
}
}
# C. Knots de exposición
knot_sets <- list(
`10_90` = c(10, 90),
`5_50_95` = c(5, 50, 95),
`10_50_90` = c(10, 50, 90),
`5_75_95` = c(5, 75, 95)
)
for (nm in names(knot_sets)) {
fit <- tryCatch(
run_dlnm_model(
cc,
outcome,
outcome_label,
argvar_knots = knot_sets[[nm]]
),
error = function(e) NULL
)
if (!is.null(fit)) {
ans[[paste0(
"knots_",
nm
)]] <- fit$summary %>%
dplyr::mutate(
Sensitivity = "Exposure knots",
Scenario = nm
)
}
}
# D. Otras exposiciones
climate_options <- tibble::tribble(
~stat, ~exposure,
"MAX", "TEMP",
"MAX", "WBGT",
"MAX", "HEAT_IND",
"MEAN", "TEMP",
"MEAN", "UTCI",
"MEAN", "WBGT",
"MEAN", "HEAT_IND"
)
for (i in seq_len(nrow(climate_options))) {
fit <- tryCatch(
run_dlnm_model(
cc,
outcome,
outcome_label,
stat = climate_options$stat[i],
exposure = climate_options$exposure[i]
),
error = function(e) NULL
)
if (!is.null(fit)) {
ans[[paste0(
climate_options$stat[i],
"_",
climate_options$exposure[i]
)]] <- fit$summary %>%
dplyr::mutate(
Sensitivity = "Climate metric",
Scenario = paste0(
climate_options$stat[i],
" ",
climate_options$exposure[i]
)
)
}
}
dplyr::bind_rows(ans) %>%
dplyr::mutate(
country_code = code,
Pais = country
)
}sensitivity_results <- tibble()
if (RUN_SENSITIVITY) {
sens_list <- list()
for (code in names(cc_by_country)) {
cc <- cc_by_country[[code]]
country <- unique(
cc$Pais
)
for (outcome in names(outcomes)) {
cat(
"Sensibilidad:",
country,
"-",
outcomes[[outcome]],
"\n"
)
sens_list[[paste(
code,
outcome,
sep = "_"
)]] <- run_sensitivity_set(
cc,
country,
code,
outcome,
outcomes[[outcome]]
)
}
}
sensitivity_results <- dplyr::bind_rows(
sens_list
)
}## Sensibilidad: Argentina - Preterm Birth
## Sensibilidad: Argentina - Still Birth
## Sensibilidad: Argentina - APGAR Score
## Sensibilidad: Argentina - Any CS
## Sensibilidad: Argentina - Composite
## Sensibilidad: Bolivia - Preterm Birth
## Sensibilidad: Bolivia - Still Birth
## Sensibilidad: Bolivia - APGAR Score
## Sensibilidad: Bolivia - Any CS
## Sensibilidad: Bolivia - Composite
## Sensibilidad: Guatemala - Preterm Birth
## Sensibilidad: Guatemala - Still Birth
## Sensibilidad: Guatemala - APGAR Score
## Sensibilidad: Guatemala - Any CS
## Sensibilidad: Guatemala - Composite## Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,
## : Loglik converged before variable 5 ; beta may be infinite.## Sensibilidad: Honduras - Preterm Birth
## Sensibilidad: Honduras - Still Birth
## Sensibilidad: Honduras - APGAR Score
## Sensibilidad: Honduras - Any CS
## Sensibilidad: Honduras - Composite
## Sensibilidad: Dominican Republic - Preterm Birth
## Sensibilidad: Dominican Republic - Still Birth
## Sensibilidad: Dominican Republic - APGAR Score
## Sensibilidad: Dominican Republic - Any CS## Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,
## : Loglik converged before variable 12,16 ; beta may be infinite.## Sensibilidad: Dominican Republic - Composite
## Sensibilidad: Uruguay - Preterm Birth
## Sensibilidad: Uruguay - Still Birth## Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,
## : Ran out of iterations and did not converge## Warning in coxexact.fit(X, Y, istrat, offset, init, control, weights = weights,
## : Ran out of iterations and did not converge## Sensibilidad: Uruguay - APGAR Score
## Sensibilidad: Uruguay - Any CS
## Sensibilidad: Uruguay - Composite26 25. Tabla de sensibilidad
if (
RUN_SENSITIVITY &&
nrow(sensitivity_results) > 0
) {
sensitivity_table <- sensitivity_results %>%
dplyr::mutate(
OR = round(OR, 2),
CILow = round(CILow, 2),
CIHigh = round(CIHigh, 2),
`OR (IC95%)` = paste0(
OR,
" (",
CILow,
"–",
CIHigh,
")"
)
) %>%
dplyr::select(
Pais,
OUTCOME,
Sensitivity,
Scenario,
Exposure,
LAG,
FUN,
ARGVAR_KNOTS,
`OR (IC95%)`,
AIC
)
DT::datatable(
sensitivity_table,
rownames = FALSE,
filter = "top",
extensions = "Buttons",
options = list(
pageLength = 25,
scrollX = TRUE,
dom = "Bfrtip",
buttons = c(
"copy",
"csv",
"excel"
)
)
)
}27 26. Forest plot de resultados principales
forest_data <- main_results_df %>%
dplyr::mutate(
label = paste(
Pais,
OUTCOME,
sep = " — "
)
)
ggplot(
forest_data,
aes(
x = OR,
y = reorder(
label,
OR
)
)
) +
geom_errorbarh(
aes(
xmin = CILow,
xmax = CIHigh
),
height = 0.2
) +
geom_point() +
geom_vline(
xintercept = 1,
linetype = 2
) +
labs(
title = "Resultados principales: P99 vs P75 de MAX UTCI",
x = "Odds Ratio (IC95%)",
y = NULL
) +
theme_minimal()## Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
## ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `height` was translated to `width`.
28 27. Metaanálisis LAC para todos los desenlaces
El segundo nivel del análisis se generaliza ahora a todos los desenlaces estimados en los seis países:
- Parto pretérmino (
PTB) - Muerte fetal (
SB) - Apgar < 7 (
APGAR) - Cesárea (
ANY_CS) - Desenlace compuesto (
COMPOSITE)
Para cada desenlace se producen dos síntesis:
- metaanálisis P99 vs P75 de MAX UTCI;
- curva exposición–respuesta pooled OR vs MAX UTCI.
El procedimiento es idéntico para todos los desenlaces. Cuando un país no tiene un modelo válido o la estimación no es finita, ese país se excluye únicamente de ese desenlace.
29 27.1 Preparar efectos país-específicos P99 vs P75
meta_all_input <- main_results_df %>%
dplyr::filter(
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh),
OR > 0,
CILow > 0,
CIHigh > 0
) %>%
dplyr::mutate(
logOR = log(OR),
SE_logOR =
(
log(CIHigh) -
log(CILow)
) /
(
2 *
stats::qnorm(0.975)
),
VAR_logOR =
SE_logOR^2
)
DT::datatable(
meta_all_input %>%
dplyr::select(
Pais,
OUTCOME,
OR,
CILow,
CIHigh,
logOR,
SE_logOR
),
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
),
caption = "Efectos país-específicos utilizados en los metaanálisis"
)30 27.2 Función de metaanálisis P99 vs P75
pool_univariate_outcome <- function(data_outcome) {
data_outcome <- data_outcome %>%
dplyr::filter(
is.finite(logOR),
is.finite(VAR_logOR),
VAR_logOR > 0
)
k <- nrow(data_outcome)
if (k < 2) {
return(
tibble::tibble(
Countries = k,
OR = NA_real_,
CILow = NA_real_,
CIHigh = NA_real_,
I2 = NA_real_,
Q = NA_real_,
Tau2 = NA_real_
)
)
}
fit <- tryCatch(
mixmeta::mixmeta(
logOR ~ 1,
S = VAR_logOR,
data = data_outcome,
method = "reml"
),
error = function(e) NULL
)
if (is.null(fit)) {
return(
tibble::tibble(
Countries = k,
OR = NA_real_,
CILow = NA_real_,
CIHigh = NA_real_,
I2 = NA_real_,
Q = NA_real_,
Tau2 = NA_real_
)
)
}
beta <- as.numeric(
stats::coef(fit)[1]
)
se <- sqrt(
as.numeric(
stats::vcov(fit)[1, 1]
)
)
w_fixed <- 1 /
data_outcome$VAR_logOR
beta_fixed <- sum(
w_fixed *
data_outcome$logOR
) /
sum(w_fixed)
Q <- sum(
w_fixed *
(
data_outcome$logOR -
beta_fixed
)^2
)
df_Q <- k - 1
I2 <- if (
is.finite(Q) &&
Q > 0
) {
max(
0,
(
Q -
df_Q
) /
Q
) * 100
} else {
0
}
tau2 <- tryCatch(
as.numeric(
fit$Psi[1, 1]
),
error = function(e) NA_real_
)
tibble::tibble(
Countries = k,
OR = exp(beta),
CILow = exp(
beta -
1.96 * se
),
CIHigh = exp(
beta +
1.96 * se
),
I2 = I2,
Q = Q,
Tau2 = tau2
)
}31 27.3 Ejecutar metaanálisis P99 vs P75 para todos los desenlaces
meta_all_summary <- meta_all_input %>%
dplyr::group_by(
OUTCOME
) %>%
tidyr::nest() %>%
dplyr::mutate(
pooled = purrr::map(
data,
pool_univariate_outcome
)
) %>%
dplyr::select(
-data
) %>%
tidyr::unnest(
pooled
) %>%
dplyr::ungroup()
meta_all_summary %>%
dplyr::mutate(
dplyr::across(
where(is.numeric),
~ round(.x, 3)
)
)32 27.4 Forest plots por desenlace
for (outcome_name in unique(meta_all_input$OUTCOME)) {
dat <- meta_all_input %>%
dplyr::filter(
OUTCOME == outcome_name
)
pooled <- meta_all_summary %>%
dplyr::filter(
OUTCOME == outcome_name
)
if (
nrow(dat) < 2 ||
nrow(pooled) == 0 ||
!is.finite(pooled$OR)
) {
next
}
forest_data <- dat %>%
dplyr::select(
Pais,
OR,
CILow,
CIHigh
) %>%
dplyr::bind_rows(
tibble::tibble(
Pais = "LAC pooled",
OR = pooled$OR,
CILow = pooled$CILow,
CIHigh = pooled$CIHigh
)
) %>%
dplyr::mutate(
Pais = factor(
Pais,
levels = rev(
c(
dat$Pais,
"LAC pooled"
)
)
)
)
cat(
"\n\n## ",
outcome_name,
"\n\n",
sep = ""
)
p <- ggplot(
forest_data,
aes(
x = OR,
y = Pais
)
) +
geom_errorbarh(
aes(
xmin = CILow,
xmax = CIHigh
),
height = 0.18
) +
geom_point(
aes(
size =
Pais == "LAC pooled"
)
) +
scale_size_manual(
values = c(
`FALSE` = 2.5,
`TRUE` = 4
),
guide = "none"
) +
geom_vline(
xintercept = 1,
linetype = 2
) +
labs(
title = paste0(
outcome_name,
": metaanálisis LAC"
),
subtitle = paste0(
"MAX UTCI P99 vs P75; ",
pooled$Countries,
" países; I² = ",
round(
pooled$I2,
1
),
"%"
),
x = "Odds Ratio (IC95%)",
y = NULL
) +
theme_minimal()
print(p)
}32.1 Preterm Birth
## `height` was translated to `width`.
32.2 Still Birth
## `height` was translated to `width`.
32.3 APGAR Score
## `height` was translated to `width`.
32.4 Any CS
## `height` was translated to `width`.
32.5 Composite
## `height` was translated to `width`.
33 28. Curvas exposición–respuesta metaanalizadas para todos los desenlaces
Para la curva pooled se utiliza una referencia térmica común LAC. Esta se define como el percentil 75 de MAX UTCI lag 0 de todos los datos analíticos combinados.
34 28.1 Distribución conjunta de MAX UTCI
lac_utci <- purrr::imap_dfr(
cc_by_country,
function(cc, code) {
tibble::tibble(
country_code = code,
Pais = unique(
cc$Pais
),
MAX_UTCI_0 =
cc$MAX_UTCI_0
)
}
) %>%
dplyr::filter(
is.finite(
MAX_UTCI_0
)
)
lac_percentiles <- stats::quantile(
lac_utci$MAX_UTCI_0,
probs = c(
.01,
.05,
.10,
.25,
.50,
.75,
.90,
.95,
.99
),
na.rm = TRUE
)
LAC_CENTER <- as.numeric(
lac_percentiles["75%"]
)
LAC_GRID_MIN <- floor(
as.numeric(
lac_percentiles["1%"]
)
)
LAC_GRID_MAX <- ceiling(
as.numeric(
lac_percentiles["99%"]
)
)
LAC_GRID <- seq(
LAC_GRID_MIN,
LAC_GRID_MAX,
by = 0.1
)
lac_percentiles## 1% 5% 10% 25% 50% 75% 90% 95%
## 8.589901 16.536791 20.636595 26.122746 31.088047 34.217313 36.314245 37.609166
## 99%
## 39.544032cat(
"\nReferencia térmica común LAC:",
round(
LAC_CENTER,
2
),
"°C\n"
)##
## Referencia térmica común LAC: 34.22 °C35 28.2 Función para recalcular una curva país-específica con referencia común
make_common_center_curve <- function(
fit_obj,
country,
code,
outcome_name
) {
cb <- fit_obj$crossbasis
model <- fit_obj$model
coef_model <- stats::coef(
model
)
coef_cb <- coef_model[
grepl(
"^cb",
names(coef_model)
)
]
cat(
"\n",
country,
" / ",
outcome_name,
": columnas basis = ",
ncol(cb),
"; coeficientes DLNM = ",
length(coef_cb),
"\n",
sep = ""
)
if (
ncol(cb) !=
length(coef_cb)
) {
warning(
country,
" / ",
outcome_name,
": crossbasis y coeficientes no coinciden."
)
return(NULL)
}
exposure_country <- as.numeric(
fit_obj$Tavs
)
p01 <- as.numeric(
stats::quantile(
exposure_country,
.01,
na.rm = TRUE
)
)
p99 <- as.numeric(
stats::quantile(
exposure_country,
.99,
na.rm = TRUE
)
)
pred <- tryCatch(
dlnm::crosspred(
cb,
model,
cen = LAC_CENTER,
at = LAC_GRID
),
error = function(e) {
warning(
country,
" / ",
outcome_name,
": ",
conditionMessage(e)
)
NULL
}
)
if (is.null(pred)) {
return(NULL)
}
tibble::tibble(
Temperature =
as.numeric(
pred$predvar
),
OR =
as.numeric(
pred$allRRfit
),
Low =
as.numeric(
pred$allRRlow
),
High =
as.numeric(
pred$allRRhigh
)
) %>%
dplyr::mutate(
country_code = code,
Pais = country,
OUTCOME = outcome_name,
Supported =
Temperature >= p01 &
Temperature <= p99
) %>%
dplyr::filter(
Supported,
is.finite(OR),
is.finite(Low),
is.finite(High),
OR > 0,
Low > 0,
High > 0
) %>%
dplyr::mutate(
logOR = log(OR),
SE_logOR =
(
log(High) -
log(Low)
) /
(
2 *
stats::qnorm(0.975)
),
VAR_logOR =
SE_logOR^2
)
}36 28.3 Recalcular las curvas de todos los países y desenlaces
country_curves_all <- list()
for (code in names(main_models)) {
country <- unique(
cc_by_country[[code]]$Pais
)
for (outcome in names(outcomes)) {
if (
is.null(
main_models[[code]][[outcome]]
)
) {
next
}
outcome_name <-
outcomes[[outcome]]
tmp <- make_common_center_curve(
main_models[[code]][[outcome]],
country,
code,
outcome_name
)
if (is.null(tmp)) {
next
}
country_curves_all[[paste(
code,
outcome,
sep = "_"
)]] <- tmp
}
}##
## Argentina / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Argentina / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Argentina / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Argentina / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Argentina / Composite: columnas basis = 16; coeficientes DLNM = 16
##
## Bolivia / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Bolivia / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Bolivia / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Bolivia / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Bolivia / Composite: columnas basis = 16; coeficientes DLNM = 16
##
## Guatemala / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Guatemala / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Guatemala / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Guatemala / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Guatemala / Composite: columnas basis = 16; coeficientes DLNM = 16
##
## Honduras / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Honduras / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Honduras / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Honduras / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Honduras / Composite: columnas basis = 16; coeficientes DLNM = 16
##
## Dominican Republic / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Dominican Republic / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Dominican Republic / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Dominican Republic / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Dominican Republic / Composite: columnas basis = 16; coeficientes DLNM = 16
##
## Uruguay / Preterm Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Uruguay / Still Birth: columnas basis = 16; coeficientes DLNM = 16
##
## Uruguay / APGAR Score: columnas basis = 16; coeficientes DLNM = 16
##
## Uruguay / Any CS: columnas basis = 16; coeficientes DLNM = 16
##
## Uruguay / Composite: columnas basis = 16; coeficientes DLNM = 16country_curves_all_df <-
dplyr::bind_rows(
country_curves_all
)
dplyr::glimpse(
country_curves_all_df
)## Rows: 5,852
## Columns: 11
## $ Temperature <dbl> 8.0, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9.0, 9.…
## $ OR <dbl> 1.468174, 1.454959, 1.441933, 1.429095, 1.416443, 1.40397…
## $ Low <dbl> 0.8276352, 0.8250015, 0.8223217, 0.8195958, 0.8168239, 0.…
## $ High <dbl> 2.604450, 2.565941, 2.528415, 2.491854, 2.456236, 2.42154…
## $ country_code <chr> "ARG", "ARG", "ARG", "ARG", "ARG", "ARG", "ARG", "ARG", "…
## $ Pais <chr> "Argentina", "Argentina", "Argentina", "Argentina", "Arge…
## $ OUTCOME <chr> "Preterm Birth", "Preterm Birth", "Preterm Birth", "Prete…
## $ Supported <lgl> TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRU…
## $ logOR <dbl> 0.3840194, 0.3749775, 0.3659846, 0.3570415, 0.3481491, 0.…
## $ SE_logOR <dbl> 0.2924555, 0.2894684, 0.2865401, 0.2836713, 0.2808628, 0.…
## $ VAR_logOR <dbl> 0.08553020, 0.08379194, 0.08210521, 0.08046940, 0.0788839…37 28.4 Función de metaanálisis punto a punto
Cada temperatura se combina únicamente si existe información válida de al menos tres países.
pool_temperature_point <- function(
data_temperature
) {
data_temperature <- data_temperature %>%
dplyr::filter(
is.finite(logOR),
is.finite(VAR_logOR),
VAR_logOR > 0
)
k <- nrow(
data_temperature
)
if (k < 3) {
return(
tibble::tibble(
OR = NA_real_,
Low = NA_real_,
High = NA_real_,
N_countries = k,
Tau2 = NA_real_
)
)
}
fit_meta <- tryCatch(
mixmeta::mixmeta(
logOR ~ 1,
S = VAR_logOR,
data = data_temperature,
method = "reml"
),
error = function(e) NULL
)
if (is.null(fit_meta)) {
return(
tibble::tibble(
OR = NA_real_,
Low = NA_real_,
High = NA_real_,
N_countries = k,
Tau2 = NA_real_
)
)
}
beta <- as.numeric(
stats::coef(
fit_meta
)[1]
)
se <- sqrt(
as.numeric(
stats::vcov(
fit_meta
)[1, 1]
)
)
tau2 <- tryCatch(
as.numeric(
fit_meta$Psi[1, 1]
),
error = function(e) NA_real_
)
tibble::tibble(
OR = exp(beta),
Low = exp(
beta -
1.96 * se
),
High = exp(
beta +
1.96 * se
),
N_countries = k,
Tau2 = tau2
)
}38 28.5 Ejecutar el metaanálisis dosis–respuesta para todos los desenlaces
pooled_curves_all <- country_curves_all_df %>%
dplyr::group_by(
OUTCOME,
Temperature
) %>%
tidyr::nest() %>%
dplyr::mutate(
meta =
purrr::map(
data,
pool_temperature_point
)
) %>%
dplyr::select(
-data
) %>%
tidyr::unnest(
meta
) %>%
dplyr::ungroup() %>%
dplyr::filter(
N_countries >= 3,
is.finite(OR),
is.finite(Low),
is.finite(High)
)
dplyr::glimpse(
pooled_curves_all
)## Rows: 1,150
## Columns: 7
## $ Temperature <dbl> 16.8, 16.9, 17.0, 17.1, 17.2, 17.3, 17.4, 17.5, 17.6, 17.7…
## $ OUTCOME <chr> "Preterm Birth", "Preterm Birth", "Preterm Birth", "Preter…
## $ OR <dbl> 1.106709, 1.107205, 1.107776, 1.108421, 1.109136, 1.109919…
## $ Low <dbl> 0.6813719, 0.6819140, 0.6825782, 0.6833592, 0.6842523, 0.6…
## $ High <dbl> 1.797559, 1.797738, 1.797843, 1.797877, 1.797849, 1.797762…
## $ N_countries <int> 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3…
## $ Tau2 <dbl> 0.07023516, 0.07047623, 0.07066018, 0.07078958, 0.07086699…39 28.6 Curvas pooled consolidadas en una sola figura
ggplot(
pooled_curves_all,
aes(
x = Temperature,
y = OR
)
) +
geom_ribbon(
aes(
ymin = Low,
ymax = High
),
alpha = 0.20
) +
geom_line(
linewidth = 0.9
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
geom_vline(
xintercept =
LAC_CENTER,
linetype = 3
) +
facet_wrap(
~ OUTCOME,
scales = "free_y"
) +
labs(
title =
"Curvas exposición–respuesta metaanalizadas LAC",
subtitle =
paste0(
"MAX UTCI, lag acumulado 0–6; referencia común P75 = ",
round(
LAC_CENTER,
1
),
" °C"
),
x = "MAX UTCI (°C)",
y = "Odds Ratio",
caption =
"Cada punto pooled requiere información válida de al menos 3 países."
) +
theme_minimal()
40 28.7 Curvas pooled individuales por desenlace
for (
outcome_name in
unique(
pooled_curves_all$OUTCOME
)
) {
dat <- pooled_curves_all %>%
dplyr::filter(
OUTCOME ==
outcome_name
)
if (nrow(dat) == 0) {
next
}
cat(
"\n\n## ",
outcome_name,
"\n\n",
sep = ""
)
p <- ggplot(
dat,
aes(
x = Temperature,
y = OR
)
) +
geom_ribbon(
aes(
ymin = Low,
ymax = High
),
alpha = 0.20
) +
geom_line(
linewidth = 1
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
geom_vline(
xintercept =
LAC_CENTER,
linetype = 3
) +
labs(
title = paste0(
outcome_name,
": curva pooled LAC"
),
subtitle =
"MAX UTCI, efecto acumulado lag 0–6",
x = "MAX UTCI (°C)",
y = "Odds Ratio"
) +
theme_minimal()
print(p)
}40.1 Preterm Birth
40.2 Still Birth
40.3 APGAR Score
40.4 Any CS
40.5 Composite
41 28.8 Curvas país-específicas y pooled por desenlace
ggplot() +
geom_line(
data =
country_curves_all_df,
aes(
x = Temperature,
y = OR,
group = Pais
),
alpha = 0.25,
linewidth = 0.45
) +
geom_ribbon(
data =
pooled_curves_all,
aes(
x = Temperature,
ymin = Low,
ymax = High
),
alpha = 0.18
) +
geom_line(
data =
pooled_curves_all,
aes(
x = Temperature,
y = OR
),
linewidth = 1
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
geom_vline(
xintercept =
LAC_CENTER,
linetype = 3
) +
facet_wrap(
~ OUTCOME,
scales = "free_y"
) +
labs(
title =
"Curvas país-específicas y metaanalizadas",
subtitle =
"Líneas finas: países; línea principal: pooled LAC",
x = "MAX UTCI (°C)",
y = "Odds Ratio"
) +
theme_minimal()
42 28.9 Número de países que contribuyen por desenlace y temperatura
ggplot(
pooled_curves_all,
aes(
x = Temperature,
y = N_countries
)
) +
geom_step() +
facet_wrap(
~ OUTCOME
) +
scale_y_continuous(
breaks = 3:6
) +
labs(
title =
"Número de países que contribuyen a cada curva pooled",
x = "MAX UTCI (°C)",
y = "Países"
) +
theme_minimal()
43 29. Leave-one-country-out para todos los desenlaces
leave_one_out_all_all <- purrr::map_dfr(
unique(
meta_all_input$OUTCOME
),
function(outcome_name) {
dat <- meta_all_input %>%
dplyr::filter(
OUTCOME ==
outcome_name
)
if (nrow(dat) < 3) {
return(NULL)
}
purrr::map_dfr(
dat$Pais,
function(country_removed) {
x <- dat %>%
dplyr::filter(
Pais !=
country_removed
)
if (nrow(x) < 2) {
return(NULL)
}
fit <- tryCatch(
mixmeta::mixmeta(
logOR ~ 1,
S = VAR_logOR,
data = x,
method = "reml"
),
error = function(e) NULL
)
if (is.null(fit)) {
return(NULL)
}
beta <- as.numeric(
stats::coef(
fit
)[1]
)
se <- sqrt(
as.numeric(
stats::vcov(
fit
)[1, 1]
)
)
tibble::tibble(
OUTCOME =
outcome_name,
Excluded =
country_removed,
OR =
exp(beta),
CILow =
exp(
beta -
1.96 * se
),
CIHigh =
exp(
beta +
1.96 * se
)
)
}
)
}
)
DT::datatable(
leave_one_out_all_all,
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
)
)44 30. Resumen final del metaanálisis por desenlace
meta_final_summary <- meta_all_summary %>%
dplyr::mutate(
Contrast =
"MAX UTCI P99 vs P75",
Method =
"Random-effects REML (mixmeta)",
Dose_response_reference =
LAC_CENTER
) %>%
dplyr::select(
OUTCOME,
Countries,
Contrast,
OR,
CILow,
CIHigh,
I2,
Q,
Tau2,
Method,
Dose_response_reference
)
meta_final_summary %>%
dplyr::mutate(
dplyr::across(
where(is.numeric),
~ round(.x, 3)
)
)45 31. Sensibilidad metaanalizada por contraste percentilar
El modelo principal compara el P99 frente al P75 de MAX UTCI. El manuscrito también presenta como análisis de sensibilidad el contraste P95 frente al P50.
Para ampliar la interpretación, este Rmd calcula además tres contrastes exploratorios:
- P99 vs P50
- P95 vs P75
- P90 vs P50
Es importante distinguirlos:
P99 vs P75 → análisis principal del manuscrito
P95 vs P50 → sensibilidad reportada en el manuscrito
P99 vs P50 → exploratorio adicional
P95 vs P75 → exploratorio adicional
P90 vs P50 → exploratorio adicionalCambiar el percentil de referencia o el percentil comparador
no requiere volver a estimar el clogit. Se utiliza el
mismo modelo DLNM ya ajustado para cada país y se recalculan las
predicciones mediante crosspred().
46 31.1 Definir los contrastes
percentile_contrasts <- tibble::tribble(
~Contrast, ~Center_p, ~Comparator_p, ~Evidence,
"P99 vs P75", 0.75, 0.99, "Manuscrito: principal",
"P95 vs P50", 0.50, 0.95, "Manuscrito: sensibilidad",
"P99 vs P50", 0.50, 0.99, "Exploratorio adicional",
"P95 vs P75", 0.75, 0.95, "Exploratorio adicional",
"P90 vs P50", 0.50, 0.90, "Exploratorio adicional"
)
percentile_contrasts47 31.2 Función para obtener un contraste desde un modelo país-específico
estimate_percentile_contrast <- function(
fit_obj,
center_p,
comparator_p,
contrast_label,
evidence_label,
country,
code,
outcome_name
) {
cb <- fit_obj$crossbasis
model <- fit_obj$model
exposure0 <- fit_obj$data[[paste0(
STAT_MAIN,
"_",
EXPOSURE_MAIN,
"_0"
)]]
exposure0 <- as.numeric(
exposure0
)
center_value <- as.numeric(
stats::quantile(
exposure0,
probs = center_p,
na.rm = TRUE
)
)
comparator_value <- as.numeric(
stats::quantile(
exposure0,
probs = comparator_p,
na.rm = TRUE
)
)
pred <- tryCatch(
dlnm::crosspred(
cb,
model,
cen = center_value,
at = comparator_value
),
error = function(e) {
warning(
country,
" / ",
outcome_name,
" / ",
contrast_label,
": ",
conditionMessage(e)
)
NULL
}
)
if (is.null(pred)) {
return(NULL)
}
OR <- as.numeric(
pred$allRRfit[1]
)
CILow <- as.numeric(
pred$allRRlow[1]
)
CIHigh <- as.numeric(
pred$allRRhigh[1]
)
if (
!all(
is.finite(
c(
OR,
CILow,
CIHigh
)
)
) ||
OR <= 0 ||
CILow <= 0 ||
CIHigh <= 0
) {
return(NULL)
}
SE_logOR <- (
log(CIHigh) -
log(CILow)
) /
(
2 *
stats::qnorm(0.975)
)
tibble::tibble(
country_code = code,
Pais = country,
OUTCOME = outcome_name,
Contrast = contrast_label,
Evidence = evidence_label,
Center_p = center_p,
Comparator_p = comparator_p,
Center_value = center_value,
Comparator_value = comparator_value,
OR = OR,
CILow = CILow,
CIHigh = CIHigh,
logOR = log(OR),
SE_logOR = SE_logOR,
VAR_logOR = SE_logOR^2
)
}48 31.3 Calcular todos los contrastes en todos los países y desenlaces
contrast_country_list <- list()
for (code in names(main_models)) {
country <- unique(
cc_by_country[[code]]$Pais
)
for (outcome in names(outcomes)) {
fit_obj <- main_models[[code]][[outcome]]
if (is.null(fit_obj)) {
next
}
outcome_name <- outcomes[[outcome]]
for (
i in seq_len(
nrow(
percentile_contrasts
)
)
) {
tmp <- estimate_percentile_contrast(
fit_obj = fit_obj,
center_p =
percentile_contrasts$Center_p[i],
comparator_p =
percentile_contrasts$Comparator_p[i],
contrast_label =
percentile_contrasts$Contrast[i],
evidence_label =
percentile_contrasts$Evidence[i],
country = country,
code = code,
outcome_name = outcome_name
)
if (is.null(tmp)) {
next
}
key <- paste(
code,
outcome,
percentile_contrasts$Contrast[i],
sep = "_"
)
contrast_country_list[[key]] <- tmp
}
}
}
contrast_country_results <- dplyr::bind_rows(
contrast_country_list
)
DT::datatable(
contrast_country_results,
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
),
caption =
"Contrastes percentilares calculados por país y desenlace"
)49 31.4 Metaanalizar cada contraste para cada desenlace
contrast_meta_results <- contrast_country_results %>%
dplyr::group_by(
OUTCOME,
Contrast,
Evidence
) %>%
tidyr::nest() %>%
dplyr::mutate(
pooled = purrr::map(
data,
pool_univariate_outcome
)
) %>%
dplyr::select(
-data
) %>%
tidyr::unnest(
pooled
) %>%
dplyr::ungroup() %>%
dplyr::mutate(
Contrast = factor(
Contrast,
levels = percentile_contrasts$Contrast
)
)
DT::datatable(
contrast_meta_results %>%
dplyr::mutate(
dplyr::across(
c(
OR,
CILow,
CIHigh,
I2,
Q,
Tau2
),
~ round(.x, 3)
)
),
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
),
caption =
"Metaanálisis de contrastes percentilares por desenlace"
)50 31.5 Gráfica consolidada de contrastes metaanalizados
ggplot(
contrast_meta_results %>%
dplyr::filter(
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh)
),
aes(
x = OR,
y = Contrast
)
) +
geom_errorbarh(
aes(
xmin = CILow,
xmax = CIHigh
),
height = 0.18
) +
geom_point(
size = 2.5
) +
geom_vline(
xintercept = 1,
linetype = 2
) +
facet_wrap(
~ OUTCOME,
scales = "free_x"
) +
labs(
title =
"Sensibilidad del metaanálisis según contraste térmico",
subtitle =
"MAX UTCI, efecto acumulado lag 0–6 días",
x = "Odds Ratio pooled (IC95%)",
y = NULL
) +
theme_minimal()## `height` was translated to `width`.
51 31.6 Forest plots de contrastes por desenlace
for (
outcome_name in
unique(
contrast_meta_results$OUTCOME
)
) {
dat <- contrast_meta_results %>%
dplyr::filter(
OUTCOME == outcome_name,
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh)
)
if (nrow(dat) == 0) {
next
}
cat(
"\n\n## ",
outcome_name,
"\n\n",
sep = ""
)
p <- ggplot(
dat,
aes(
x = OR,
y = Contrast
)
) +
geom_errorbarh(
aes(
xmin = CILow,
xmax = CIHigh
),
height = 0.18
) +
geom_point(
size = 2.8
) +
geom_vline(
xintercept = 1,
linetype = 2
) +
labs(
title = paste0(
outcome_name,
": contrastes percentilares"
),
subtitle =
"MAX UTCI, lag acumulado 0–6 días",
x = "Odds Ratio pooled (IC95%)",
y = NULL
) +
theme_minimal()
print(p)
}51.1 Preterm Birth
## `height` was translated to `width`.
51.2 Still Birth
## `height` was translated to `width`.
51.3 APGAR Score
## `height` was translated to `width`.
51.4 Any CS
## `height` was translated to `width`.
51.5 Composite
## `height` was translated to `width`.
52 32. Sensibilidad metaanalizada por ventana de lag
El manuscrito compara la exposición durante:
lag 0
lag 0–1
lag 0–2
lag 0–3
lag 0–4
lag 0–5
lag 0–6 ← modelo principalLos modelos de lag 0–5 ya fueron estimados en la sección de sensibilidad. El modelo principal aporta lag 0–6.
53 32.1 Reunir resultados país-específicos por lag
lag_country_sensitivity <- sensitivity_results %>%
dplyr::filter(
Sensitivity == "Lag"
) %>%
dplyr::mutate(
LAG = as.integer(
LAG
)
)
lag_country_main <- main_results_df %>%
dplyr::mutate(
Sensitivity = "Lag",
Scenario = paste0(
"Lag 0–",
LAG
)
)
lag_country_results <- dplyr::bind_rows(
lag_country_sensitivity,
lag_country_main
) %>%
dplyr::filter(
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh),
OR > 0,
CILow > 0,
CIHigh > 0
) %>%
dplyr::mutate(
logOR = log(OR),
SE_logOR =
(
log(CIHigh) -
log(CILow)
) /
(
2 *
stats::qnorm(0.975)
),
VAR_logOR =
SE_logOR^2,
Lag_label = dplyr::case_when(
LAG == 0 ~ "Lag 0",
TRUE ~ paste0(
"Lag 0–",
LAG
)
)
)
DT::datatable(
lag_country_results %>%
dplyr::select(
Pais,
OUTCOME,
LAG,
Lag_label,
OR,
CILow,
CIHigh
),
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
),
caption =
"Efectos país-específicos según ventana de lag"
)54 32.2 Metaanalizar cada ventana de lag
lag_meta_results <- lag_country_results %>%
dplyr::group_by(
OUTCOME,
LAG,
Lag_label
) %>%
tidyr::nest() %>%
dplyr::mutate(
pooled = purrr::map(
data,
pool_univariate_outcome
)
) %>%
dplyr::select(
-data
) %>%
tidyr::unnest(
pooled
) %>%
dplyr::ungroup() %>%
dplyr::arrange(
OUTCOME,
LAG
)
DT::datatable(
lag_meta_results %>%
dplyr::mutate(
dplyr::across(
c(
OR,
CILow,
CIHigh,
I2,
Q,
Tau2
),
~ round(.x, 3)
)
),
rownames = FALSE,
filter = "top",
options = list(
pageLength = 25,
scrollX = TRUE
),
caption =
"Metaanálisis por ventana de lag y desenlace"
)55 32.3 Gráfica consolidada de lags metaanalizados
ggplot(
lag_meta_results %>%
dplyr::filter(
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh)
),
aes(
x = LAG,
y = OR
)
) +
geom_errorbar(
aes(
ymin = CILow,
ymax = CIHigh
),
width = 0.12
) +
geom_point(
size = 2.5
) +
geom_line(
aes(
group = 1
)
) +
geom_hline(
yintercept = 1,
linetype = 2
) +
scale_x_continuous(
breaks = 0:6,
labels = c(
"0",
"0–1",
"0–2",
"0–3",
"0–4",
"0–5",
"0–6"
)
) +
facet_wrap(
~ OUTCOME,
scales = "free_y"
) +
labs(
title =
"Sensibilidad del metaanálisis según ventana de exposición",
subtitle =
"Contraste P99 vs P75 de MAX UTCI",
x = "Ventana de lag (días)",
y = "Odds Ratio pooled (IC95%)"
) +
theme_minimal()
56 32.4 Forest plots de lag por desenlace
for (
outcome_name in
unique(
lag_meta_results$OUTCOME
)
) {
dat <- lag_meta_results %>%
dplyr::filter(
OUTCOME == outcome_name,
is.finite(OR),
is.finite(CILow),
is.finite(CIHigh)
) %>%
dplyr::mutate(
Lag_label = factor(
Lag_label,
levels = rev(
c(
"Lag 0",
"Lag 0–1",
"Lag 0–2",
"Lag 0–3",
"Lag 0–4",
"Lag 0–5",
"Lag 0–6"
)
)
)
)
if (nrow(dat) == 0) {
next
}
cat(
"\n\n## ",
outcome_name,
"\n\n",
sep = ""
)
p <- ggplot(
dat,
aes(
x = OR,
y = Lag_label
)
) +
geom_errorbarh(
aes(
xmin = CILow,
xmax = CIHigh
),
height = 0.18
) +
geom_point(
size = 2.8
) +
geom_vline(
xintercept = 1,
linetype = 2
) +
labs(
title = paste0(
outcome_name,
": sensibilidad por lag"
),
subtitle =
"MAX UTCI, contraste P99 vs P75",
x = "Odds Ratio pooled (IC95%)",
y = NULL
) +
theme_minimal()
print(p)
}56.1 APGAR Score
## `height` was translated to `width`.
56.2 Any CS
## `height` was translated to `width`.
56.3 Composite
## `height` was translated to `width`.
56.4 Preterm Birth
## `height` was translated to `width`.
56.5 Still Birth
## `height` was translated to `width`.
57 33. Exportación opcional de resultados recién calculados
Los archivos de esta sección son nuevos y proceden exclusivamente de los objetos calculados durante el Knit.
La exportación se hace de forma segura: si un objeto opcional no fue generado porque un análisis no pudo estimarse, el Knit continúa y simplemente informa que ese archivo no se exportó.
output_dir <- here::here(
"Rmd_Resultados_Recalculados"
)
dir.create(
output_dir,
showWarnings = FALSE
)
export_if_exists <- function(
object_name,
filename
) {
if (
exists(
object_name,
envir = knitr::knit_global()
)
) {
obj <- get(
object_name,
envir = knitr::knit_global()
)
utils::write.csv(
obj,
file.path(
output_dir,
filename
),
row.names = FALSE
)
cat(
"Exportado:",
filename,
"\n"
)
} else {
cat(
"No exportado:",
filename,
"— objeto",
object_name,
"no existe en esta sesión.\n"
)
}
invisible(NULL)
}
export_if_exists(
"main_results_df",
"resultados_principales_recalculados.csv"
)## Exportado: resultados_principales_recalculados.csvexport_if_exists(
"meta_all_input",
"meta_todos_desenlaces_paises.csv"
)## Exportado: meta_todos_desenlaces_paises.csvexport_if_exists(
"pooled_curves_all",
"curvas_meta_lac_todos_desenlaces_max_utci.csv"
)## Exportado: curvas_meta_lac_todos_desenlaces_max_utci.csvexport_if_exists(
"leave_one_out_all",
"meta_todos_desenlaces_leave_one_out.csv"
)## No exportado: meta_todos_desenlaces_leave_one_out.csv — objeto leave_one_out_all no existe en esta sesión.export_if_exists(
"contrast_country_results",
"contrastes_percentilares_por_pais.csv"
)## Exportado: contrastes_percentilares_por_pais.csvexport_if_exists(
"contrast_meta_results",
"meta_contrastes_percentilares.csv"
)## Exportado: meta_contrastes_percentilares.csvexport_if_exists(
"lag_country_results",
"lags_por_pais.csv"
)## Exportado: lags_por_pais.csvexport_if_exists(
"lag_meta_results",
"meta_sensibilidad_lags.csv"
)## Exportado: meta_sensibilidad_lags.csvexport_if_exists(
"percentile_table",
"percentiles_exposicion_recalculados.csv"
)## Exportado: percentiles_exposicion_recalculados.csvexport_if_exists(
"sensitivity_results",
"analisis_sensibilidad_recalculado.csv"
)## Exportado: analisis_sensibilidad_recalculado.csvcat(
"\nResultados escritos en:\n",
output_dir,
"\n"
)##
## Resultados escritos en:
## C:/Users/GERMAN/Dropbox/Clima (1)/Rmd_Resultados_Recalculados58 34. Interpretación de las curvas
La curva representa el efecto acumulado de la exposición durante lag 0–6.
El eje X corresponde a MAX UTCI en °C.
El eje Y corresponde al Odds Ratio respecto del percentil 75.
- OR = 1: misma odds que la exposición de referencia.
- OR > 1: mayor odds que la referencia.
- OR < 1: menor odds que la referencia.
- La banda corresponde al IC95%.
- El análisis principal resume P99 frente a P75.
La curva no debe interpretarse como un efecto lineal “por cada grado”, porque la exposición se modela mediante spline natural.
59 35. Trazabilidad de objetos calculados
tibble::tibble(
Objeto = c(
"df_raw",
"cc_by_country",
"main_models",
"main_results_df",
"main_curves_df",
"percentile_table",
"sensitivity_results",
"contrast_country_results",
"contrast_meta_results",
"lag_country_results",
"lag_meta_results"
),
Existe = c(
exists("df_raw"),
exists("cc_by_country"),
exists("main_models"),
exists("main_results_df"),
exists("main_curves_df"),
exists("percentile_table"),
exists("sensitivity_results"),
exists("contrast_country_results"),
exists("contrast_meta_results"),
exists("lag_country_results"),
exists("lag_meta_results")
)
)60 36. Información de sesión
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8 LC_CTYPE=Spanish_Spain.utf8
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] htmltools_0.5.8.1 DT_0.33 scales_1.4.0 ggplot2_4.0.3
## [5] mixmeta_1.2.2 dlnm_2.4.10 survival_3.8-3 ncdf4_1.24
## [9] raster_3.6-32 sp_2.2-0 data.table_1.17.2 lubridate_1.9.4
## [13] stringr_1.5.1 tibble_3.2.1 purrr_1.0.4 tidyr_1.3.1
## [17] dplyr_1.1.4 here_1.0.2
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.52 bslib_0.9.0 htmlwidgets_1.6.4
## [5] lattice_0.22-6 vctrs_0.6.5 tools_4.5.0 crosstalk_1.2.1
## [9] generics_0.1.4 proxy_0.4-27 pacman_0.5.1 pkgconfig_2.0.3
## [13] Matrix_1.7-3 KernSmooth_2.23-26 RColorBrewer_1.1-3 S7_0.2.2
## [17] lifecycle_1.0.4 compiler_4.5.0 farver_2.1.2 terra_1.8-60
## [21] codetools_0.2-20 class_7.3-23 sass_0.4.10 yaml_2.3.10
## [25] crayon_1.5.3 pillar_1.10.2 jquerylib_0.1.4 classInt_0.4-11
## [29] cachem_1.1.0 nlme_3.1-168 tidyselect_1.2.1 digest_0.6.37
## [33] stringi_1.8.7 sf_1.0-21 labeling_0.4.3 rprojroot_2.1.1
## [37] fastmap_1.2.0 grid_4.5.0 cli_3.6.5 magrittr_2.0.3
## [41] dichromat_2.0-0.1 e1071_1.7-16 withr_3.0.2 timechange_0.3.0
## [45] rmarkdown_2.29 tsModel_0.6-2 evaluate_1.0.3 knitr_1.50
## [49] mgcv_1.9-1 rlang_1.1.6 Rcpp_1.1.2 glue_1.8.0
## [53] DBI_1.2.3 rstudioapi_0.17.1 jsonlite_2.0.0 R6_2.6.1
## [57] units_0.8-7