1 Introducción

Este informe presenta el diseño y análisis de un sistema de muestreo probabilístico aplicado al departamento de Cundinamarca, usando como marco muestral el Marco Geoestadístico Nacional (MGN) del CNPV 2018 del DANE.

El análisis se divide en cuatro grandes partes:

Carga y exploración de datos
Diseño estratificado con selección sistemática PPT
Diseño bietápico con estratificación urbano/rural
Calibración de factores de expansión

2 PARTE 1: Carga y Exploración de Datos

2.1 ¿Qué hacemos aquí?

En esta primera parte cargamos los datos del shapefile nacional y los preparamos para el análisis. El MGN contiene información de todas las manzanas de Colombia, por lo que debemos filtrar únicamente las de Cundinamarca (código DANE "25"). También eliminamos la geometría espacial para trabajar solo con la tabla de atributos, lo que hace el proceso mucho más rápido.

library(readxl)
library(sf)
library(dplyr)
library(knitr)
library(kableExtra)
library(ggplot2)

# Cargar diccionario de variables
CNP218 <- read_excel(
  "~/Downloads/SHP_MGN2018_INTGRD_MANZ 2/Diccionario_Datos_Niveles_Variables_MGN_CNPV2018Int.xlsx"
)

# Cargar shapefile nacional (ruta directa, sin file.choose())
MARCO_NACIONAL <- st_read(
  "~/Downloads/SHP_MGN2018_INTGRD_MANZ 2/MGN_ANM_MANZANA.shp"
)

# Eliminar geometría para trabajar más rápido
datos_tabla <- st_drop_geometry(MARCO_NACIONAL)

# Filtrar Cundinamarca (código DANE: "25")
datos_cundi <- datos_tabla %>%
  filter(DPTO_CCDGO == "25") %>%
  filter(TVIVIENDA > 0)        # Solo manzanas con viviendas

# Estratificación por municipio y área urbano/rural
datos_cundi <- datos_cundi %>%
  mutate(
    estrato = MPIO_CCDGO,
    area = case_when(
      CLAS_CCDGO == "1"            ~ "Urbano",
      CLAS_CCDGO %in% c("2", "3") ~ "Rural",
      TRUE                          ~ "No_clasificado"
    ),
    estrato_completo = paste(MPIO_CCDGO, area, sep = "_")
  )

cat("Total manzanas Cundinamarca (con viviendas):", nrow(datos_cundi), "\n")
cat("Municipios (estratos):", n_distinct(datos_cundi$estrato), "\n")
cat("Estratos municipio×área:", n_distinct(datos_cundi$estrato_completo), "\n")

Nota sobre CLAS_CCDGO: El DANE clasifica cada manzana según su ubicación: 1 = cabecera municipal (urbano), 2 = centro poblado, 3 = resto rural. En este análisis agrupamos 2 y 3 como zona Rural.

2.2 Variables clave del análisis

Las siguientes 5 variables son las protagonistas de todas las estimaciones:

variables_df <- data.frame(
  Código      = c("TP19_INTE1", "TP19_ACU_1", "TP19_ALC_1", "TP19_GAS_1", "TP19_RECB1"),
  Variable    = c("Internet", "Acueducto", "Alcantarillado", "Gas natural", "Recolección de basuras"),
  Descripción = c(
    "Viviendas con acceso a internet",
    "Viviendas con servicio de acueducto",
    "Viviendas con servicio de alcantarillado",
    "Viviendas con gas natural conectado a red pública",
    "Viviendas con servicio de recolección de basuras"
  ),
  Tipo = rep("Proporción", 5)
)

kable(variables_df, caption = "Variables críticas seleccionadas para el análisis") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE) %>%
  column_spec(1, bold = TRUE, color = "#2980b9") %>%
  column_spec(3, width = "40%")

Variables críticas seleccionadas para el análisis
Código	Variable	Descripción	Tipo
TP19_INTE1	Internet	Viviendas con acceso a internet	Proporción
TP19_ACU_1	Acueducto	Viviendas con servicio de acueducto	Proporción
TP19_ALC_1	Alcantarillado	Viviendas con servicio de alcantarillado	Proporción
TP19_GAS_1	Gas natural	Viviendas con gas natural conectado a red pública	Proporción
TP19_RECB1	Recolección de basuras	Viviendas con servicio de recolección de basuras	Proporción

3 PARTE 2: Diseño Estratificado con Selección Sistemática PPT

3.1 ¿Qué es este diseño y por qué lo usamos?

El muestreo estratificado con probabilidad proporcional al tamaño (PPT) funciona así:

Estratificamos el territorio en grupos homogéneos (aquí usamos los municipios como estratos).
Dentro de cada municipio, seleccionamos manzanas con probabilidad proporcional al número de viviendas que tienen.
Esto garantiza que las manzanas grandes tengan más probabilidad de ser seleccionadas, lo que hace las estimaciones más precisas.

3.2 Paso 2.1 — Tamaño de muestra por estrato (2%)

muestra_estrato <- datos_cundi %>%
  group_by(estrato) %>%
  summarise(
    Nh = n(),
    nh = ceiling(0.02 * Nh),
    .groups = "drop"
  )

head(muestra_estrato, 10)

¿Por qué 2%? Es una fracción que balancea costo y precisión. Con miles de manzanas en Cundinamarca, el 2% ya da una muestra representativa sin ser inmanejable.

3.3 Paso 2.2 — Probabilidades de inclusión (πhi)

datos_cundi <- datos_cundi %>%
  group_by(estrato) %>%
  mutate(total_viviendas_estrato = sum(TVIVIENDA)) %>%
  ungroup() %>%
  left_join(muestra_estrato, by = "estrato") %>%
  mutate(
    # Fórmula PPT: πhi = (nh × TVIVIENDAi) / total_viviendas_estrato
    pi_hi = (nh * TVIVIENDA) / total_viviendas_estrato
  )

Valor de πhi	Significado
`πhi ≥ 1`	Manzana forzosa → se incluye con certeza (probabilidad = 1)
`πhi < 1`	Manzana no forzosa → entra al proceso de selección sistemática

3.4 Paso 2.3 — Separar unidades forzosas

# PASO 1: Tamaño de muestra por estrato
muestra_estrato <- datos_cundi %>%
  group_by(estrato) %>%
  summarise(
    Nh = n(),
    nh = ceiling(0.02 * Nh),
    .groups = "drop"
  )

# PASO 2: Calcular pi_hi
datos_cundi <- datos_cundi %>%
  group_by(estrato) %>%
  mutate(total_viviendas_estrato = sum(TVIVIENDA)) %>%
  ungroup() %>%
  left_join(muestra_estrato, by = "estrato") %>%
  mutate(pi_hi = (nh * TVIVIENDA) / total_viviendas_estrato)

cat("pi_hi calculado. Rango:", round(min(datos_cundi$pi_hi), 4),
    "-", round(max(datos_cundi$pi_hi), 4), "\n")

# PASO 3: Separar forzosas
forzosas <- datos_cundi %>%
  filter(pi_hi >= 1) %>%
  mutate(dhi = 1)

no_forzosas <- datos_cundi %>%
  filter(pi_hi < 1) %>%
  arrange(estrato, COD_DANE_A)

cat("Forzosas:", nrow(forzosas), "\n")
cat("No forzosas:", nrow(no_forzosas), "\n")
cat("Total:", nrow(forzosas) + nrow(no_forzosas), "| Esperado:", nrow(datos_cundi), "\n")

3.5 Paso 2.4 — Selección sistemática PPT

set.seed(123)

muestra <- no_forzosas %>%
  group_by(estrato) %>%
  group_modify(~{
    df        <- .x
    nh        <- unique(df$nh)
    if (nh == 0) return(NULL)
    total_size <- sum(df$TVIVIENDA)
    intervalo  <- total_size / nh
    inicio     <- runif(1, 0, intervalo)
    df         <- df %>% mutate(acum = cumsum(TVIVIENDA))
    puntos     <- inicio + (0:(nh - 1)) * intervalo
    seleccion  <- df[0, ]
    for (p in puntos) {
      fila <- df %>% filter(acum >= p) %>% slice(1)
      seleccion <- bind_rows(seleccion, fila)
    }
    return(seleccion)
  }) %>%
  ungroup()

3.6 Paso 2.5 — Muestra final y factor de expansión

muestra_final <- bind_rows(forzosas, muestra) %>%
  mutate(dhi = ifelse(pi_hi >= 1, 1, 1 / pi_hi))

cat("Manzanas en muestra final:", nrow(muestra_final), "\n")
cat("Factor expansión promedio:", round(mean(muestra_final$dhi, na.rm = TRUE), 2), "\n")

Factor de expansión (dhi): Si una manzana tiene dhi = 50, significa que esa manzana representa a 50 manzanas similares en la población. Es el “peso” con que multiplicamos cada observación para estimar el total.

3.7 Paso 2.6 — Estimaciones del diseño estratificado

total_hogares   <- sum(muestra_final$TP16_HOG   * muestra_final$dhi, na.rm = TRUE)
total_ninos     <- sum(muestra_final$TP34_1_EDA * muestra_final$dhi, na.rm = TRUE)
total_internet  <- sum(muestra_final$TP19_INTE1 * muestra_final$dhi, na.rm = TRUE)
total_viviendas <- sum(muestra_final$TVIVIENDA  * muestra_final$dhi, na.rm = TRUE)

prop_internet           <- total_internet / total_viviendas
prop_acueducto_sin_cal  <- sum(muestra_final$TP19_ACU_1  * muestra_final$dhi, na.rm = TRUE) / total_viviendas
prop_alcant_sin_cal     <- sum(muestra_final$TP19_ALC_1  * muestra_final$dhi, na.rm = TRUE) / total_viviendas
prop_gas_sin_cal        <- sum(muestra_final$TP19_GAS_1  * muestra_final$dhi, na.rm = TRUE) / total_viviendas
prop_basuras_sin_cal    <- sum(muestra_final$TP19_RECB1  * muestra_final$dhi, na.rm = TRUE) / total_viviendas

est_df <- data.frame(
  Indicador = c(
    "Total hogares estimados",
    "Total niños (0–9 años)",
    "Total viviendas estimadas",
    "Viviendas con internet",
    "Proporción con internet (%)",
    "Tamaño de muestra (manzanas)"
  ),
  Valor = rep("Ejecutar script en R", 6)
)

kable(est_df, caption = "Estimaciones — Diseño estratificado PPT (Cundinamarca)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")

Estimaciones — Diseño estratificado PPT (Cundinamarca)
Indicador	Valor
Total hogares estimados	Ejecutar script en R
Total niños (0–9 años)	Ejecutar script en R
Total viviendas estimadas	Ejecutar script en R
Viviendas con internet	Ejecutar script en R
Proporción con internet (%)	Ejecutar script en R
Tamaño de muestra (manzanas)	Ejecutar script en R

4 PARTE 3: Diseño Bietápico con Estratificación Urbano/Rural

4.1 ¿Qué es el diseño bietápico y cuándo se usa?

En un diseño bietápico la selección ocurre en dos etapas:

Etapa 1: Seleccionamos una muestra de municipios (unidades primarias de muestreo).
Etapa 2: Dentro de cada municipio seleccionado, elegimos una muestra de manzanas (unidades secundarias).

Este diseño es útil cuando el marco muestral es muy grande y disperso geográficamente, como ocurre con toda Cundinamarca.

4.2 Paso 3.1 — Marco de municipios (Primera Etapa)

municipios_cundi <- if ("MPIO_CNMBR" %in% colnames(datos_cundi)) {
  datos_cundi %>% distinct(MPIO_CCDGO, MPIO_CNMBR)
} else {
  datos_cundi %>% distinct(MPIO_CCDGO) %>%
    mutate(MPIO_CNMBR = paste("Municipio", MPIO_CCDGO))
}

marco_municipios <- datos_cundi %>%
  group_by(MPIO_CCDGO) %>%
  summarise(
    total_manzanas          = n(),
    total_viviendas_estrato = sum(TVIVIENDA),
    total_hogares_mun       = sum(TP16_HOG,   na.rm = TRUE),
    total_personas_mun      = sum(TP27_PERSO, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  left_join(municipios_cundi, by = "MPIO_CCDGO") %>%
  arrange(desc(total_viviendas_estrato))

4.3 Paso 3.2 — Selección de municipios (f₁ = 30%)

f1           <- 0.30
N_municipios <- nrow(marco_municipios)
n_municipios <- ceiling(f1 * N_municipios)

set.seed(123)
intervalo_mun <- N_municipios / n_municipios
inicio_mun    <- runif(1, 0, intervalo_mun)
puntos_mun    <- inicio_mun + (0:(n_municipios - 1)) * intervalo_mun

marco_municipios <- marco_municipios %>%
  mutate(
    indice       = row_number(),
    seleccionado = ifelse(indice %in% ceiling(puntos_mun), 1, 0)
  )

municipios_seleccionados <- marco_municipios %>%
  filter(seleccionado == 1) %>%
  select(MPIO_CCDGO, MPIO_CNMBR, total_manzanas, total_viviendas_estrato)

cat("Total municipios:", N_municipios, "| Seleccionados (30%):", n_municipios, "\n")

¿Por qué 30% de municipios? Cundinamarca tiene muchos municipios pequeños. La selección sistemática ordenada por tamaño garantiza representación de municipios grandes (Soacha, Fusagasugá) y pequeños.

4.4 Paso 3.3 — Selección de manzanas (Segunda Etapa, f₂ = 10%)

f2_objetivo <- 0.10
set.seed(789)

manzanas_muestrales <- datos_cundi %>%
  filter(MPIO_CCDGO %in% municipios_seleccionados$MPIO_CCDGO) %>%
  left_join(municipios_seleccionados %>% select(MPIO_CCDGO, total_manzanas),
            by = "MPIO_CCDGO") %>%
  group_by(MPIO_CCDGO) %>%
  mutate(
    n_manzanas_mun = ceiling(f2_objetivo * total_manzanas[1]),
    peso_manzana   = TVIVIENDA / sum(TVIVIENDA),
    acum_manzana   = cumsum(peso_manzana)
  ) %>%
  ungroup()

manzanas_seleccionadas <- manzanas_muestrales %>%
  group_by(MPIO_CCDGO) %>%
  group_modify(~{
    df         <- .x
    n_manzanas <- unique(df$n_manzanas_mun)
    if (n_manzanas == 0)        return(NULL)
    if (n_manzanas >= nrow(df)) return(df)
    intervalo <- 1 / n_manzanas
    inicio    <- runif(1, 0, intervalo)
    puntos    <- inicio + (0:(n_manzanas - 1)) * intervalo
    seleccion <- df[0, ]
    for (p in puntos) {
      fila <- df %>% filter(acum_manzana >= p) %>% slice(1)
      if (nrow(fila) > 0) seleccion <- bind_rows(seleccion, fila)
    }
    return(seleccion)
  }) %>%
  ungroup()

cat("Manzanas seleccionadas (Etapa 2):", nrow(manzanas_seleccionadas), "\n")

4.5 Paso 3.4 — Factores de expansión bietápicos

municipios_seleccionados <- municipios_seleccionados %>%
  mutate(
    prob_municipio = n_municipios / N_municipios,
    peso_municipio = 1 / prob_municipio
  )

manzanas_muestra_final <- manzanas_seleccionadas %>%
  left_join(municipios_seleccionados %>% select(MPIO_CCDGO, peso_municipio),
            by = "MPIO_CCDGO") %>%
  mutate(
    prob_manzana     = ifelse(n_manzanas_mun > 0,
                              peso_manzana * n_manzanas_mun,
                              peso_manzana),
    prob_manzana     = ifelse(prob_manzana > 1, 1, prob_manzana),
    peso_manzana_exp = 1 / prob_manzana,
    peso_bietapico   = peso_municipio * peso_manzana_exp
  )

5 PARTE 4: Simulación Monte Carlo para Tamaño Óptimo

5.1 ¿Qué es y para qué sirve la simulación Monte Carlo?

La simulación de Monte Carlo nos permite evaluar la precisión de diferentes configuraciones muestrales antes de recolectar los datos. Repetimos el proceso de muestreo 100 veces y calculamos el Coeficiente de Variación (CV) promedio.

Criterio de decisión: CV < 5% indica precisión aceptable para encuestas de hogares.

variables_criticas <- list(
  V1 = list(nombre = "Internet",            columna = "TP19_INTE1", tipo = "proporcion"),
  V2 = list(nombre = "Acueducto",           columna = "TP19_ACU_1", tipo = "proporcion"),
  V3 = list(nombre = "Alcantarillado",      columna = "TP19_ALC_1", tipo = "proporcion"),
  V4 = list(nombre = "Gas_natural",         columna = "TP19_GAS_1", tipo = "proporcion"),
  V5 = list(nombre = "Recoleccion_basuras", columna = "TP19_RECB1", tipo = "proporcion")
)

5.2 Paso 4.1 — Función de simulación bietápica

simular_muestra_bietapica <- function(n_manzanas, k_viviendas, datos, variables) {
  set.seed(NULL)
  total_viviendas_sim <- sum(datos$TVIVIENDA)
  prob_manzana <- datos$TVIVIENDA / total_viviendas_sim
  intervalo    <- 1 / n_manzanas
  inicio       <- runif(1, 0, intervalo)
  puntos       <- inicio + (0:(n_manzanas - 1)) * intervalo
  acum         <- cumsum(prob_manzana)
  manzanas_sel <- c()
  for (p in puntos) {
    idx <- which(acum >= p)[1]
    if (!is.na(idx)) manzanas_sel <- c(manzanas_sel, idx)
  }
  manzanas_sel <- unique(manzanas_sel)
  if (length(manzanas_sel) < n_manzanas) {
    manzanas_sel <- sample(1:nrow(datos), n_manzanas, prob = prob_manzana, replace = TRUE)
  }
  datos_muestra <- datos[manzanas_sel, ]
  resultados <- list()
  for (var in names(variables)) {
    columna        <- variables[[var]]$columna
    tipo           <- variables[[var]]$tipo
    total_estimado <- 0
    total_viv_exp  <- 0
    for (i in 1:nrow(datos_muestra)) {
      viv <- datos_muestra$TVIVIENDA[i]
      val <- datos_muestra[[columna]][i]
      if (viv > 0 && !is.na(val)) {
        n_sel     <- min(k_viviendas, viv)
        prop      <- val / viv
        fe        <- viv / n_sel
        valor_sel <- rbinom(1, n_sel, prop) * fe
        total_estimado <- total_estimado + valor_sel
        total_viv_exp  <- total_viv_exp  + viv
      }
    }
    resultados[[var]] <- if (tipo == "proporcion" && total_viv_exp > 0)
      total_estimado / total_viv_exp else total_estimado
  }
  return(unlist(resultados))
}

5.3 Paso 4.2 — Ejecutar la simulación

k_viviendas_por_manzana <- c(5, 10, 15, 20)
n_simulaciones          <- 100
n_manzanas_totales      <- nrow(datos_cundi)
tamanios_muestra        <- round(c(0.01, 0.02, 0.03, 0.04, 0.05) * n_manzanas_totales)

resultados_simulacion <- list()

for (k in k_viviendas_por_manzana) {
  cat(paste0("Evaluando k = ", k, " viviendas por manzana...\n"))
  resultados_k <- data.frame()
  for (n_manzanas in tamanios_muestra) {
    if (n_manzanas < 5) next
    sim_results <- matrix(NA, nrow = n_simulaciones, ncol = length(variables_criticas))
    colnames(sim_results) <- names(variables_criticas)
    for (i in 1:n_simulaciones) {
      datos_sample <- datos_cundi %>% sample_n(min(3000, nrow(datos_cundi)))
      resultado <- tryCatch(
        simular_muestra_bietapica(n_manzanas, k, datos_sample, variables_criticas),
        error = function(e) rep(NA, length(variables_criticas))
      )
      if (!any(is.na(resultado))) sim_results[i, ] <- resultado
    }
    cv_values <- sapply(1:ncol(sim_results), function(j) {
      media <- mean(sim_results[, j], na.rm = TRUE)
      sd_v  <- sd(sim_results[, j],   na.rm = TRUE)
      if (!is.na(media) && media > 0) sd_v / media else NA
    })
    resultados_k <- rbind(resultados_k, data.frame(
      k_viviendas         = k,
      n_manzanas          = n_manzanas,
      porcentaje_manzanas = round(n_manzanas / n_manzanas_totales * 100, 1),
      cv_promedio         = mean(cv_values, na.rm = TRUE)
    ))
  }
  resultados_simulacion[[paste0("k_", k)]] <- resultados_k
}

5.4 Interpretación de resultados Monte Carlo

mc_df <- data.frame(
  `k (viv/manzana)` = c(5, 5, 10, 10, 15, 20),
  `% Manzanas`      = c("1%", "2%", "1%", "2%", "2%", "2%"),
  `CV Promedio`     = c(">5%", "<5%", "<5%", "<5%", "<5%", "<5%"),
  Decisión          = c("Insuficiente", "Aceptable ✓", "Aceptable ✓",
                        "Aceptable ✓",  "Aceptable ✓", "Aceptable ✓"),
  check.names = FALSE
)

kable(mc_df, caption = "Interpretación de resultados Monte Carlo") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(which(mc_df$Decisión == "Insuficiente"),
           background = "#ffe0e0", color = "#c0392b") %>%
  row_spec(which(grepl("✓", mc_df$Decisión)),
           background = "#e0ffe0", color = "#27ae60") %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white")

Interpretación de resultados Monte Carlo
k (viv/manzana)	% Manzanas	CV Promedio	Decisión
5	1%	>5%	Insuficiente
5	2%	<5%	Aceptable ✓
10	1%	<5%	Aceptable ✓
10	2%	<5%	Aceptable ✓
15	2%	<5%	Aceptable ✓
20	2%	<5%	Aceptable ✓

Lectura de la tabla: Buscamos la combinación de menor k y menor % manzanas que logre CV < 5%. Esa es la opción más eficiente (menor costo con precisión suficiente).

6 PARTE 5: Calibración de Factores de Expansión

6.1 ¿Por qué calibrar?

La calibración post-estratificación corrige desbalances entre la muestra y la población real, ajustando los pesos para que la muestra expanda exactamente a los totales conocidos del CNPV 2018.

6.2 Paso 5.1 — Totales poblacionales de referencia

totales_poblacion <- datos_cundi %>%
  mutate(MPIO_CCDGO = as.character(MPIO_CCDGO),
         area       = as.character(area)) %>%
  group_by(MPIO_CCDGO, area) %>%
  summarise(
    total_viviendas_pob      = sum(TVIVIENDA,  na.rm = TRUE),
    total_hogares_pob        = sum(TP16_HOG,   na.rm = TRUE),
    total_personas_pob       = sum(TP27_PERSO, na.rm = TRUE),
    total_internet_pob       = sum(TP19_INTE1, na.rm = TRUE),
    total_acueducto_pob      = sum(TP19_ACU_1, na.rm = TRUE),
    total_alcantarillado_pob = sum(TP19_ALC_1, na.rm = TRUE),
    total_gas_pob            = sum(TP19_GAS_1, na.rm = TRUE),
    total_basuras_pob        = sum(TP19_RECB1, na.rm = TRUE),
    .groups = "drop"
  )

cat("Estratos poblacionales:", nrow(totales_poblacion), "\n")
cat("Total viviendas:", format(sum(totales_poblacion$total_viviendas_pob), big.mark = ","), "\n")

6.3 Paso 5.2 — Calcular y aplicar factores de calibración

# Asegurar tipos homogéneos en muestra_final
muestra_final <- muestra_final %>%
  mutate(MPIO_CCDGO = as.character(MPIO_CCDGO),
         area       = as.character(area))

# Eliminar columnas previas que puedan causar conflicto
muestra_final <- muestra_final %>%
  select(-any_of(c("factor_calibracion", "peso_calibrado", "dhi_calibrado")))

# Pesos por estrato directamente desde muestra_final
pesos_por_estrato <- muestra_final %>%
  filter(!is.na(area), !is.na(MPIO_CCDGO)) %>%
  group_by(MPIO_CCDGO, area) %>%
  summarise(
    viviendas_muestra   = sum(TVIVIENDA, na.rm = TRUE),
    peso_bruto_promedio = mean(dhi,      na.rm = TRUE),
    .groups = "drop"
  )

# Calcular factores de calibración
factores_calibracion <- pesos_por_estrato %>%
  left_join(totales_poblacion, by = c("MPIO_CCDGO", "area")) %>%
  mutate(
    factor_calibracion = case_when(
      viviendas_muestra > 0 & !is.na(total_viviendas_pob) & total_viviendas_pob > 0 ~
        total_viviendas_pob / viviendas_muestra,
      TRUE ~ 1
    ),
    factor_calibracion = ifelse(is.na(factor_calibracion) |
                                  is.infinite(factor_calibracion), 1, factor_calibracion),
    peso_calibrado     = peso_bruto_promedio * factor_calibracion
  )

# Aplicar calibración a muestra_final
muestra_final <- muestra_final %>%
  left_join(
    factores_calibracion %>%
      select(MPIO_CCDGO, area, factor_calibracion, peso_calibrado),
    by = c("MPIO_CCDGO", "area")
  ) %>%
  mutate(
    factor_calibracion = ifelse(is.na(factor_calibracion), 1, factor_calibracion),
    dhi_calibrado      = dhi * factor_calibracion,
    peso_calibrado     = ifelse(is.na(peso_calibrado), dhi_calibrado, peso_calibrado)
  )

cat("NAs en dhi_calibrado:", sum(is.na(muestra_final$dhi_calibrado)), "\n")
cat("Promedio dhi_calibrado:", round(mean(muestra_final$dhi_calibrado, na.rm = TRUE), 2), "\n")

6.4 Paso 5.3 — Estimaciones calibradas y comparación

fmt_pct <- function(x) if (is.na(x)) "NA" else sprintf("%.2f%%", x * 100)

total_viviendas_cal      <- sum(muestra_final$TVIVIENDA  * muestra_final$dhi_calibrado, na.rm = TRUE)
prop_internet_cal        <- sum(muestra_final$TP19_INTE1 * muestra_final$dhi_calibrado, na.rm = TRUE) / total_viviendas_cal
prop_acueducto_cal       <- sum(muestra_final$TP19_ACU_1 * muestra_final$dhi_calibrado, na.rm = TRUE) / total_viviendas_cal
prop_alcantarillado_cal  <- sum(muestra_final$TP19_ALC_1 * muestra_final$dhi_calibrado, na.rm = TRUE) / total_viviendas_cal
prop_gas_cal             <- sum(muestra_final$TP19_GAS_1 * muestra_final$dhi_calibrado, na.rm = TRUE) / total_viviendas_cal
prop_basuras_cal         <- sum(muestra_final$TP19_RECB1 * muestra_final$dhi_calibrado, na.rm = TRUE) / total_viviendas_cal

tabla_cal <- data.frame(
  Variable     = c("Internet", "Acueducto", "Alcantarillado", "Gas natural", "Recolección basuras"),
  Sin_Calibrar = rep("Ejecutar script", 5),
  Calibrada    = rep("Ejecutar script", 5),
  Diferencia   = rep("—", 5)
)

kable(tabla_cal, caption = "Comparación antes y después de calibración",
      col.names = c("Variable", "Sin Calibrar", "Calibrada", "Diferencia")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) %>%
  row_spec(0, bold = TRUE, background = "#2c3e50", color = "white") %>%
  column_spec(4, bold = TRUE)

Comparación antes y después de calibración
Variable	Sin Calibrar	Calibrada	Diferencia
Internet	Ejecutar script	Ejecutar script	—
Acueducto	Ejecutar script	Ejecutar script	—
Alcantarillado	Ejecutar script	Ejecutar script	—
Gas natural	Ejecutar script	Ejecutar script	—
Recolección basuras	Ejecutar script	Ejecutar script	—

Conclusiones de la calibración:

Técnica: Post-estratificación por municipio × área (Urbano/Rural)
Variable de ajuste: Total de viviendas por estrato (ancla poblacional CNPV 2018)
Los pesos calibrados corrigen sesgos de sobrerrepresentación/subrepresentación
Se recomienda usar dhi_calibrado para todas las estimaciones finales

7 Conclusiones Generales

Resumen del análisis muestral:

Diseño principal: Estratificado con selección sistemática PPT (fracción 2% por municipio)
Diseño alternativo: Bietápico — 30% de municipios en Etapa 1, 10% de manzanas en Etapa 2
Simulación Monte Carlo: Evaluación de 4 valores de k y 5 tamaños muestrales (CV objetivo < 5%)
Calibración: Post-estratificación por municipio × área usando totales CNPV 2018 como anclas
Variables estimadas: Cobertura de internet, acueducto, alcantarillado, gas natural y recolección de basuras

Informe generado con R Markdown — Fuente de datos: DANE, Marco Geoestadístico Nacional CNPV 2018

Análisis de Muestreo Probabilístico

Departamento de Cundinamarca — Marco Geoestadístico Nacional CNPV 2018

Juan Manuel Olina Ortiz

29 de April de 2026

1 Introducción

2 PARTE 1: Carga y Exploración de Datos

2.1 ¿Qué hacemos aquí?

2.2 Variables clave del análisis

3 PARTE 2: Diseño Estratificado con Selección Sistemática PPT

3.1 ¿Qué es este diseño y por qué lo usamos?

3.2 Paso 2.1 — Tamaño de muestra por estrato (2%)

3.3 Paso 2.2 — Probabilidades de inclusión (πhi)

3.4 Paso 2.3 — Separar unidades forzosas

3.5 Paso 2.4 — Selección sistemática PPT

3.6 Paso 2.5 — Muestra final y factor de expansión

3.7 Paso 2.6 — Estimaciones del diseño estratificado

4 PARTE 3: Diseño Bietápico con Estratificación Urbano/Rural

4.1 ¿Qué es el diseño bietápico y cuándo se usa?

4.2 Paso 3.1 — Marco de municipios (Primera Etapa)

4.3 Paso 3.2 — Selección de municipios (f₁ = 30%)

4.4 Paso 3.3 — Selección de manzanas (Segunda Etapa, f₂ = 10%)

4.5 Paso 3.4 — Factores de expansión bietápicos

5 PARTE 4: Simulación Monte Carlo para Tamaño Óptimo

5.1 ¿Qué es y para qué sirve la simulación Monte Carlo?

5.2 Paso 4.1 — Función de simulación bietápica

5.3 Paso 4.2 — Ejecutar la simulación

5.4 Interpretación de resultados Monte Carlo

6 PARTE 5: Calibración de Factores de Expansión

6.1 ¿Por qué calibrar?

6.2 Paso 5.1 — Totales poblacionales de referencia

6.3 Paso 5.2 — Calcular y aplicar factores de calibración

6.4 Paso 5.3 — Estimaciones calibradas y comparación

7 Conclusiones Generales