Auditoría de aleatoriedad por posiciones y precisión por tamaño muestral

Objetivo del documento

Este documento ejecuta, sin interfaz gráfica, dos módulos:

Aleatoriedad de la selección por posiciones (IDs en el orden del marco): se compara la muestra observada con una SRSWOR del mismo tamaño usando simulación Monte Carlo (bloques chi², KS de posiciones, rachas en 0/1 y brecha máxima).

Precisión de la media bajo corrección por población finita (CPF), efecto de diseño (EFD) y tasa de respuesta (TR), sin requerir valores de la variable de estudio en toda la población. Si existe score en la prueba piloto, se usa su varianza; de lo contrario, se emplea el acotamiento conservador S²_ref = (U−L)²/4.

Librerías y utilidades

# Librerías (solo necesarias para lectura rápida de archivos)
suppressPackageStartupMessages({
  if (!requireNamespace("readr", quietly = TRUE))  message("Sugerencia: instalar 'readr' para lectura rápida de CSV/TXT.")
  if (!requireNamespace("readxl", quietly = TRUE)) message("Sugerencia: instalar 'readxl' para leer XLSX.")
})

# Lectura flexible por extensión
leer_archivo <- function(path, delim = ",", sheet = 1L) {
  ext <- tolower(tools::file_ext(path))
  if (ext %in% c("csv")) {
    if (requireNamespace("readr", quietly = TRUE)) {
      as.data.frame(readr::read_csv(path, show_col_types = FALSE, progress = FALSE))
    } else {
      utils::read.csv(path, stringsAsFactors = FALSE)
    }
  } else if (ext %in% c("txt")) {
    if (requireNamespace("readr", quietly = TRUE)) {
      as.data.frame(readr::read_delim(path, delim = delim, show_col_types = FALSE, progress = FALSE))
    } else {
      utils::read.table(path, header = TRUE, sep = delim, stringsAsFactors = FALSE)
    }
  } else if (ext %in% c("xlsx")) {
    if (!requireNamespace("readxl", quietly = TRUE)) stop("Para leer XLSX instale 'readxl' (install.packages('readxl')).")
    as.data.frame(readxl::read_excel(path, sheet = sheet))
  } else {
    stop("Formato no soportado: ", ext)
  }
}

# ---------------- Funciones estadísticas (base R) ----------------

# Mapea IDs seleccionados a posiciones 1..N del marco
map_a_pos <- function(ids_sel, ids_marco) {
  pos <- match(ids_sel, ids_marco)
  pos <- pos[!is.na(pos) & pos >= 1 & pos <= length(ids_marco)]
  sort(unique(pos))
}

# Nº de rachas en vector binario (1 si seleccionado, 0 en caso contrario)
runs_count <- function(idx, N) {
  z <- integer(N); z[idx] <- 1L
  sum(c(1L, diff(z)) != 0L)
}

# Máxima brecha (nº consecutivo de no-seleccionados) entre seleccionados
max_gap <- function(idx, N) {
  gaps <- diff(c(0, sort(idx), N+1L)) - 1L
  max(gaps)
}

# KS sobre posiciones normalizadas u = pos/N; compara ECDF(u) vs Uniforme(0,1)
ks_pos_stat <- function(pos, N) {
  u <- sort(pos) / N
  n <- length(u)
  max( max(abs((1:n)/n - u)), max(abs(((1:n)-1)/n - u)) )
}

# Chi² por bloques: divide 1..N en K bloques y compara frecuencias observadas vs esperadas
chi_blocks_stat <- function(pos, N, K) {
  n <- length(pos)
  bloques <- cut(1:N, breaks = K, labels = FALSE)
  obs <- tabulate(bloques[pos], nbins = K)
  esp <- rep(n/K, K)
  sum((obs - esp)^2 / esp)
}

# Combinación por Fisher de p-valores independientes (aprox. aquí)
p_fisher <- function(pvals) {
  X <- -2 * sum(log(pvals))
  1 - pchisq(X, df = 2*length(pvals))
}

# Diagnóstico completo por Monte Carlo sobre posiciones
diag_posiciones <- function(N, pos, K = 10L, B = 5000L, semilla = NULL) {
  if (!is.null(semilla)) set.seed(semilla)
  n <- length(pos); pos <- sort(pos)
  T_chi <- chi_blocks_stat(pos, N, K)
  T_ks  <- ks_pos_stat(pos, N)
  T_run <- runs_count(pos, N)
  T_gap <- max_gap(pos, N)

  chi_b <- numeric(B); ks_b <- numeric(B); run_b <- numeric(B); gap_b <- numeric(B)
  for (b in seq_len(B)) {
    pb <- sort(sample.int(N, n, replace = FALSE))
    chi_b[b] <- chi_blocks_stat(pb, N, K)
    ks_b[b]  <- ks_pos_stat(pb, N)
    run_b[b] <- runs_count(pb, N)
    gap_b[b] <- max_gap(pb, N)
  }

  p_chi <- mean(chi_b >= T_chi)
  p_ks  <- mean(ks_b  >= T_ks)
  p_run <- 2*min(mean(run_b <= T_run), mean(run_b >= T_run)); p_run <- min(1, p_run)
  p_gap <- mean(gap_b >= T_gap)

  list(
    n = n, N = N, K = K,
    T_obs = c(chi = T_chi, ks = T_ks, runs = T_run, gap = T_gap),
    pvals = c(p_chi = p_chi, p_ks = p_ks, p_runs = p_run, p_gap = p_gap),
    p_fisher = p_fisher(c(p_chi, p_ks, p_run, p_gap)),
    null = list(chi = chi_b, ks = ks_b, runs = run_b, gap = gap_b),
    bloques_obs = tabulate(cut(1:N, K, labels = FALSE)[pos], nbins = K),
    pos = pos
  )
}

# Planificación de n para la media (CPF + EFD + TR)
plan_n <- function(N, S2, EA, alpha = 0.05, EFD = 1, TR = 1) {
  z  <- stats::qnorm(1 - alpha/2)
  num <- N * EFD * (z^2) * S2
  den <- (EA^2) * N + EFD * (z^2) * S2
  n   <- num / den
  n_obj <- ceiling(n)
  n_ctc <- ceiling(n / TR)
  semiancho <- z * sqrt( (1 - n_obj/N) * (EFD * S2) / n_obj )
  list(n_objetivo = n_obj, n_contacto = n_ctc, z = z, semiancho_teorico = semiancho)
}

# Curva de semiancho del IC95% (media) en función de n
curva_semiancho <- function(N, S2, alpha = 0.05, EFD = 1, n_max = N-1) {
  z  <- stats::qnorm(1 - alpha/2)
  ns <- unique(c(seq(20, min(600, n_max), by=10), seq(min(620, n_max), n_max, by=20)))
  sem <- z * sqrt( (1 - ns/N) * (EFD * S2) / ns )
  data.frame(n = ns, semiancho = sem)
}

Carga/Simulación de datos

set.seed(params$semilla)

if (params$modo == "simular") {
  # Simula IDs del marco y selecciona posiciones de piloto y muestra
  N <- as.integer(params$N)
  ids_marco <- sprintf("ESC-%04d", 1:N)
  id_pil <- sort(sample.int(N, size = as.integer(params$n_pil), replace = FALSE))
  id_mue <- sort(sample.int(N, size = as.integer(params$n_mue), replace = FALSE))

  if (isTRUE(params$incluye_score)) {
    y_all <- round(100 * stats::rbeta(N, 2.2, 2.8), 2)  # 0..100
    piloto  <- data.frame(id = ids_marco[id_pil],  score = y_all[id_pil])
    muestra <- data.frame(id = ids_marco[id_mue],  score = y_all[id_mue])
  } else {
    piloto  <- data.frame(id = ids_marco[id_pil])
    muestra <- data.frame(id = ids_marco[id_mue])
  }
  marco <- data.frame(id = ids_marco)
  id_col <- "id"; score_col <- "score"
} else {
  stopifnot(!is.null(params$file_marco), !is.null(params$file_piloto), !is.null(params$file_muestra))
  marco   <- leer_archivo(params$file_marco)
  piloto  <- leer_archivo(params$file_piloto)
  muestra <- leer_archivo(params$file_muestra)
  id_col <- params$id_col; score_col <- params$score_col
}

# Verificación de columnas
stopifnot(id_col %in% names(marco), id_col %in% names(piloto), id_col %in% names(muestra))

# Vector de IDs del marco y posiciones para piloto/muestra
ids_marco   <- as.character(marco[[id_col]])
ids_piloto  <- as.character(piloto[[id_col]])
ids_muestra <- as.character(muestra[[id_col]])

N <- length(ids_marco)
pos_pil <- map_a_pos(ids_piloto, ids_marco)
pos_mue <- map_a_pos(ids_muestra, ids_marco)

cat("N (marco):", N, "\n")

## N (marco): 800

cat("n_piloto :", length(pos_pil), "\n")

## n_piloto : 60

cat("n_muestra:", length(pos_mue), "\n")

## n_muestra: 220

Aleatoriedad por posiciones (Monte Carlo)

K <- as.integer(params$K); B <- as.integer(params$B)

diag_pil <- diag_posiciones(N, pos_pil, K = K, B = B, semilla = params$semilla)
diag_mue <- diag_posiciones(N, pos_mue, K = K, B = B, semilla = params$semilla)

# Tabla de p-valores
tab_p <- rbind(
  cbind(Conjunto="Piloto",  t(round(diag_pil$pvals,4)), p_Fisher = round(diag_pil$p_fisher,4)),
  cbind(Conjunto="Muestra", t(round(diag_mue$pvals,4)), p_Fisher = round(diag_mue$p_fisher,4))
)
tab_p

##      Conjunto  p_chi    p_ks     p_runs   p_gap    p_Fisher
## [1,] "Piloto"  "0.5972" "0.9156" "1"      "0.8094" "0.9903"
## [2,] "Muestra" "0.3486" "0.1288" "0.5056" "0.0526" "0.0969"

4.1 Cobertura de posiciones

par(mfrow=c(1,2), mar=c(4,3,3,1))
plot(1:N, rep(0, N), type="n", yaxt="n", xlab="ID (1..N)", ylab="", main="Cobertura PILOTO")
points(pos_pil, rep(0, length(pos_pil)), pch=15, col="red")

plot(1:N, rep(0, N), type="n", yaxt="n", xlab="ID (1..N)", ylab="", main="Cobertura MUESTRA")
points(pos_mue, rep(0, length(pos_mue)), pch=15, col="red")

Objetivo del documento

4.2 Bloques (chi²)

par(mfrow=c(1,2), mar=c(4,3,3,1))
barplot(diag_pil$bloques_obs, main="Bloques PILOTO", xlab="Bloque", ylab="Frecuencia")
abline(h = diag_pil$n / diag_pil$K, col="blue", lty=2, lwd=2)

barplot(diag_mue$bloques_obs, main="Bloques MUESTRA", xlab="Bloque", ylab="Frecuencia")
abline(h = diag_mue$n / diag_mue$K, col="blue", lty=2, lwd=2)

4.3 Nulas de rachas y brecha máxima

par(mfrow=c(1,2), mar=c(4,3,3,1))
hist(diag_mue$null$runs, breaks=40, col="gray85", border="white",
     main="Nula de rachas (MUESTRA)", xlab="N° rachas")
abline(v = diag_mue$T_obs["runs"], col="red", lwd=2)

hist(diag_mue$null$gap, breaks=40, col="gray85", border="white",
     main="Nula de brecha (MUESTRA)", xlab="Máxima brecha")
abline(v = diag_mue$T_obs["gap"], col="red", lwd=2)

4.4 ECDF de posiciones vs uniforme y QQ de posiciones

par(mfrow=c(1,2), mar=c(4,3,3,1))
plot(ecdf(pos_pil/N), verticals=TRUE, do.points=FALSE, col="black", lwd=2,
     main="ECDF PILOTO vs Uniforme", xlab="u = pos/N", ylab="F(u)")
abline(0,1,col="red",lty=2,lwd=2)

plot(ecdf(pos_mue/N), verticals=TRUE, do.points=FALSE, col="black", lwd=2,
     main="ECDF MUESTRA vs Uniforme", xlab="u = pos/N", ylab="F(u)")
abline(0,1,col="red",lty=2,lwd=2)

par(mfrow=c(1,2), mar=c(4,3,3,1))
k <- seq_along(pos_pil); exp_pos <- (k * (N + 1)) / (length(pos_pil) + 1)
plot(exp_pos, pos_pil, pch=19, cex=.6, xlab="Posición esperada", ylab="Observada",
     main="QQ posiciones (PILOTO)")
abline(0,1,col="blue",lty=2,lwd=2)

k <- seq_along(pos_mue); exp_pos <- (k * (N + 1)) / (length(pos_mue) + 1)
plot(exp_pos, pos_mue, pch=19, cex=.6, xlab="Posición esperada", ylab="Observada",
     main="QQ posiciones (MUESTRA)")
abline(0,1,col="blue",lty=2,lwd=2)

Precisión de la media (CPF + EFD + TR)

Idea: con la varianza de la prueba piloto (S²_pil) o bien con el acotamiento conservador S²_ref=(U−L)²/4 si no hay piloto, se calcula:

n_objetivo para alcanzar un semiancho objetivo EA (IC95%),

n_contacto = n_objetivo / TR (inflado por no-respuesta),

el semiancho logrado con la muestra actual n_actual para contrastar con el objetivo.

# S^2 de referencia
has_score <- (!is.null(score_col)) && (score_col %in% names(piloto)) && (score_col %in% names(muestra))

if (has_score) {
  Y_pil <- suppressWarnings(as.numeric(piloto[[score_col]])); Y_pil <- Y_pil[is.finite(Y_pil)]
  S2_ref <- tryCatch(stats::var(Y_pil), error=function(e) NA_real_)
  fuente_S2 <- "Varianza de la prueba piloto"
} else {
  S2_ref <- NA_real_
  fuente_S2 <- "Acotamiento [L,U]"
}
if (!is.finite(S2_ref) || S2_ref <= 0) {
  L <- params$L; U <- params$U
  if (!is.finite(L) || !is.finite(U) || U <= L) {
    S2_ref <- 1.0
    fuente_S2 <- "Fallback: S2_ref=1 (parámetros inválidos)"
  } else {
    S2_ref <- (U - L)^2 / 4
    fuente_S2 <- sprintf("Acotamiento [%.3g, %.3g]", L, U)
  }
}

alpha <- max(min(params$alpha, 0.49), 1e-6)
EFD   <- max(params$EFD, 1.0)
TR    <- min(max(params$TR, 1e-6), 0.999999)
EA    <- max(params$EA, 1e-6)

plan <- plan_n(N, S2 = S2_ref, EA = EA, alpha = alpha, EFD = EFD, TR = TR)

n_actual <- length(pos_mue)
semiancho_actual <- plan$z * sqrt( (1 - n_actual/N) * (EFD * S2_ref) / n_actual )

tab_plan <- data.frame(
  N = N, S2_referencia = round(S2_ref,3), fuente_S2 = fuente_S2,
  alpha = alpha, EFD = EFD, TR = TR, EA = EA,
  n_actual = n_actual, n_objetivo = plan$n_objetivo, n_contacto = plan$n_contacto,
  semiancho_en_n_actual = round(semiancho_actual,3),
  semiancho_en_n_obj = round(plan$semiancho_teorico,3),
  check.names = FALSE
)
tab_plan

##     N S2_referencia                    fuente_S2 alpha EFD  TR EA n_actual
## 1 800       391.384 Varianza de la prueba piloto  0.05 1.3 0.8  2      220
##   n_objetivo n_contacto semiancho_en_n_actual semiancho_en_n_obj
## 1        304        380                 2.538              1.997

5.1 Curva precisión–n y comparación con n_actual / n_objetivo

stopifnot(is.finite(S2_ref), S2_ref > 0, N > 2, is.finite(EFD), EFD >= 1, is.finite(EA), EA > 0)

df <- curva_semiancho(N, S2_ref, alpha, EFD, n_max = N-1)
stopifnot(all(is.finite(df$semiancho)))

plot(df$n, df$semiancho, type="l", lwd=2, xlab="n", ylab="Semiancho IC95% (media)",
     main = sprintf("Curva precisión–n — S2_ref: %s", fuente_S2))
abline(h = EA, col="red", lty=2, lwd=2)
abline(v = plan$n_objetivo, col="blue", lty=3, lwd=2)
abline(v = n_actual, col="darkgreen", lty=3, lwd=2)
points(plan$n_objetivo, plan$semiancho_teorico, pch=19, cex=1.1)
points(n_actual, semiancho_actual, pch=19, cex=1.1, col="darkgreen")
legend("topright",
       legend=c("Semiancho", "EA objetivo", "n objetivo", "n actual"),
       col=c("black","red","blue","darkgreen"),
       lty=c(1,2,3,3), lwd=2, pch=c(NA,NA,19,19), bty="n")
mtext(sprintf("sem(n_actual)=%.3f | sem(n_obj)=%.3f", semiancho_actual, plan$semiancho_teorico),
      side=3, adj=1, cex=0.9)

5.2 Efecto sobre el IC: cuantificación

delta <- semiancho_actual - plan$semiancho_teorico
ratio <- semiancho_actual / plan$semiancho_teorico
tab_ci <- data.frame(
  Indicador = c("Semiancho actual", "Semiancho en n_obj", "Diferencia (actual - objetivo)", "Razón (actual/objetivo)"),
  Valor = round(c(semiancho_actual, plan$semiancho_teorico, delta, ratio), 3)
)
tab_ci

##                        Indicador Valor
## 1               Semiancho actual 2.538
## 2             Semiancho en n_obj 1.997
## 3 Diferencia (actual - objetivo) 0.541
## 4        Razón (actual/objetivo) 1.271

if (semiancho_actual <= EA) {
  cat("\nConclusión: Con n_actual ya se cumple el criterio de precisión (semiancho ≤ EA).\n")
} else {
  cat("\nConclusión: Aún no se cumple el criterio de precisión; aumentar n o revisar EFD/TR/S².\n")
}

## 
## Conclusión: Aún no se cumple el criterio de precisión; aumentar n o revisar EFD/TR/S².