🍲 Cazuela

El Ajiaco2 aumentado con tendencias y efectos de casas encuestadoras

Autor/a

Afiliación

Recetas Electorales

Análisis independiente

“Don’t ask why you should use a multilevel model. Ask instead, Why not?.”
–Richard McElreath

Introduciendo la Cazuela

La Cazuela es una extensión del Ajiaco2 que incorpora covariables al modelo Dirichlet-Multinomial, de manera similar a un modelo lineal generalizado (GLM). Mientras que el Ajiaco2 estima una sola proporción pooled \(\mathbf{p}\) común a todas las encuestas, la Cazuela permite que las proporciones varíen en función de:

1. Tendencia temporal: ¿Cómo ha cambiado el apoyo a cada candidato a lo largo del tiempo?
2. Efecto casa encuestadora: ¿Cada encuestadora tiene un sesgo sistemático hacia o en contra de ciertos candidatos?

El resultado es un modelo que no solo estima el estado actual de las preferencias electorales sino que también captura la dinámica temporal y corrige por las diferencias sistemáticas entre encuestadoras.

Veamos los resultados de la Cazuela Figura 1:

# Extraer p[j,k] para la encuesta más reciente ####
J <- cazuela_data$J
K <- cazuela_data$K
cand_names <- colnames(conteos_2026)

cols_last <- paste0("p[", J, ",", 1:K, "]")

draws_last <- cazuela_fit |>
  posterior::as_draws_df() |>
  tibble::as_tibble() |>
  dplyr::select(dplyr::all_of(cols_last)) |>
  dplyr::rename_with(~ paste0("p_", cand_names), .cols = dplyr::everything()) |>
  tidyr::pivot_longer(cols = dplyr::everything(),
                      names_to = "candidato", values_to = "prop") |>
  dplyr::mutate(candidato = stringr::str_sub(candidato, start = 3L)) |>
  dplyr::group_by(candidato) |>
  dplyr::mutate(candidato_m = mean(prop), .group_by = "drop") |>
  dplyr::left_join(
    candidatos_2026 |> dplyr::mutate(candidato = paste0("cand_", cod)),
    by = "candidato"
  ) |>
  dplyr::filter(!stringr::str_detect(candidato, "ruido"))

nombre_ord_caz <- draws_last |>
  dplyr::distinct(nombre, candidato_m) |>
  dplyr::group_by(nombre) |>
  dplyr::summarise(m = mean(candidato_m), .groups = "drop") |>
  dplyr::arrange(m) |>
  dplyr::pull(nombre)

resultados_caz_todos <- tibble::tribble(
  ~nombre,                       ~color_cand,             ~resultado,
  "Iván Cepeda",                 colores_2026[["ic"]],    0.4092,
  "Abelardo de la Espriella",    colores_2026[["adle"]],  0.4372,
  "Paloma Valencia",              colores_2026[["pv"]],    0.0692,
  "Sergio Fajardo",               colores_2026[["sf"]],    0.0426,
  "Claudia López",                colores_2026[["cl"]],    0.0095
) |>
  dplyr::mutate(
    nombre = factor(nombre, levels = nombre_ord_caz),
    y_num  = as.numeric(factor(nombre, levels = nombre_ord_caz))
  )

draws_last |>
  dplyr::mutate(nombre = factor(nombre, levels = nombre_ord_caz)) |>
  ggplot2::ggplot(ggplot2::aes(x = prop, y = nombre)) +
  ggdist::stat_dist_slabinterval(
    ggplot2::aes(
      fill = color_cand,
      fill_ramp = ggplot2::after_stat(level)
    ),
    .width = c(0.5, 0.8, 0.95, 0.99),
    interval_alpha = 0.95,
    show.legend = c(fill = FALSE, fill_ramp = TRUE, size = FALSE)
  ) +
  ggdist::stat_pointinterval(
    ggplot2::aes(
      label = scales::percent(candidato_m, accuracy = 0.1),
      color = oscurecer(color_cand)
    ),
    geom = "label",
    .width = 0,
    vjust = 1.5,
    fontface = "bold",
    fill = "white",
    label.size = 0.4
  ) +
  ggplot2::geom_segment(
    data = resultados_caz_todos,
    ggplot2::aes(
      x = resultado, xend = resultado,
      y = y_num - 0.4, yend = y_num + 0.4,
      color = color_cand
    ),
    linewidth = 2.5, inherit.aes = FALSE
  ) +
  ggplot2::geom_label(
    data = resultados_caz_todos,
    ggplot2::aes(
      x     = resultado,
      y     = y_num + 0.6,
      label = scales::percent(resultado, accuracy = 0.1),
      fill  = color_cand
    ),
    color         = "white",
    size          = 5,
    fontface      = "bold",
    label.size    = 0,
    label.padding = ggplot2::unit(0.25, "lines"),
    inherit.aes   = FALSE
  ) +
  ggplot2::scale_color_identity(guide = "none") +
  ggplot2::scale_x_continuous(
    labels = scales::percent_format(accuracy = 1),
    expand = ggplot2::expansion(mult = c(0.05, 0.18))
  ) +
  ggplot2::scale_y_discrete(
    labels = scales::label_wrap(20),
    expand = ggplot2::expansion(add = c(0.85, 0.5))   # más aire abajo para el label de la media
  ) +
  ggplot2::scale_fill_identity(na.translate = FALSE) +
  ggdist::scale_fill_ramp_discrete(
    name   = "Intervalo creíble",
    range  = c(0.25, 1),
    breaks = c(0.5, 0.8, .95, 0.99),
    labels = c("50%", "80%", "95%", "99%"),
    na.translate = FALSE
  ) +
  ggplot2::labs(
    title    = "Cazuela: Tendencias y efectos casa encuestadora",
    subtitle = paste0("Modelo probabilístico estimado ", Sys.Date(), " con ", dplyr::n_distinct(encuestas_2026$encuesta_id), " encuestas"),
    caption  = "Línea vertical coloreada ◆: resultado real 31 mayo. Fuente: recetas-electorales.com",
    x        = "Distribución posterior estimada, % intención de voto",
    y        = NULL
  ) +
  ggplot2::theme_minimal(base_size = 18) +
  ggplot2::theme(
    text = ggplot2::element_text(family = "news-cycle"),
    panel.grid.major.y = ggplot2::element_blank(),
    panel.grid.minor   = ggplot2::element_blank(),
    legend.position    = "none"
  )

1

Recuperamos las dimensiones del modelo: número de encuestas J, de candidatos K y los nombres de columna.

2

Construimos los nombres de los parámetros p[J,k] correspondientes a la encuesta más reciente (J).

3

Extraemos esas columnas del posterior, las pasamos a formato largo y calculamos la media posterior de cada candidato, uniendo nombre y color.

4

Ordenamos los candidatos por su media posterior.

5

Definimos los resultados reales de la primera vuelta para superponerlos.

6

Construimos el gráfico base (proporción posterior en X, candidato en Y).

7

Capa principal: distribuciones posteriores como slabs con intervalos creíbles (50/80/95/99%).

8

Etiqueta con la media posterior de cada candidato.

9

Marcador vertical con el resultado real de cada candidato.

10

Etiqueta con el porcentaje real obtenido.

11

Escalas: colores de marca, eje X en porcentaje, nombres ajustados y rampa de intervalos creíbles.

12

Títulos, subtítulo (fecha y número de encuestas) y nota de fuente.

13

Tema minimalista con tipografía de marca y sin leyenda.

Figura 1: Cazuela: Dirichlet-Multinomial con covariables

Valencia o De la Espriella? Competencia probabilística

Una ventaja del enfoque bayesiano es que las comparaciones directas entre candidatos son triviales: basta con calcular la fracción de draws en los que la proporción de un candidato supera la del otro.

Según la Cazuela a 24 de mayo de 2026, la probabilidad posterior de que Paloma Valencia obtenga más votos que De la Espriella en primera vuelta es 1.0%.

library(ggtext)

draws_comp |>
  tidyr::pivot_longer(
    cols = c(p_cand_pv, p_cand_adle),
    names_to = "candidato",
    values_to = "prop"
  ) |>
  dplyr::mutate(
    color_cand = dplyr::case_when(
      candidato == "p_cand_pv"   ~ "#63B9E9",
      candidato == "p_cand_adle" ~ "#000066"
    )
  ) |>
  ggplot2::ggplot(ggplot2::aes(x = prop, fill = color_cand)) +
  ggplot2::geom_density(alpha = 0.55, color = NA) +
  ggplot2::scale_x_continuous(labels = scales::percent_format(accuracy = 1)) +
  ggplot2::scale_fill_identity(
    guide  = "legend",
    labels = c("#63B9E9" = "Paloma Valencia", "#000066" = "De la Espriella"),
    breaks = c("#63B9E9", "#000066")
  ) +
  ggtext::geom_richtext(
    data = tibble::tibble(x = Inf, y = Inf),
    ggplot2::aes(
      x = x, y = y,
      label = sprintf(
        "<b>P(Paloma > De la Espriella) = <span style='color:#63B9E9'>%.1f%%</span></b>",
        prob_pv_gt_adle * 100
      )
    ),
    inherit.aes  = FALSE,
    hjust        = 1.1,
    vjust        = 1.5,
    size         = 7,
    family       = "news-cycle",
    fill         = "#F8F9FA",
    color        = "#212529",
    label.color  = NA,
    label.padding = ggplot2::unit(c(0.4, 0.6, 0.4, 0.6), "lines")
  ) +
  ggplot2::labs(
    title    = "Paloma Valencia vs. De la Espriella: distribuciones posteriores",
    subtitle = "Cazuela — encuesta más reciente, primera vuelta 2026",
    x        = "Distribución posterior estimada, % intención de voto",
    y        = "Densidad",
    fill     = NULL,
    caption  = "Fuente: https://recetas-electorales.netlify.app/"
  ) +
  ggplot2::theme_minimal(base_size = 24) +
  ggplot2::theme(
    text               = ggplot2::element_text(family = "news-cycle"),
    panel.grid.minor   = ggplot2::element_blank(),
    legend.position    = "top"
  )

1

Cargamos ggtext para rotular con texto enriquecido (HTML) dentro del gráfico.

2

Tomamos las muestras posteriores de Paloma Valencia y De la Espriella (encuesta más reciente), las pasamos a formato largo y les asignamos su color de marca.

3

Construimos el gráfico base y dibujamos la densidad posterior de cada candidato.

4

Eje X en porcentaje y leyenda con los nombres asociados a cada color.

5

Anotación enriquecida con la probabilidad posterior P(Paloma > De la Espriella).

6

Títulos, ejes y nota de fuente.

7

Tema minimalista con tipografía de marca y leyenda arriba.

Figura 2: Distribuciones posteriores de Paloma Valencia y De la Espriella (encuesta más reciente).

Tendencia temporal

La gran ventaja de la Cazuela sobre el Ajiaco2 es que podemos visualizar cómo han cambiado las proporciones estimadas a lo largo del tiempo, corrigiendo por el efecto casa de cada encuestadora Figura 3.

# Extraer p[j,k] para todas las encuestas ####
draws_all <- cazuela_fit |>
  posterior::as_draws_df() |>
  tibble::as_tibble()

# Construir data.frame largo con p[j,k] para cada draw y encuesta
p_draws <- purrr::map_dfr(1:J, function(j) {
  # Seleccionar columnas exactas p[j,1], p[j,2], ..., p[j,K]
  cols_j <- paste0("p[", j, ",", 1:K, "]")
  draws_all |>
    dplyr::select(dplyr::all_of(cols_j)) |>
    dplyr::rename_with(~ cand_names, .cols = dplyr::everything()) |>
    dplyr::mutate(.draw = dplyr::row_number(), .stan_idx = j)
}) |>
  dplyr::left_join(
    encuestas_meta |> dplyr::mutate(.stan_idx = dplyr::row_number()),
    by = ".stan_idx"
  ) |>
  dplyr::select(-".stan_idx") |>
  tidyr::pivot_longer(
    cols = dplyr::all_of(cand_names),
    names_to = "candidato",
    values_to = "prop"
  ) |>
  dplyr::mutate(cod = stringr::str_sub(candidato, start = 6L)) |>
  dplyr::left_join(candidatos_2026, by = c("cod")) |>
  dplyr::filter(cod != "ruido")

# Resumir por fecha y candidato
p_summary <- p_draws |>
  dplyr::group_by(fecha, nombre, color_cand, cod) |>
  dplyr::summarise(
    media = mean(prop),
    q025 = quantile(prop, 0.025),
    q975 = quantile(prop, 0.975),
    q25 = quantile(prop, 0.25),
    q75 = quantile(prop, 0.75),
    .groups = "drop"
  )

# Orden global de candidatos por media posterior (usado en todas las gráficas)
orden_candidatos <- p_summary |>
  dplyr::group_by(nombre) |>
  dplyr::summarise(media_global = mean(media), .groups = "drop") |>
  dplyr::arrange(dplyr::desc(media_global)) |>
  dplyr::pull(nombre)

resultados_caz <- tibble::tribble(
  ~nombre,                       ~color_cand,             ~resultado,
  "Iván Cepeda",                 colores_2026[["ic"]],    0.4092,
  "Abelardo de la Espriella",    colores_2026[["adle"]],  0.4372,
  "Paloma Valencia",              colores_2026[["pv"]],    0.0692,
  "Sergio Fajardo",               colores_2026[["sf"]],    0.0426,
  "Claudia López",                colores_2026[["cl"]],    0.0095
) |>
  dplyr::mutate(nombre = factor(nombre, levels = orden_candidatos))

p_summary |>
  dplyr::mutate(nombre = factor(nombre, levels = orden_candidatos)) |>
  ggplot2::ggplot(ggplot2::aes(x = fecha, y = media, color = color_cand, fill = color_cand)) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = q025, ymax = q975), alpha = 0.15, color = NA) +
  ggplot2::geom_ribbon(ggplot2::aes(ymin = q25, ymax = q75), alpha = 0.3, color = NA) +
  ggplot2::geom_line(linewidth = 1.2) +
  ggplot2::geom_point(size = 2.5) +
  ggplot2::geom_hline(
    data = resultados_caz,
    ggplot2::aes(yintercept = resultado),
    linetype = "dashed", color = "#333333", linewidth = 0.8,
    inherit.aes = FALSE
  ) +
  ggplot2::geom_label(
    data = resultados_caz,
    ggplot2::aes(
      x     = max(p_summary$fecha),
      y     = resultado,
      label = scales::percent(resultado, accuracy = 0.1),
      fill  = color_cand
    ),
    hjust         = 1.1,
    color         = "white",
    size          = 3.5,
    fontface      = "bold",
    label.size    = 0,
    label.padding = ggplot2::unit(0.2, "lines"),
    inherit.aes   = FALSE
  ) +
  ggplot2::scale_color_identity() +
  ggplot2::scale_fill_identity() +
  ggplot2::scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  ggplot2::facet_wrap(~ nombre, scales = "free_y", ncol = 3) +
  ggplot2::labs(
    title = "Cazuela: Tendencia temporal por candidato",
    subtitle = paste0("Intervalos creíbles del 50% y 95%, estimados con ", dplyr::n_distinct(encuestas_2026$encuesta_id), " encuestas"),
    caption = "Fuente: https://recetas-electorales.netlify.app/",
    x = NULL,
    y = "% intención de voto"
  ) +
  ggplot2::theme_minimal(base_size = 14) +
  ggplot2::theme(
    text = ggplot2::element_text(family = "news-cycle"),
    panel.grid.minor = ggplot2::element_blank(),
    strip.text = ggplot2::element_text(face = "bold", size = 14),
    legend.position = "none"
  )

1

Extraemos todas las muestras posteriores (draws) del modelo.

2

Para cada encuesta j, seleccionamos sus columnas p[j,k], las unimos con los metadatos (fecha, encuestadora) y pasamos a formato largo por candidato.

3

Resumimos por fecha y candidato: media posterior e intervalos creíbles del 50% y 95%.

4

Definimos un orden global de candidatos por su media posterior (común a todos los gráficos).

5

Definimos los resultados reales de la primera vuelta para superponerlos.

6

Construimos el gráfico base (fecha en X, intención media en Y).

7

Bandas de incertidumbre: intervalo creíble del 95% (claro) y del 50% (oscuro), más línea y puntos de la media.

8

Línea horizontal punteada con el resultado real de cada candidato y su etiqueta.

9

Un panel (facet) por candidato con eje Y libre.

10

Títulos, subtítulo (número de encuestas) y nota de fuente.

11

Tema minimalista con tipografía de marca y sin leyenda.

Figura 3: Tendencia temporal estimada por la Cazuela

Efectos de casa

Cada encuestadora tiene un sesgo sistemático implícito. La Cazuela modela estos efectos de casa como un modelo multinivel de efectos aleatorios: cada vector \(\boldsymbol{\delta}_h \in \mathbb{R}^{K-1}\) se extrae de una distribución normal multivariada con matriz de covarianza estimada vía descomposición LKJ. Esto permite que el modelo aprenda si una encuestadora que sobreestima a un candidato tiende a subestimar a otro Figura 4.

# Extraer delta[k, h] (nota: índices [k, h] en el modelo actualizado) ####
K_minus_1 <- K - 1
H <- cazuela_data$H

delta_draws <- draws_all |>
  dplyr::select(dplyr::matches("^delta\\[")) |>
  tidyr::pivot_longer(dplyr::everything(), names_to = "param", values_to = "value") |>
  dplyr::mutate(
    indices = stringr::str_extract(param, "\\[(.+)\\]", group = 1),
    k = as.integer(stringr::str_extract(indices, "^\\d+")),
    h = as.integer(stringr::str_extract(indices, "\\d+$")),
    encuestadora = house_levels[h],
    candidato = cand_names[k],
    cod = stringr::str_sub(candidato, start = 6L)
  ) |>
  dplyr::left_join(candidatos_2026, by = c("cod")) |>
  dplyr::filter(cod != "ruido")

delta_summary <- delta_draws |>
  dplyr::group_by(encuestadora, nombre, color_cand, cod) |>
  dplyr::summarise(
    media = mean(value),
    q025 = quantile(value, 0.025),
    q975 = quantile(value, 0.975),
    .groups = "drop"
  )

delta_summary |>
  dplyr::mutate(nombre = factor(nombre, levels = orden_candidatos)) |>
  ggplot2::ggplot(ggplot2::aes(x = media, y = encuestadora, color = color_cand)) +
  ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  ggplot2::geom_pointrange(ggplot2::aes(xmin = q025, xmax = q975), size = 0.5) +
  ggplot2::scale_color_identity() +
  ggplot2::facet_wrap(~ nombre, scales = "free_x", ncol = 2) +
  ggplot2::labs(
    title = "Cazuela: Efectos de casa por encuestadora",
    subtitle = "Desviaciones multinivel (log-ratio), intervalos creíbles del 95%",
    caption = "Fuente: https://recetas-electorales.netlify.app/",
    x = "Efecto casa (log-ratio escala)",
    y = NULL
  ) +
  ggplot2::theme_minimal(base_size = 16) +
  ggplot2::theme(
    text = ggplot2::element_text(family = "news-cycle"),
    panel.grid.minor = ggplot2::element_blank(),
    strip.text = ggplot2::element_text(face = "bold", size = 14),
    legend.position = "none"
  )

1

Recuperamos las dimensiones: K-1 log-ratios y H encuestadoras.

2

Extraemos los parámetros delta[k,h] del posterior, decodificamos los índices candidato/encuestadora y les unimos nombre y color.

3

Resumimos cada efecto de casa: media posterior e intervalo creíble del 95%.

4

Construimos el gráfico (efecto de casa en X, encuestadora en Y).

5

Línea de referencia en cero (sin sesgo) y rangos posteriores con geom_pointrange.

6

Un panel (facet) por candidato con eje X libre.

7

Títulos, ejes y nota de fuente.

8

Tema minimalista con tipografía de marca y sin leyenda.

Figura 4: Efectos de casa estimados por la Cazuela

El partial pooling en acción

La característica más poderosa del modelo multinivel es el shrinkage: cada encuestadora es un punto en el espacio de estimaciones y el modelo las “jala” hacia el centro según cuánta evidencia tiene. La animación Figura 5 muestra cómo se mueve cada casa encuestadora entre tres escenarios:

1. Pooling completo: Todas colapsan al mismo punto, el promedio global.
2. Sin pooling: Cada encuestadora se estimada de forma independiente.
3. Partial pooling (Cazuela): El modelo jerárquico regula —las encuestadoras con pocas encuestas se encogen más hacia el centro que las que tienen varias muestras.

Figura 5: Pooling en la Cazuela

Correlación entre efectos de casa

El modelo estima explícitamente la matriz de correlación \(\Omega\) entre los efectos de casa de los candidatos usando un prior LKJ. Esto captura si las encuestadoras que sobreestiman a un candidato tienden a subestimar a otro Figura 6.

# Extraer Omega[i,j] ####
omega_draws <- draws_all |>
  dplyr::select(dplyr::matches("^Omega\\[")) |>
  tidyr::pivot_longer(dplyr::everything(), names_to = "param", values_to = "value") |>
  dplyr::mutate(
    indices = stringr::str_extract(param, "\\[(.+)\\]", group = 1),
    i = as.integer(stringr::str_extract(indices, "^\\d+")),
    j = as.integer(stringr::str_extract(indices, "\\d+$"))
  ) |>
  dplyr::group_by(i, j) |>
  dplyr::summarise(media = mean(value), .groups = "drop") |>
  dplyr::mutate(
    cand_i = cand_names[i],
    cand_j = cand_names[j],
    cod_i = stringr::str_sub(cand_i, start = 6L),
    cod_j = stringr::str_sub(cand_j, start = 6L)
  ) |>
  dplyr::left_join(candidatos_2026 |> dplyr::select(cod, nombre), by = c("cod_i" = "cod")) |>
  dplyr::rename(nombre_i = nombre) |>
  dplyr::left_join(candidatos_2026 |> dplyr::select(cod, nombre), by = c("cod_j" = "cod")) |>
  dplyr::rename(nombre_j = nombre) |>
  dplyr::filter(!is.na(nombre_i), !is.na(nombre_j))

orden_corr <- orden_candidatos[orden_candidatos %in% unique(omega_draws$nombre_i)]

omega_draws |>
  dplyr::mutate(
    nombre_i = factor(nombre_i, levels = orden_corr),
    nombre_j = factor(nombre_j, levels = orden_corr)
  ) |>
  ggplot2::ggplot(ggplot2::aes(x = nombre_i, y = nombre_j, fill = media)) +
  ggplot2::geom_tile() +
  ggplot2::geom_text(ggplot2::aes(label = round(media, 2)), size = 3.5) +
  ggplot2::scale_fill_gradient2(low = "#2166AC", mid = "white", high = "#B2182B",
                       midpoint = 0, limits = c(-1, 1),
                       name = "Correlación") +
  ggplot2::scale_x_discrete(labels = scales::label_wrap(12)) +
  ggplot2::scale_y_discrete(labels = scales::label_wrap(12)) +
  ggplot2::labs(
    title = "Correlación entre efectos de casa",
    subtitle = "Matriz posterior media estimada con prior LKJ(2)",
    caption = "Fuente: https://recetas-electorales.netlify.app/",
    x = NULL, y = NULL
  ) +
  ggplot2::theme_minimal(base_size = 18) +
  ggplot2::theme(
    text = ggplot2::element_text(family = "news-cycle"),
    axis.text.x = ggplot2::element_text(angle = 45, hjust = 1),
    panel.grid = ggplot2::element_blank(),
    legend.position = "right"
  )

1

Extraemos los parámetros Omega[i,j] (matriz de correlación) del posterior, decodificamos los índices y promediamos cada celda, uniendo los nombres de candidato.

2

Reordenamos las filas/columnas según el orden global de candidatos.

3

Construimos el mapa de calor (candidato × candidato) con la correlación como relleno y su valor como texto.

4

Escala divergente azul-blanco-rojo centrada en cero, de −1 a 1.

5

Títulos y nota de fuente.

6

Tema minimalista con tipografía de marca y etiquetas del eje X rotadas.

Figura 6: Matriz de correlación posterior de efectos de casa

Diagnósticos

El muestreo MCMC:

bayesplot::mcmc_trace(cazuela_fit$draws(), pars = ggplot2::vars(matches("beta0_raw")))

1

Traceplots de las cadenas MCMC para los interceptos beta0_raw: si las cadenas se mezclan bien, el muestreo convergió.

Receta Cazuela

La Cazuela extiende el Ajiaco2 al permitir que las proporciones de candidatos varíen por encuesta como función de covariables. En lugar de un solo vector \(\mathbf{p}\) pooled, cada encuesta \(j\) tiene su propio \(\mathbf{p}_j\) construido a partir de un predictor lineal con función de enlace softmax, como en un GLM multinomial.

La novedad clave es que los efectos de casa se modelan como efectos aleatorios multinivel con una estructura de covarianza completa. Los vectores \(\boldsymbol{\delta}_h\) se extraen de una distribución normal multivariada cuya matriz de covarianza se descompone en desviaciones estándar por candidato (\(\boldsymbol{\sigma}_\delta\)) y una matriz de correlación (\(\Omega\)) con prior LKJ. Se usa la parametrización no-centrada (\(\boldsymbol{\delta}_h = \mathrm{diag}(\boldsymbol{\sigma}_\delta) \cdot L_\Omega \cdot \mathbf{z}_h\)) para mejorar la eficiencia del muestreo MCMC:

\[ \begin{aligned} &\textbf{Predictor lineal:}\\ &\eta_{jk} = \beta_{0k} + \beta_{1k} \cdot t_j + \delta_{k,h_j}, \quad k=1,\dots,K-1,\\ &\eta_{jK} = 0 \quad \text{(candidato de referencia)}.\\[6pt] &\textbf{Enlace:}\\ &\mathbf{p}_j = \mathrm{softmax}(\boldsymbol{\eta}_j).\\[6pt] &\textbf{Efectos aleatorios de casa (multinivel):}\\ &\boldsymbol{\delta}_h \sim \mathrm{Normal}_{K-1}(\mathbf{0},\; \Sigma), \quad h=1,\dots,H,\\ &\Sigma = \mathrm{diag}(\boldsymbol{\sigma}_\delta)\;\Omega\;\mathrm{diag}(\boldsymbol{\sigma}_\delta).\\[6pt] &\textbf{Parametrización no-centrada:}\\ &\boldsymbol{\delta}_h = \mathrm{diag}(\boldsymbol{\sigma}_\delta) \cdot L_\Omega \cdot \mathbf{z}_h, \quad \mathbf{z}_h \sim \mathrm{Normal}_{K-1}(\mathbf{0}, I).\\[6pt] &\textbf{Priors:}\\ &\beta_{0k} \sim \mathrm{Normal}(0, 2),\\ &\beta_{1k} \sim \mathrm{Normal}(0, 1),\\ &\sigma_{\delta,k} \sim \mathrm{Normal}^+(0, 1),\\ &L_\Omega \sim \mathrm{LKJ\_cholesky}(2),\\ &\kappa \sim \mathrm{LogNormal}(\log 200, 0.7).\\[6pt] &\textbf{Likelihood:}\\ &\mathbf{Y}_j \mid \mathbf{p}_j, \kappa \sim \mathrm{Dirichlet\text{-}Multinomial}(n_j, \kappa \mathbf{p}_j). \end{aligned} \]

Donde \(t_j\) es el tiempo centrado y escalado de la encuesta \(j\), \(h_j\) es el índice de la encuestadora, y \(L_\Omega\) es el factor Cholesky de la matriz de correlación \(\Omega\). El prior \(\mathrm{LKJ}(2)\) favorece correlaciones moderadas, lo que actúa como regularización contra matrices de correlación extremas.

Modelo

Stan

data {
  int<lower=1> J;                     // encuestas
  int<lower=2> K;                     // candidatos
  int<lower=1> H;                     // encuestadoras
  array[J] int<lower=0> n;           // tamaño muestral
  array[J, K] int<lower=0> Y;        // conteos Y[j, k]
  vector[J] time;                     // covariable temporal (centrada)
  array[J] int<lower=1, upper=H> house; // índice de encuestadora
}

parameters {
  vector[K - 1] beta0_raw;           // intercepto (log-ratio vs candidato K)
  vector[K - 1] beta_time_raw;       // pendiente temporal (log-ratio vs K)

  // Efectos aleatorios de casa: multinivel con covarianza LKJ
  matrix[K - 1, H] z_delta;          // efectos estandarizados (non-centered)
  vector<lower=0>[K - 1] sigma_delta; // sd por candidato
  cholesky_factor_corr[K - 1] L_Omega; // factor Cholesky de correlación

  real<lower=0> kappa;                // concentración DM
}

transformed parameters {
  // Reconstruir efectos de casa: delta[,h] = diag(sigma) * L * z[,h]
  matrix[K - 1, H] delta;
  delta = diag_pre_multiply(sigma_delta, L_Omega) * z_delta;

  matrix[J, K] eta;
  for (j in 1:J) {
    for (k in 1:(K - 1)) {
      eta[j, k] = beta0_raw[k] + beta_time_raw[k] * time[j] + delta[k, house[j]];
    }
    eta[j, K] = 0;
  }

  array[J] simplex[K] p;
  for (j in 1:J) {
    p[j] = softmax(eta[j]');
  }
}

model {
  // Priors efectos fijos
  beta0_raw ~ normal(0, 2);
  beta_time_raw ~ normal(0, 1);

  // Priors multinivel para efectos de casa
  sigma_delta ~ normal(0, 1);
  L_Omega ~ lkj_corr_cholesky(2);
  to_vector(z_delta) ~ std_normal();

  // Prior concentración
  kappa ~ lognormal(log(200), 0.7);

  // Likelihood
  for (j in 1:J) {
    Y[j] ~ dirichlet_multinomial(kappa * to_vector(p[j]));
  }
}

generated quantities {
  array[J, K] int Y_rep;
  for (j in 1:J) {
    Y_rep[j] = dirichlet_multinomial_rng(kappa * to_vector(p[j]), n[j]);
  }

  // Matriz de correlación de efectos de casa
  corr_matrix[K - 1] Omega;
  Omega = multiply_lower_tri_self_transpose(L_Omega);
}

Comparación con el Ajiaco2

	Ajiaco2	Cazuela
Proporciones	Un solo \(\mathbf{p}\) pooled	\(\mathbf{p}_j\) por encuesta
Tendencia temporal	No	Sí (\(\beta_{1k}\))
Efecto casa	No	Sí, multinivel con LKJ
Covarianza casa	No	\(\Sigma = \mathrm{diag}(\sigma)\,\Omega\,\mathrm{diag}(\sigma)\)
Parametrización	Directa	Non-centered (\(\mathbf{z}_h\))
Enlace	Dirichlet prior directo	Softmax (multinomial logit)
Parámetros	\(K + 1\)	\((K-1)(2 + H) + K + \binom{K-1}{2}\)
Interpretación	Estado promedio	Estado dinámico corregido

cazuela_fit$summary(variables = c("beta0_raw", "beta_time_raw", "kappa"))

1

Tabla resumen posterior de los efectos fijos (beta0_raw, beta_time_raw) y la concentración kappa, con intervalos y diagnósticos de convergencia (rhat, ess).

# A tibble: 11 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 beta0_raw…   1.12     1.12    0.144   0.132   8.82e-1  1.35e+0  1.00    4676.
 2 beta0_raw…   0.596    0.596   0.173   0.165   3.15e-1  8.79e-1  1.00    4863.
 3 beta0_raw…   0.510    0.507   0.215   0.198   1.71e-1  8.56e-1  1.00    5055.
 4 beta0_raw…  -1.41    -1.41    0.210   0.197  -1.75e+0 -1.07e+0  1.00    6432.
 5 beta0_raw…  -1.70    -1.71    0.191   0.174  -2.01e+0 -1.39e+0  1.00    8567.
 6 beta_time…   0.0111   0.0115  0.0512  0.0512 -7.24e-2  9.43e-2  1.00    8041.
 7 beta_time…   0.0867   0.0863  0.0548  0.0537 -2.79e-3  1.77e-1  1.00    8675.
 8 beta_time…  -0.0867  -0.0871  0.0571  0.0563 -1.79e-1  8.39e-3  1.00    8655.
 9 beta_time…  -0.204   -0.205   0.0911  0.0903 -3.51e-1 -5.18e-2  1.00   12919.
10 beta_time…  -0.0911  -0.0912  0.104   0.101  -2.60e-1  8.21e-2  1.00   14577.
11 kappa      310.     305.     63.8    62.2     2.14e+2  4.24e+2  1.00    6287.
# ℹ 1 more variable: ess_tail <dbl>

---

Resultado — primera vuelta 31 de mayo de 2026

NotaResultados oficiales

Candidato	% Votos válidos
Abelardo de la Espriella	43.72%
Iván Cepeda	40.92%
Paloma Valencia	6.92%
Sergio Fajardo	4.26%
Claudia López	0.95%

La Cazuela modelaba tendencias temporales y efectos de casa encuestadora. En las últimas encuestas el coeficiente de tendencia de ADLE era positivo — la señal estaba ahí — pero la extrapolación al 31 de mayo no alcanzó el resultado final de 43.72%. Para la evaluación completa de precisión ver Café.