# =============================================================================
# Validación con R survey: intervalos de confianza y pruebas con diseño complejo
# Dengue en el Perú 2026: ENDES 2025 y ENAHO 2025-2026
#
# Segunda implementación, independiente de los scripts de Python: construye los
# indicadores desde los microdatos, define el diseño sobre el marco completo y
# estima los dominios con subset() (no filtrando los datos antes del diseño).
#
# Uso:  Rscript r/validacion_survey.R [carpeta_endes] [carpeta_enaho] [carpeta_salida] [carga_dengue_depto.csv]
#       (por defecto: insumos/microdatos/endes, insumos/microdatos/enaho, salida_r)
# Requiere: R >= 4.1, survey >= 4.2, data.table
# =============================================================================

suppressPackageStartupMessages({ library(survey); library(data.table) })
invisible(Sys.setlocale("LC_CTYPE", "C.UTF-8"))   # nombres con tilde (Áncash, Junín…) iguales en todas las tablas
args <- commandArgs(trailingOnly = TRUE)
ENDES <- if (length(args) >= 1) args[1] else "insumos/microdatos/endes"
ENAHO <- if (length(args) >= 2) args[2] else "insumos/microdatos/enaho"
OUT   <- if (length(args) >= 3) args[3] else "salida_r"
CARGA <- if (length(args) >= 4) args[4] else "insumos/carga_dengue_depto.csv"
dir.create(OUT, showWarnings = FALSE, recursive = TRUE)

# Especificación principal ------------------------------------------------------
# - Diseño sobre el marco analítico completo (hogares entrevistados, mujeres
#   entrevistadas, personas del módulo 400); dominios con subset().
# - Estratos con una sola UPM: "adjust" (centrado en la media general; conservador).
# - IC de proporciones: método logit con t y gl del dominio (svyciprop).
# - Efecto de diseño: varianza frente a MAS con reemplazo (svymean, deff='replace').
# - Pruebas: Rao-Scott de segundo orden, estadístico F (svychisq).
LONELY <- Sys.getenv("SURVEY_LONELY", "adjust")   # "certainty" para comparar con samplics (single_psu = skip)
options(survey.lonely.psu = LONELY)
t0 <- Sys.time()
msg <- function(...) cat(sprintf("[%5.0f s] ", as.numeric(Sys.time() - t0, units = "secs")), ..., "\n", sep = "")

DEP <- c("Amazonas", "Áncash", "Apurímac", "Arequipa", "Ayacucho", "Cajamarca", "Callao", "Cusco", "Huancavelica",
         "Huánuco", "Ica", "Junín", "La Libertad", "Lambayeque", "Lima", "Loreto", "Madre de Dios", "Moquegua", "Pasco",
         "Piura", "Puno", "San Martín", "Tacna", "Tumbes", "Ucayali")
Encoding(DEP) <- "UTF-8"
REG <- c("Lima Metropolitana", "Resto de la Costa", "Sierra", "Selva")
DOM2REG <- c("Resto de la Costa", "Resto de la Costa", "Resto de la Costa", "Sierra", "Sierra", "Sierra", "Selva", "Lima Metropolitana")

leer <- function(path) {
  if (!file.exists(path)) {
    alternative <- file.path(dirname(dirname(path)), basename(path))
    if (file.exists(alternative)) path <- alternative else stop(paste("Archivo no encontrado:", path))
  }
  first <- readLines(path, n = 1, encoding = "latin1", warn = FALSE)
  sep <- if (lengths(regmatches(first, gregexpr(";", first))) > 5) ";" else ","
  selected <- NULL
  if (grepl("^Enaho", basename(path), ignore.case = TRUE)) {
    available <- names(fread(path, sep = sep, nrows = 0, encoding = "Latin-1"))
    wanted <- c("MES","CONGLOME","UBIGEO","DOMINIO","ESTRATO","RESULT","FACTOR07","FACTOR","P110","P110A","P110C","P110C1","P110C2","P110C3","P110E","P208A",paste0("P402",1:5),paste0("P403",1:9),"P40310","P40311","P40313","P40314",paste0("P409",1:11))
    selected <- intersect(available, wanted)
  }
  d <- fread(path, sep = sep, select = selected, colClasses = "character", encoding = "Latin-1", showProgress = FALSE)
  setnames(d, trimws(sub("^|^", "", names(d))))
  d
}
num <- function(x) suppressWarnings(as.numeric(gsub(",", ".", trimws(x), fixed = TRUE)))
rb <- function(...) rbind(..., fill = TRUE)
QA <- list()
ind <- function(cond, dominio) {
  # A missing answer is never a negative answer. Structural skips are defined
  # explicitly in each indicator's domain, not inferred from !is.na alone.
  dom <- rep_len(dominio, length(cond)) %in% TRUE
  QA[[length(QA) + 1]] <<- data.table(
    condicion = paste(deparse(substitute(cond)), collapse = " "),
    n_dominio = sum(dom), faltantes_en_dominio = sum(dom & is.na(cond)))
  out <- rep(NA_real_, length(cond))
  ok <- dom & !is.na(cond)
  out[ok] <- as.numeric(cond[ok])
  out
}

# ---------------------------------------------------------------------------
# 1. Datos y diseños
# ---------------------------------------------------------------------------
msg("ENDES 2025: hogares")
h0  <- leer(file.path(ENDES, "1038-Modulo1629/1036-Modulo1629/RECH0_2025.csv"))
h23 <- leer(file.path(ENDES, "1038-Modulo1630/1036-Modulo1630/RECH23_2025.csv"))
h0[, HHID := trimws(HHID)]; h23[, HHID := trimws(HHID)]
stopifnot(!anyDuplicated(h0$HHID), !anyDuplicated(h23$HHID))
h <- merge(h0, h23, by = "HHID")
for (v in c("HV001", "HV005", "HV022", "HV024", "HV201", "HV237", "SHREGION", "SH42", "SH48", "SH49", "SH50")) set(h, j = v, value = num(h[[v]]))
h <- h[HV005 > 0]                                    # hogares entrevistados (marco analítico)
h[, `:=`(w = HV005 / 1e6, region = REG[SHREGION], dep = DEP[HV024])]
h[, `:=`(
  fuente_dentro = ind(HV201 == 11, !is.na(HV201) & !HV201 %in% c(98, 99)),
  trata_agua    = ind(HV237 == 1,  HV237 %in% 0:1),
  continuidad   = ind(SH42 == 1,   SH42 %in% 0:1),
  almacena      = ind(SH48 == 1,   SH48 %in% 0:1),
  olla_balde    = ind(SH48 == 1 & SH49 %in% 1:2, SH48 == 0 | (SH48 == 1 & SH49 %in% c(1:4, 96))),
  sin_tapa      = ind(SH48 == 1 & SH50 == 0, SH48 == 0 | (SH48 == 1 & SH50 %in% 0:1))
)]
d_h <- svydesign(ids = ~HV001, strata = ~HV022, weights = ~w, nest = TRUE, data = h)

msg("ENDES 2025: mujeres")
w0  <- leer(file.path(ENDES, "1038-Modulo1631/1036-Modulo1631/REC0111_2025.csv"))
r91 <- leer(file.path(ENDES, "1038-Modulo1631/1036-Modulo1631/REC91_2025.csv"))
r75 <- leer(file.path(ENDES, "1038-Modulo1636/1036-Modulo1636/RE758081_2025.csv"))
w0[, CASEID := trimws(CASEID)]; r91[, CASEID := trimws(CASEID)]; r75[, CASEID := trimws(CASEID)]
m <- merge(merge(w0, r91, by = "CASEID", all.x = TRUE), r75, by = "CASEID", all.x = TRUE)
for (v in c("V001", "V005", "V022", "V024", "SREGION", "S481AA", "S481AB", "S481AC", "V751", "V754JP")) set(m, j = v, value = num(m[[v]]))
m <- m[V005 > 0]
m[, `:=`(w = V005 / 1e6, region = REG[SREGION], dep = DEP[V024])]
m[, `:=`(
  S481AB       = ind(S481AB == 1, S481AB %in% 0:1),
  S481AC       = ind(S481AC == 1, S481AC %in% 0:1),
  vih_mosquito = ind(V754JP == 1, V751 == 1 & V754JP %in% c(0, 1, 8)), # «no sabe» (8) es una categoría explícita
  vih_mosquito_sin_ns = ind(V754JP == 1, V751 %in% 1 & V754JP %in% 0:1) # sensibilidad: excluye «no sabe»
)]
d_m <- svydesign(ids = ~V001, strata = ~V022, weights = ~w, nest = TRUE, data = m)

msg("ENAHO: hogares y personas")
prep_h <- function(d, peso) {
  for (v in c("MES", "CONGLOME", "DOMINIO", "ESTRATO", "RESULT", "P110", "P110A", "P110C", "P110C1", "P110C2", "P110C3", "P110E", peso)) set(d, j = v, value = num(d[[v]]))
  d <- d[RESULT %in% 1:2]
  d[, `:=`(w = get(peso), strat = DOMINIO * 10 + ESTRATO, region = DOM2REG[DOMINIO],
           dep = DEP[as.integer(substr(sprintf("%06s", trimws(UBIGEO)), 1, 2))],
           strat_dep = as.integer(substr(sprintf("%06s", trimws(UBIGEO)), 1, 2)) * 10 + ESTRATO,
           trim = (MES - 1) %/% 3 + 1)]
  red <- d$P110 %in% 1:2 & d$P110C %in% 1:2
  hrs <- ifelse(d$P110C == 1, d$P110C1, d$P110C3); dias <- ifelse(d$P110C == 1, 7, d$P110C2)
  d[, `:=`(
    red_dentro       = ind(P110 == 1, P110 %in% 1:8),
    todos_dias       = ind(P110C == 1, red),
    agua_24h         = ind(P110C == 1 & P110C1 == 24, red & (P110C == 2 | P110C1 %in% 1:24)),
    horas_dia_media  = ifelse(red & hrs >= 1 & hrs <= 24 & dias >= 1 & dias <= 7, hrs * dias / 7, NA_real_),
    cloro_seguro     = ind(P110A == 1, P110A %in% 1:3),
    sin_cloro        = ind(P110A == 3, P110A %in% 1:3),
    muestra_deposito = ind(P110E %in% 2:4, P110E %in% 1:7)
  )]
  d
}
RZ <- c("P4091", "P4092", "P4093", "P4094", "P4095", "P4096", "P4097", "P4098", "P4099", "P40910", "P40911")
prep_p <- function(d, peso) {
  vs <- c("MES", "CONGLOME", "DOMINIO", "ESTRATO", "P208A", paste0("P402", 1:5), paste0("P403", 1:9), "P40310", "P40311", "P40313", "P40314", RZ, peso)
  for (v in vs) set(d, j = v, value = num(d[[v]]))
  d <- d[get(peso) > 0]
  d[, `:=`(w = get(peso), strat = DOMINIO * 10 + ESTRATO, region = DOM2REG[DOMINIO],
           dep = DEP[as.integer(substr(sprintf("%06s", trimws(UBIGEO)), 1, 2))], trim = (MES - 1) %/% 3 + 1)]
  ill   <- rowSums(d[, paste0("P402", 1:4), with = FALSE] == 1, na.rm = TRUE) > 0
  estab <- rowSums(d[, paste0("P403", 1:9), with = FALSE] == 1, na.rm = TRUE) > 0
  anyc  <- rowSums(d[, c(paste0("P403", 1:9), "P40310", "P40311", "P40313", "P40314"), with = FALSE] == 1, na.rm = TRUE) > 0
  nino_enf <- d$P208A <= 11 & ill & anyc
  d[, `:=`(
    acudio_establecimiento = ind(estab, nino_enf),
    farmacia               = ind(P40310 == 1, nino_enf),
    no_busco               = ind(P40314 == 1, nino_enf)
  )]
  anyr <- rowSums(d[, RZ, with = FALSE] == 1, na.rm = TRUE) > 0
  for (r in RZ) set(d, j = paste0("r_", r), value = ind(d[[r]] == 1, nino_enf & !estab & anyr))
  d
}
e25 <- prep_h(leer(file.path(ENAHO, "1031-Modulo01/1031-Modulo01/Enaho01-2025-100.csv")), "FACTOR07"); e25[, anio := 2025]
q25 <- prep_p(leer(file.path(ENAHO, "1031-Modulo04/1031-Modulo04/Enaho01A-2025-400.csv")), "FACTOR07"); q25[, anio := 2025]
e26 <- rbind(prep_h(leer(file.path(ENAHO, "1039-Modulo01/1039-Modulo01/Enaho01-2026-100.csv")), "FACTOR"),
             prep_h(leer(file.path(ENAHO, "1048-Modulo01/1048-Modulo01/Enaho01-2026-100.csv")), "FACTOR")); e26[, anio := 2026]
q26 <- rbind(prep_p(leer(file.path(ENAHO, "1039-Modulo04/1039-Modulo04/Enaho01A-2026-400.csv")), "FACTOR"),
             prep_p(leer(file.path(ENAHO, "1048-Modulo04/1048-Modulo04/Enaho01A-2026-400.csv")), "FACTOR")); q26[, anio := 2026]
dz <- function(d, strata = ~strat) svydesign(ids = ~CONGLOME, strata = strata, weights = ~w, nest = TRUE, data = d)
d_e25 <- dz(e25); d_q25 <- dz(q25)
# Agrupación exploratoria de dos trimestres con FACTOR/2. La escala constante
# no cambia proporciones/medias; NO acredita calibración semestral ni el diseño.
# Debe cotejarse con la documentación oficial antes de inferencia definitiva.
e26[, w := w / 2]; q26[, w := w / 2]
d_e26 <- dz(e26); d_q26 <- dz(q26)

# Exportación individual opcional SOLO para trabajo local. Nunca se publica.
if (identical(Sys.getenv("EXPORT_ANALITICO"), "1")) {
AN <- file.path(OUT, "analitico"); dir.create(AN, showWarnings = FALSE)
fwrite(h[, c("HHID", "HV001", "HV022", "w", "region", "dep", HOG_ <- c("fuente_dentro", "continuidad", "almacena", "olla_balde", "sin_tapa", "trata_agua")), with = FALSE], file.path(AN, "endes25_hogares.csv"))
fwrite(m[, c("CASEID", "V001", "V022", "w", "region", "dep", "S481AB", "S481AC", "vih_mosquito", "vih_mosquito_sin_ns"), with = FALSE], file.path(AN, "endes25_mujeres.csv"))
EHc <- c("CONGLOME", "strat", "w", "region", "dep", "trim", "red_dentro", "todos_dias", "agua_24h", "horas_dia_media", "cloro_seguro", "sin_cloro", "muestra_deposito")
EPc <- c("CONGLOME", "strat", "w", "region", "dep", "trim", "acudio_establecimiento", "farmacia", "no_busco")
fwrite(e25[, EHc, with = FALSE], file.path(AN, "enaho25_hogares.csv")); fwrite(q25[, EPc, with = FALSE], file.path(AN, "enaho25_personas.csv"))
fwrite(e26[, EHc, with = FALSE], file.path(AN, "enaho26_hogares.csv")); fwrite(q26[, EPc, with = FALSE], file.path(AN, "enaho26_personas.csv"))
}

# ---------------------------------------------------------------------------
# 2. Estimación por dominio
# ---------------------------------------------------------------------------
est1 <- function(des, y, media = FALSE) {
  f <- as.formula(paste0("~", y))
  n <- sum(!is.na(des$variables[[y]]))
  if (n < 2) return(data.table(n = n))
  sm <- svymean(f, des, na.rm = TRUE, deff = "replace")
  p <- as.numeric(coef(sm))[1]; se <- as.numeric(SE(sm))[1]; gl <- degf(des)
  tq <- qt(0.975, max(gl, 1))
  lw <- p - tq * se; hw <- p + tq * se
  if (media) { lo <- lw; hi <- hw } else {
    if (p <= 0 || p >= 1) { lo <- NA_real_; hi <- NA_real_ } else {
      ci <- attr(svyciprop(f, des, method = "logit", df = gl, na.rm = TRUE), "ci")
      lo <- as.numeric(ci)[1]; hi <- as.numeric(ci)[2]
    }
    lw <- max(0, lw); hw <- min(1, hw)
  }
  data.table(n = n, p = p, se = se, deff = as.numeric(deff(sm))[1], gl = gl, lo = lo, hi = hi, lo_wald = lw, hi_wald = hw)
}
estimar <- function(des, vars, encuesta, periodo, medias = character(), niveles = c("region", "dep", "nacional"), dom = NULL) {
  if (!is.null(dom)) des <- subset(des, eval(dom, des$variables))
  res <- list()
  for (y in vars) {
    dy <- subset(des, !is.na(des$variables[[y]]))
    for (nv in niveles) {
      grupos <- if (nv == "nacional") "Perú" else sort(unique(na.omit(dy$variables[[nv]])))
      for (g in grupos) {
        dg <- if (nv == "nacional") dy else subset(dy, dy$variables[[nv]] == g)
        res[[length(res) + 1]] <- cbind(data.table(encuesta = encuesta, periodo = periodo, indicador = y,
          nivel = c(region = "Región natural", dep = "Departamento", nacional = "Nacional")[nv], grupo = g),
          est1(dg, y, media = y %in% medias))
      }
    }
  }
  rbindlist(res, fill = TRUE)
}
msg("Estimaciones ENDES")
HOG <- c("fuente_dentro", "continuidad", "almacena", "olla_balde", "sin_tapa", "trata_agua")
MUJ <- c("S481AB", "S481AC", "vih_mosquito", "vih_mosquito_sin_ns")
E <- rb(estimar(d_h, HOG, "ENDES 2025", "Año 2025"), estimar(d_m, MUJ, "ENDES 2025", "Año 2025"))
msg("Estimaciones ENAHO 2025")
EH <- c("red_dentro", "todos_dias", "agua_24h", "cloro_seguro", "sin_cloro", "muestra_deposito", "horas_dia_media")
EP <- c("acudio_establecimiento", "farmacia", "no_busco")
N25 <- rb(estimar(d_e25, EH, "ENAHO 2025", "Año 2025", medias = "horas_dia_media"),
             estimar(d_q25, EP, "ENAHO 2025", "Año 2025"))
# No se estiman trimestres 2025 con el factor anual.
msg("Estimaciones ENAHO 2026")
N26 <- rb(estimar(d_e26, EH, "ENAHO 2026", "Semestre 2026", medias = "horas_dia_media"),
             estimar(d_q26, EP, "ENAHO 2026", "Semestre 2026"))
for (q in 1:2) {
  dq_e <- dz(e26[trim == q][, w := w * 2]); dq_q <- dz(q26[trim == q][, w := w * 2])   # cada trimestre con su propio diseño
  N26 <- rb(N26, estimar(dq_e, EH, "ENAHO 2026", paste0("T", q, " 2026"), medias = "horas_dia_media", niveles = c("region", "nacional")),
                    estimar(dq_q, EP, "ENAHO 2026", paste0("T", q, " 2026"), niveles = c("region", "nacional")))
}
RZT <- rb(estimar(d_q25, paste0("r_", RZ), "ENAHO 2025", "Año 2025", niveles = "nacional"),
             estimar(d_q26, paste0("r_", RZ), "ENAHO 2026", "Semestre 2026", niveles = "nacional"))
EST <- rbind(E, N25, N26, RZT, fill = TRUE)
# Precisión: coeficiente de variación (EE / estimación). Criterio INEI: CV > 15% = cifra referencial.
EST[, cv := se / p]
EST[, precision := fifelse(n < 30, "insuficiente (n < 30)", fifelse(cv > 0.30, "muy imprecisa (CV > 30%)", fifelse(cv > 0.15, "referencial (CV 15-30%)", "confiable (CV <= 15%)")))]
fwrite(EST, file.path(OUT, "r_estimaciones.csv"))

# ---------------------------------------------------------------------------
# 3. Rao-Scott (segundo orden, F) y V de Cramér
# ---------------------------------------------------------------------------
msg("Rao-Scott")
rs <- function(des, y, g, fuente, psu, est) {
  dy <- subset(des, !is.na(des$variables[[y]]) & !is.na(des$variables[[g]]))
  # Variante con los datos filtrados antes del diseño (como Python y samplics), para comparar implementaciones
  df_f <- des$variables[!is.na(des$variables[[y]]) & !is.na(des$variables[[g]]), ]
  df_f$.y <- factor(df_f[[y]]); df_f$.g <- factor(df_f[[g]])
  dff <- svydesign(ids = reformulate(psu), strata = reformulate(est), weights = ~w, nest = TRUE, data = df_f)
  tf <- svychisq(~.y + .g, dff, statistic = "F")
  dy$variables$.y <- factor(dy$variables[[y]]); dy$variables$.g <- factor(dy$variables[[g]])
  tst <- svychisq(~.y + .g, dy, statistic = "F")
  n <- sum(weights(dy, "sampling") > 0)
  tab <- svytable(~.y + .g, dy, Ntotal = n)
  x2 <- suppressWarnings(chisq.test(tab, correct = FALSE)$statistic)
  data.table(fuente = fuente, indicador = y, nivel = c(region = "Región natural", dep = "Departamento")[g], n = n,
             F = unname(tst$statistic), gl1 = unname(tst$parameter[1]), gl2 = unname(tst$parameter[2]), p = tst$p.value,
             cramer_v = sqrt(unname(x2) / (n * (min(dim(tab)) - 1))),
             F_filtrado = unname(tf$statistic), gl1_filtrado = unname(tf$parameter[1]), gl2_filtrado = unname(tf$parameter[2]), p_filtrado = tf$p.value)
}
RS <- rbindlist(c(
  lapply(HOG, function(y) rbind(rs(d_h, y, "region", "ENDES 2025", "HV001", "HV022"), rs(d_h, y, "dep", "ENDES 2025", "HV001", "HV022"))),
  lapply(MUJ, function(y) rbind(rs(d_m, y, "region", "ENDES 2025", "V001", "V022"), rs(d_m, y, "dep", "ENDES 2025", "V001", "V022"))),
  lapply(setdiff(EH, "horas_dia_media"), function(y) rbind(rs(d_e25, y, "region", "ENAHO 2025", "CONGLOME", "strat"), rs(d_e25, y, "dep", "ENAHO 2025", "CONGLOME", "strat"))),
  lapply(EP, function(y) rbind(rs(d_q25, y, "region", "ENAHO 2025", "CONGLOME", "strat"), rs(d_q25, y, "dep", "ENAHO 2025", "CONGLOME", "strat")))))
RS[, efecto := cut(cramer_v, c(-Inf, 0.1, 0.3, 0.5, Inf), right = FALSE, labels = c("despreciable", "pequeño", "moderado", "grande"))]
fwrite(RS, file.path(OUT, "r_rao_scott.csv"))

# ---------------------------------------------------------------------------
# 4. Comparación temporal: ejecutar contrastes_temporales.R.
# Usa referencia anual 2025 y trimestres 2026, con ajuste por mes.

# Trimestre I frente a II de 2026: sin conglomerados comunes, son muestras independientes.
msg("Trimestre I frente a II de 2026")
tq <- list()
for (y in c("agua_24h", "todos_dias", "cloro_seguro", "horas_dia_media")) for (g in c(REG, "Perú")) {
  dd <- subset(d_e26, !is.na(d_e26$variables[[y]]) & (g == "Perú" | d_e26$variables$region == g))
  b <- svyby(as.formula(paste0("~", y)), ~trim, dd, svymean, covmat = TRUE); V <- vcov(b); est <- coef(b)
  d <- unname(est[2] - est[1]); se <- sqrt(V[1, 1] + V[2, 2] - 2 * V[1, 2]); gl <- degf(dd)
  tq[[length(tq) + 1]] <- data.table(indicador = y, region = g, trim1 = unname(est[1]), trim2 = unname(est[2]), diferencia = d, se = se,
    ic_bajo = d - qt(0.975, gl) * se, ic_alto = d + qt(0.975, gl) * se, p_valor = 2 * pt(-abs(d / se), gl))
}
fwrite(rbindlist(tq), file.path(OUT, "r_trimestres_2026.csv"))

# ---------------------------------------------------------------------------
# 5. Análisis ecológico con las estimaciones departamentales de R
# ---------------------------------------------------------------------------
if (identical(Sys.getenv("INCLUIR_ECOLOGICO"), "1")) {
msg("Análisis ecológico exploratorio, no conciliado")
carga <- fread(CARGA, encoding = "UTF-8")
stopifnot(all(DEP %in% carga$departamento))
IND_ECO <- list(c("ENDES 2025", "continuidad"), c("ENDES 2025", "olla_balde"), c("ENDES 2025", "fuente_dentro"),
                c("ENAHO 2025", "agua_24h"), c("ENAHO 2025", "cloro_seguro"), c("ENAHO 2025", "sin_cloro"),
                c("ENAHO 2025", "muestra_deposito"), c("ENAHO 2025", "no_busco"))
eco <- list()
for (k in IND_ECO) {
  dd <- merge(EST[encuesta == k[1] & periodo == "Año 2025" & indicador == k[2] & nivel == "Departamento", .(departamento = grupo, p)], carga, by = "departamento")
  stopifnot(nrow(dd) == 25)
  for (gr in c("tia25_alta", "top26")) {
    a <- dd[get(gr) == 1, p]; b <- dd[get(gr) == 0, p]
    wt <- wilcox.test(a, b, exact = TRUE)
    eco[[length(eco) + 1]] <- data.table(fuente = k[1], indicador = k[2], grupo = gr, n_grupo = length(a), n_resto = length(b),
      U = unname(wt$statistic), p = wt$p.value, r_biserial = 2 * unname(wt$statistic) / (length(a) * length(b)) - 1)
  }
}
ECO <- rbindlist(eco)[, p_holm := p.adjust(p, "holm")]
fwrite(ECO, file.path(OUT, "r_ecologico.csv"))
}

# ---------------------------------------------------------------------------
# 6. Sensibilidad de la especificación
# ---------------------------------------------------------------------------
msg("Sensibilidad")
sens <- list()
# 6a. Política para estratos con una sola UPM (departamentos)
for (pol in c("adjust", "certainty", "average")) {
  options(survey.lonely.psu = pol)
  for (y in c("continuidad", "olla_balde")) for (g in DEP) {
    dg <- subset(d_h, !is.na(d_h$variables[[y]]) & d_h$variables$dep == g)
    sm <- svymean(as.formula(paste0("~", y)), dg)
    sens[[length(sens) + 1]] <- data.table(prueba = "UPM única", variante = pol, encuesta = "ENDES 2025", indicador = y, grupo = g, p = as.numeric(coef(sm))[1], se = as.numeric(SE(sm))[1])
  }
}
options(survey.lonely.psu = LONELY)
# 6b. Filtrar antes del diseño (como los scripts de Python) frente a subset() del marco completo
for (y in HOG) for (g in REG) {
  sub <- h[!is.na(get(y)) & region == g]
  dpre <- svydesign(ids = ~HV001, strata = ~HV022, weights = ~w, nest = TRUE, data = sub)
  options(survey.lonely.psu = "certainty"); smp <- svymean(as.formula(paste0("~", y)), dpre); options(survey.lonely.psu = LONELY)
  sens[[length(sens) + 1]] <- data.table(prueba = "Filtrar antes del diseño", variante = "filtrado, UPM única con varianza cero (Python)", encuesta = "ENDES 2025", indicador = y, grupo = g, p = as.numeric(coef(smp))[1], se = as.numeric(SE(smp))[1])
}
# 6c. Estratos ENAHO: DOMINIO×ESTRATO (principal) frente a departamento×ESTRATO
d_e25b <- dz(e25, ~strat_dep)
for (y in c("agua_24h", "cloro_seguro")) for (g in REG) {
  dg <- subset(d_e25b, !is.na(d_e25b$variables[[y]]) & d_e25b$variables$region == g)
  sm <- svymean(as.formula(paste0("~", y)), dg)
  sens[[length(sens) + 1]] <- data.table(prueba = "Estratos ENAHO", variante = "departamento × ESTRATO", encuesta = "ENAHO 2025", indicador = y, grupo = g, p = as.numeric(coef(sm))[1], se = as.numeric(SE(sm))[1])
}
fwrite(rbindlist(sens), file.path(OUT, "r_sensibilidad.csv"))
fwrite(rbindlist(QA), file.path(OUT, "r_control_faltantes.csv"))

writeLines(c(R.version.string, paste("survey", packageVersion("survey")), paste("data.table", packageVersion("data.table")),
             paste("lonely.psu:", LONELY), format(Sys.time(), "%Y-%m-%d %H:%M %Z")), file.path(OUT, "r_sesion.txt"))
msg("Listo: ", OUT)