MRP y raking: postestratificar sin tener la conjunta

2026

estadística

encuestas

MRP

modelos mixtos

Multilevel Regression and Poststratification con la encuesta del CIS, y raking para fabricar la tabla de postestratificación que nadie publica.

Author

José Luis Cañadas Reche

Published

October 4, 2026

Este post está escrito entre Claude y yo. Ha llevado varias horas, porque le he guiado en cada paso y comprobado lo que hace en cada paso. También hay cosas que he reescrito yo, e intentar que fuera pedagógico en la medida de lo posible.
Creo que este post puede ser muy útil a todo aquel que trabaje con encuestas y tenga que hacer calibración y postestratificación, y además combinarlo con un modelo bayesiano, que funciona mejor en combinaciones de variables dónde apenas hay dato ( o incluso sin dato) y tener en cuenta la incertidumbre.

Listening

La idea

MRP son las siglas de Multilevel Regression and Poststratification. Son dos pasos:

Multilevel Regression: un modelo que predice la intención de voto en función de unas cuantas variables (comunidad, sexo, edad, estudios, recuerdo de voto).
Poststratification: aplicar esas predicciones a la población real, con cuánta gente hay de verdad en cada combinación de esas variables.

Uso el barómetro del CIS de septiembre (estudio 3577). La descarga y la preparación de los datos están en scripts aparte, y en el post solo leo lo que generan.

Reproducir el análisis

Los scripts están en el repositorio del blog. Se ejecutan en este orden:

00_partidos.R: no se ejecuta solo. Define los partidos, sus territorios y la recodificación de las etiquetas del CIS, y lo cargan los demás scripts.
01_descarga.R: descarga la población del INE (ECP y EPA) y los resultados del 23-J.
02_harmoniza.R: recodifica la encuesta del CIS y la deja con las mismas categorías que el INE y el Ministerio del Interior.
03_targets.R: construye los objetivos de población para el raking.
04_raking.R: imputa el recuerdo de voto que falta y hace el raking. Genera la tabla de postestratificación.
05_mrp.R: el modelo y la postestratificación del post, en un script, más una comparación con el raking clásico de pesos.

Antes de empezar:

Descarga los microdatos del estudio 3577 desde la web del CIS y guarda 3577.sav en datos/. El script no los descarga.
Ejecuta los scripts desde la raíz del repositorio. Las rutas son relativas a ella.
Paquetes: ineapir, infoelectoral, haven, mipfp, nnet, survey, brms y cmdstanr.

Si no quieres ejecutar los scripts, los ficheros que lee el post, incluidos los modelos ya ajustados, están en esta release. Para bajarlos todos a datos/:

dir_datos <- "2026/10/mrp-raking/datos"
url_base <- "https://github.com/joscani/blog_quarto/releases/download/mrp-raking-datos/"
ficheros <- c(
  "cis_recodificado.rds", "cis_imputado.rds", "targets.rds",
  "tabla_postestratificacion.rds", "mod_sin_rec.rds", "mod_mrp.rds"
)

dir.create(dir_datos, recursive = TRUE, showWarnings = FALSE)
for (f in ficheros) {
  download.file(paste0(url_base, f), file.path(dir_datos, f), mode = "wb")
}

Fuente de los datos de la encuesta: Centro de Investigaciones Sociológicas, estudio 3577, descargado el 22 de septiembre de 2026. Los ficheros de la release son recodificaciones propias. El CIS no participa en este análisis ni lo respalda.

Primero la M: por qué multinivel

Mostrar código

library(tidyverse)
library(brms)
library(nnet)
library(gt)

dir_datos <- here::here("2026/10/mrp-raking/datos")
voto_lv <- c(
  "PP", "PSOE", "VOX", "SUMAR",
  "ERC", "JUNTS", "BILDU", "PNV", "BNG", "CCA",
  "OTROS"
)
voto5_lv <- c("PP", "PSOE", "VOX", "SUMAR", "OTROS") # para la primera parte

# Comunidades en las que se presenta cada partido regional (23-J)
territorio <- list(
  ERC = "Cataluña",
  JUNTS = "Cataluña",
  BILDU = c("País Vasco", "Navarra"),
  PNV = "País Vasco",
  BNG = "Galicia",
  CCA = "Canarias"
)

cis <- readRDS(file.path(dir_datos, "cis_imputado.rds"))

dat <- cis |>
  filter(!is.na(voto), !is.na(edu), !is.na(edad), !is.na(rec))

nrow(dat)
#> [1] 3128

De los 4018 entrevistados, solo hay 3128 que declaran intención de voto válido. Si estuviéramos haciendo la estimación en serio , habría que hacer ciertas hipótesis o ver si esa falta de respuesta es aleatoria o viene determinada por otras variables. No estamos ahí, así que me los cargo.

Viendo algunas filas de los datos.

Mostrar código

dat |>
  select(ccaa, sexo, edad, edu, rec, voto) |>
  head()
#> # A tibble: 6 × 6
#>   ccaa      sexo  edad  edu   rec   voto 
#>   <fct>     <fct> <fct> <fct> <fct> <fct>
#> 1 Andalucía H     45-54 media PSOE  SUMAR
#> 2 Andalucía M     45-54 alta  PSOE  VOX  
#> 3 Andalucía M     65+   media PSOE  PSOE 
#> 4 Andalucía M     25-34 alta  PP    VOX  
#> 5 Andalucía M     65+   media SUMAR PSOE 
#> 6 Andalucía M     45-54 alta  PSOE  PSOE

El recuerdo de voto es lo que más afecta a la intención, pero para ejemplificar mejor lo que hace mejor un modelo multinivel vs uno clásico lo voy a eliminar en esta primera parte del post

Cómo agrupo los partidos

Tanto en la intención de voto (voto) como en el recuerdo (rec) uso los mismos grupos:

PP, PSOE y VOX.
SUMAR, que en realidad es SUMAR y aliados: Sumar, Podemos, IU, Compromís, Más Madrid y CHA. Es decir, los partidos que fueron juntos en la candidatura de Sumar el 23-J. Así el recuerdo y la intención miden lo mismo.
Los nacionalistas y regionalistas, cada uno por separado: ERC y JUNTS (Cataluña), BILDU (País Vasco y Navarra), PNV (País Vasco), BNG (Galicia) y CCA (Coalición Canaria).
OTROS: el resto de partidos (Se Acabó la Fiesta, PACMA, CUP, UPN…) y el voto en blanco.

Los partidos regionales solo se pueden votar en su comunidad. Las pocas respuestas de fuera de su territorio (por ejemplo, alguien en Baleares que dice ERC) las paso a OTROS.

El recuerdo tiene además dos categorías: ABST (se abstuvo) y NO_PODIA (tenía menos de 18 años en julio de 2023).

No se llega a tener datos de todas las combinaciones

Si cruzo comunidad, sexo, edad y estudios salen 19 × 2 × 7 × 3 = 798 combinaciones posibles o celdas, para entendernos. Veamos cuántas personas hay en cada celda con datos.

En esta primera parte agrupo los partidos regionales dentro de OTROS y me quedo con cinco categorías. Para ver qué hace un modelo multinivel no hacen falta más, y los gráficos se leen mejor:

Mostrar código

dat_dem <- dat |>
  mutate(voto = fct_other(voto, keep = voto5_lv[-5], other_level = "OTROS"))

count(dat_dem, voto)
#> # A tibble: 5 × 2
#>   voto      n
#>   <fct> <int>
#> 1 PP      804
#> 2 PSOE    889
#> 3 VOX     452
#> 4 SUMAR   333
#> 5 OTROS   650

Mostrar código

celdas_dem <- dat_dem |>
  count(ccaa, sexo, edad, edu, voto) |>
  pivot_wider(names_from = voto, values_from = n, values_fill = 0)

celdas_dem$y <- as.matrix(celdas_dem[, voto5_lv])
celdas_dem$n <- rowSums(celdas_dem$y)

nrow(celdas_dem)
#> [1] 544

table(cut(celdas_dem$n, c(0, 1, 2, 3, 5, 10, Inf)))
#> 
#>    (0,1]    (1,2]    (2,3]    (3,5]   (5,10] (10,Inf] 
#>      134       96       70       75       94       75

Solo 544 de las 798 celdas tienen alguna entrevista. De esas, 230 tienen una o dos personas. Estimar la intención de voto celda a celda, con su proporción bruta, no tiene sentido. Hay que ajustar un modelo que use información de unas celdas para otras.

Por comunidad pasa lo mismo

Mostrar código

n_ccaa <- dat |> count(ccaa, name = "n_ccaa") |> arrange(n_ccaa)
n_ccaa
#> # A tibble: 19 × 2
#>    ccaa                 n_ccaa
#>    <fct>                 <int>
#>  1 Ceuta                    16
#>  2 Melilla                  17
#>  3 Navarra                  61
#>  4 Cantabria                71
#>  5 Extremadura              73
#>  6 La Rioja                 75
#>  7 Aragón                   77
#>  8 Asturias                 78
#>  9 Baleares                 83
#> 10 Murcia                   88
#> 11 País Vasco              133
#> 12 Castilla-La Mancha      134
#> 13 Canarias                148
#> 14 Castilla y León         156
#> 15 Galicia                 176
#> 16 Comunitat Valenciana    300
#> 17 Madrid                  450
#> 18 Cataluña                458
#> 19 Andalucía               534

Ceuta y Melilla tienen 16 y 17 entrevistas. Andalucía, 534.

Dos modelos

Do modelos multinomiales con las mismas variables. Lo único que cambia es cómo tratan las variables categóricas.

Efectos fijos (nnet::multinom): cada comunidad tiene su coeficiente y se estima solo con sus datos. Es lo que se conoce como no pooling.
Multinivel (brms): los coeficientes de comunidad, edad y estudios salen de una distribución normal común, con su desviación típica estimada. Es el partial pooling.

Mostrar código

mod_fe <- multinom(
  voto ~ sexo + edad + edu + ccaa,
  data = dat_dem,
  trace = FALSE,
  maxit = 500
)

El modelo multinivel se ajusta sobre los datos agregados por celda, con trials(n). Da lo mismo que ajustarlo sobre cada persona y es mucho más rápido. Ver post antiguo

Mostrar código

# Las mismas priors para cada partido (cada uno es un "mu" respecto al PP)
priors_multinom <- function(partidos) {
  Reduce(
    `+`,
    lapply(paste0("mu", partidos[-1]), function(dp) {
      prior_string("normal(0, 1.5)", class = "Intercept", dpar = dp) +
        prior_string("normal(0, 1)", class = "b", dpar = dp) +
        prior_string("exponential(1)", class = "sd", dpar = dp)
    })
  )
}

priors_dem <- priors_multinom(voto5_lv)
priors_dem
#>           prior     class coef group resp    dpar nlpar   lb   ub tag source
#>  normal(0, 1.5) Intercept                  muPSOE       <NA> <NA>       user
#>    normal(0, 1)         b                  muPSOE       <NA> <NA>       user
#>  exponential(1)        sd                  muPSOE       <NA> <NA>       user
#>  normal(0, 1.5) Intercept                   muVOX       <NA> <NA>       user
#>    normal(0, 1)         b                   muVOX       <NA> <NA>       user
#>  exponential(1)        sd                   muVOX       <NA> <NA>       user
#>  normal(0, 1.5) Intercept                 muSUMAR       <NA> <NA>       user
#>    normal(0, 1)         b                 muSUMAR       <NA> <NA>       user
#>  exponential(1)        sd                 muSUMAR       <NA> <NA>       user
#>  normal(0, 1.5) Intercept                 muOTROS       <NA> <NA>       user
#>    normal(0, 1)         b                 muOTROS       <NA> <NA>       user
#>  exponential(1)        sd                 muOTROS       <NA> <NA>       user

Mostrar código


mod_dem <- brm(
  y | trials(n) ~ sexo + (1 | edad) + (1 | edu) + (1 | ccaa),
  data = celdas_dem,
  family = multinomial(),
  prior = priors_dem,
  chains = 4,
  cores = 4,
  iter = 2000,
  refresh = 0,
  silent = 2,
  control = list(adapt_delta = 0.95),
  backend = "cmdstanr",
  seed = 2026,
  file = file.path(dir_datos, "mod_sin_rec")
)

La categoría de referencia es el PP. Así que cada partido tiene su propio conjunto de coeficientes, en escala logit respecto al PP.

El efecto de cada comunidad

Veamos el coeficiente de cada comunidad en los dos modelos. En el de efectos fijos los coeficientes van respecto a Andalucía, así que los centro restando su mediana para que sean comparables con los efectos aleatorios del multinivel, que ya están centrados en 0.

Mostrar código

coef_fe <- coef(mod_fe)[, grep("^ccaa", colnames(coef(mod_fe)))]
coef_fe <- cbind(ccaaAndalucía = 0, coef_fe) # Andalucía es la referencia
colnames(coef_fe) <- sub("^ccaa", "", colnames(coef_fe))

efectos_fe <- as_tibble(coef_fe, rownames = "partido") |>
  pivot_longer(-partido, names_to = "ccaa", values_to = "efecto") |>
  group_by(partido) |>
  mutate(efecto = efecto - median(efecto), modelo = "efectos fijos") |>
  ungroup()

re_ccaa <- ranef(mod_dem)$ccaa
efectos_ml <- map_dfr(voto5_lv[-1], function(p) {
  tibble(
    partido = p,
    ccaa = rownames(re_ccaa),
    efecto = re_ccaa[, "Estimate", paste0("mu", p, "_Intercept")],
    modelo = "multinivel"
  )
})

efectos <- bind_rows(efectos_fe, efectos_ml) |>
  left_join(n_ccaa |> mutate(ccaa = as.character(ccaa)), by = "ccaa") |>
  mutate(partido = factor(partido, levels = voto5_lv[-1]))

efectos |>
  filter(ccaa %in% c("Ceuta", "Melilla")) |>
  select(ccaa, partido, modelo, efecto) |>
  pivot_wider(names_from = modelo, values_from = efecto)
#> # A tibble: 8 × 4
#>   ccaa    partido `efectos fijos` multinivel
#>   <chr>   <fct>             <dbl>      <dbl>
#> 1 Ceuta   PSOE            -1.06      -0.189 
#> 2 Melilla PSOE             0          0.0462
#> 3 Ceuta   VOX              1.33       0.0943
#> 4 Melilla VOX             -0.770     -0.0175
#> 5 Ceuta   SUMAR           -0.0236    -0.0549
#> 6 Melilla SUMAR           -9.81      -0.150 
#> 7 Ceuta   OTROS            0.841      0.175 
#> 8 Melilla OTROS           -0.726     -0.316

En Melilla no hay ningún votante de SUMAR entre sus 17 entrevistas. El modelo de efectos fijos hace lo que le pedimos: busca el coeficiente que mejor ajusta esos datos. Y el que mejor ajusta un 0 es “menos infinito”. Se queda en −9,8 porque el optimizador se para. El multinivel no se cree que un 0 de 17 signifique que en Melilla nadie vota a SUMAR, y deja el efecto cerca de 0.

Ceuta tiene el problema contrario: 7 de sus 16 entrevistados votarían a VOX. Con efectos fijos, su efecto en VOX es 1,33. El multinivel lo deja en 0,09.

Lo vemos para todas las comunidades, ordenadas por tamaño muestral. Recorto el eje horizontal para que no se coma el gráfico el valor de Melilla.

Mostrar código

efectos |>
  mutate(ccaa = fct_reorder(ccaa, n_ccaa)) |>
  ggplot(aes(x = efecto, y = ccaa, color = modelo)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
  geom_line(aes(group = ccaa), color = "grey70") +
  geom_point(size = 2) +
  facet_wrap(~partido, nrow = 1) +
  coord_cartesian(xlim = c(-2.2, 2.2)) +
  scale_color_manual(
    values = c("efectos fijos" = "firebrick", "multinivel" = "steelblue")
  ) +
  labs(
    x = "Efecto de la comunidad (logit respecto al PP, centrado)",
    y = NULL,
    color = NULL,
    title = "Efecto de cada comunidad en los dos modelos",
    subtitle = "Ordenadas de menos (abajo) a más entrevistas (arriba)"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

Los puntos azules están más cerca del 0 que los rojos. Eso es el partial pooling, también llamado shrinkage o contracción: cada efecto de comunidad se acerca a la media de todas. Y no todas se acercan lo mismo. En el gráfico las comunidades con pocas entrevistas se quedan casi en el 0. Arriba, Cataluña y País Vasco mantienen buena parte de su efecto en OTROS, donde están ERC, Junts, PNV y Bildu.

Llamamos contracción a 1 − efecto multinivel / efecto fijos. Vale 0 si el multinivel deja el efecto como está y 1 si lo lleva a 0. Para los efectos fijos mayores de 0,2 en valor absoluto; con efectos casi nulos el cociente no significa nada.

Mostrar código

efectos |>
  select(ccaa, partido, modelo, efecto, n_ccaa) |>
  pivot_wider(names_from = modelo, values_from = efecto) |>
  filter(abs(`efectos fijos`) > 0.2) |>
  mutate(
    contraccion = 1 - multinivel / `efectos fijos`,
    entrevistas = cut(
      n_ccaa,
      c(0, 20, 100, 200, Inf),
      labels = c("≤ 20", "21-100", "101-200", "> 200")
    )
  ) |>
  group_by(entrevistas) |>
  summarise(
    contraccion_mediana = median(contraccion),
    n_efectos = n()
  ) |>
  gt() |>
  cols_label(
    entrevistas = "Entrevistas en la comunidad",
    contraccion_mediana = "Contracción mediana",
    n_efectos = "Nº de efectos"
  ) |>
  fmt_percent(contraccion_mediana, decimals = 0, locale = "es") |>
  cols_align("center", -entrevistas)

Entrevistas en la comunidad	Contracción mediana	Nº de efectos
≤ 20	88%	6
21-100	77%	12
101-200	74%	13
> 200	51%	8

Cuantas menos entrevistas, más contracción: casi un 90 % en Ceuta y Melilla, y la mitad en las comunidades con más de 200 entrevistas. El modelo decide cuánto contraer cada efecto con dos cosas:

Cuántos datos tiene la comunidad. Con pocas entrevistas, su estimación propia es poco fiable y el modelo se fía más de la media.
Cuánto varían las comunidades entre sí. Es la desviación típica de los efectos de comunidad, y también la estima el modelo.

Mostrar código

VarCorr(mod_dem)$ccaa$sd |> round(2)
#>                   Estimate Est.Error Q2.5 Q97.5
#> muPSOE_Intercept      0.26      0.07 0.14  0.43
#> muVOX_Intercept       0.13      0.09 0.01  0.33
#> muSUMAR_Intercept     0.29      0.11 0.10  0.53
#> muOTROS_Intercept     0.71      0.14 0.49  1.04

En OTROS la desviación típica es la mayor, 0,71. Las comunidades difieren de verdad en el voto a partidos regionales, y el modelo deja que se separen más.

Y en una celda concreta

Al final lo que vamos a usar son predicciones por celda. Tomo dos celdas de Ceuta en las que no hay ninguna entrevista:

Mostrar código

nd <- tibble(
  ccaa = c("Ceuta", "Ceuta"),
  sexo = c("H", "M"),
  edad = c("65+", "25-34"),
  edu = c("baja", "media"),
  n = 1
) |>
  mutate(
    ccaa = factor(ccaa, levels = levels(dat$ccaa)),
    sexo = factor(sexo, levels = levels(dat$sexo)),
    edad = factor(edad, levels = levels(dat$edad)),
    edu = factor(edu, levels = levels(dat$edu))
  )

nd
#> # A tibble: 2 × 5
#>   ccaa  sexo  edad  edu       n
#>   <fct> <fct> <fct> <fct> <dbl>
#> 1 Ceuta H     65+   baja      1
#> 2 Ceuta M     25-34 media     1

# ¿Hay alguien en estas celdas?
dat |> semi_join(nd, by = c("ccaa", "sexo", "edad", "edu")) |> nrow()
#> [1] 0

Ninguna de las dos tiene datos, pero los dos modelos pueden predecir. Combinan lo que saben de cada sexo, cada tramo de edad, cada nivel de estudios y de la gente de Ceuta.

Mostrar código

pred_fe <- predict(mod_fe, newdata = nd, type = "probs")
pred_ml <- apply(posterior_epred(mod_dem, newdata = nd), c(2, 3), mean)

nombre_celda <- c(
  "Ceuta · hombre · 65+ · estudios bajos",
  "Ceuta · mujer · 25-34 · estudios medios"
)

bind_rows(
  as_tibble(pred_fe) |> mutate(celda = nombre_celda, modelo = "efectos fijos"),
  as_tibble(pred_ml) |> mutate(celda = nombre_celda, modelo = "multinivel")
) |>
  arrange(celda) |>
  mutate(across(all_of(voto5_lv), ~ 100 * .x)) |>
  gt(groupname_col = "celda", rowname_col = "modelo") |>
  tab_spanner("% de voto predicho", columns = all_of(voto5_lv)) |>
  fmt_number(all_of(voto5_lv), decimals = 1, locale = "es") |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = c(VOX, SUMAR))
  )

	% de voto predicho
	PP	PSOE	VOX	SUMAR	OTROS
Ceuta · hombre · 65+ · estudios bajos
efectos fijos	28,0	10,9	34,8	5,8	20,4
multinivel	30,6	33,8	11,5	8,0	16,1
Ceuta · mujer · 25-34 · estudios medios
efectos fijos	7,9	5,3	58,7	5,2	22,9
multinivel	15,8	24,0	27,0	10,0	23,2

Las diferencias son grandes. El modelo de efectos fijos se cree del todo los 7 votantes de VOX de 16 entrevistas, y le da a VOX un 35 % en la primera celda y un 59 % en la segunda. El multinivel lo modera: VOX baja al 12 % y al 27 %.

¿Cuál está bien? No lo sabemos, porque en esas celdas no hay nadie. Pero un 59 % de VOX en una celda sin datos, sacado de 16 entrevistas de toda Ceuta, no es creíble.

Cuando sumemos cientos de celdas en la postestratificación, esos ceros y esos extremos se acumulan. Por eso la M de MRP es multinivel.

¿Y el recuerdo de voto?

El modelo que voy a usar para la estimación final sí lleva el recuerdo de voto. Hay dos razones:

Es, con diferencia, la variable que más explica la intención de voto actual.
Solo se puede postestratificar por las variables que están en el modelo. Si el modelo no usa el recuerdo, corregir en la población cuánta gente votó a cada partido en 2023 no cambia nada en la estimación.

La segunda razón es la importante para lo que viene: la encuesta del CIS tiene el recuerdo de voto bastante desviado respecto a lo que salió en las urnas.

La encuesta no se parece a la población

Una encuesta es una muestra, y nunca sale exactamente como la población. Unos grupos contestan más que otros, y el diseño muestral también reparte las entrevistas a su manera. Para ver cuánto se desvía, comparo la distribución de cada variable en la encuesta con la de la población.

Para la población uso tres fuentes:

Variable	Fuente
Sexo, edad y comunidad	Estadística Continua de Población (INE)
Nivel de estudios por sexo y comunidad	Encuesta de Población Activa (INE)
Recuerdo de voto por comunidad	Resultados del 23-J por municipio (Ministerio del Interior), más los que han cumplido 18 años desde entonces

Para el recuerdo uso los resultados por municipio y no por provincia. Los provinciales incluyen el voto de los residentes en el extranjero (CERA): 2,3 millones de personas en el censo, con una participación del 8,9 %. El CIS no los entrevista, y si los meto la abstención del objetivo sube casi cuatro puntos, del 28,2 % al 32,0 %.

Cada fuente tiene su universo: residentes, mayores de 16, censo electoral. En 03_targets.R las reescalo todas al mismo total dentro de cada comunidad: el electorado de hoy, 36,83 millones de personas.

Mostrar código

tg <- readRDS(file.path(dir_datos, "targets.rds"))
names(tg)
#> [1] "sexo_edad_ccaa" "edu_sexo_ccaa"  "rec_ccaa"       "total_ccaa"

head(tg$rec_ccaa)
#> # A tibble: 6 × 3
#>   ccaa      rec         N
#>   <fct>     <fct>   <dbl>
#> 1 Andalucía OTROS  179883
#> 2 Andalucía PP    1588654
#> 3 Andalucía PSOE  1459537
#> 4 Andalucía SUMAR  520834
#> 5 Andalucía VOX    668358
#> 6 Aragón    OTROS   46735

Cada elemento de tg es una tabla con cuánta gente hay en cada combinación de comunidad y otra variable. Por ejemplo, en Andalucía hay 1,6 millones de personas que votaron al PP en el 23-J.

Ahora comparo. Para la encuesta uso las 4.018 entrevistas, no solo las de voto válido, porque lo que me interesa aquí es a quién ha entrevistado el CIS. Añado también la distribución con los pesos que publica el CIS en la variable PESO.

Mostrar código

marginal_pob <- function(d, v) {
  d |>
    group_by(categoria = as.character(.data[[v]])) |>
    summarise(N = sum(N), .groups = "drop") |>
    mutate(poblacion = 100 * N / sum(N)) |>
    select(-N)
}

marginal_enc <- function(v) {
  cis |>
    group_by(categoria = as.character(.data[[v]])) |>
    summarise(
      encuesta = n(),
      ponderada_cis = sum(peso),
      .groups = "drop"
    ) |>
    mutate(
      encuesta = 100 * encuesta / sum(encuesta),
      ponderada_cis = 100 * ponderada_cis / sum(ponderada_cis)
    )
}

fuentes <- list(
  ccaa = tg$sexo_edad_ccaa,
  sexo = tg$sexo_edad_ccaa,
  edad = tg$sexo_edad_ccaa,
  edu = tg$edu_sexo_ccaa,
  rec = tg$rec_ccaa
)

marginales <- imap_dfr(fuentes, function(d, v) {
  marginal_enc(v) |>
    left_join(marginal_pob(d, v), by = "categoria") |>
    mutate(variable = v, .before = 1)
})

marginales |>
  filter(variable %in% c("edu", "rec")) |>
  mutate(
    variable = if_else(
      variable == "edu",
      "Nivel de estudios",
      "Recuerdo de voto"
    )
  ) |>
  gt(groupname_col = "variable", rowname_col = "categoria") |>
  cols_label(
    encuesta = "Encuesta (bruta)",
    ponderada_cis = "Encuesta (pesos CIS)",
    poblacion = "Población"
  ) |>
  tab_spanner("%", columns = c(encuesta, ponderada_cis, poblacion)) |>
  fmt_number(decimals = 1, locale = "es") |>
  tab_style(
    style = cell_fill(color = "#fde0dd"),
    locations = cells_body(
      columns = c(encuesta, ponderada_cis, poblacion),
      rows = abs(encuesta - poblacion) > 5
    )
  ) |>
  tab_footnote(
    "Resaltadas, las categorías en las que la encuesta bruta se aleja más de 5 puntos de la población."
  )

	%
	Encuesta (bruta)	Encuesta (pesos CIS)	Población
Nivel de estudios
alta	64,5	37,2	34,2
baja	14,0	39,7	42,5
media	21,5	23,1	23,3
Recuerdo de voto
ABST	12,4	15,6	28,2
BILDU	1,2	0,9	0,9
BNG	0,8	0,7	0,4
CCA	0,3	0,3	0,3
ERC	1,9	1,7	1,3
JUNTS	1,2	0,9	1,1
NO_PODIA	2,7	4,2	4,6
OTROS	6,8	5,7	3,0
PNV	0,6	0,4	0,7
PP	22,0	19,5	22,0
PSOE	29,5	31,0	21,1
SUMAR	12,1	10,5	8,2
VOX	8,7	8,6	8,2
Resaltadas, las categorías en las que la encuesta bruta se aleja más de 5 puntos de la población.

Mostrar código

orden <- c(
  "H",
  "M",
  "18-20",
  "21-24",
  "25-34",
  "35-44",
  "45-54",
  "55-64",
  "65+",
  "baja",
  "media",
  "alta",
  "PP",
  "PSOE",
  "VOX",
  "SUMAR",
  "ERC",
  "JUNTS",
  "BILDU",
  "PNV",
  "BNG",
  "CCA",
  "OTROS",
  "ABST",
  "NO_PODIA"
)

marginales |>
  filter(variable != "ccaa") |>
  pivot_longer(
    c(encuesta, ponderada_cis, poblacion),
    names_to = "fuente",
    values_to = "pct"
  ) |>
  mutate(
    categoria = factor(categoria, levels = orden),
    variable = factor(variable, levels = c("sexo", "edad", "edu", "rec")),
    fuente = factor(
      fuente,
      levels = c("encuesta", "ponderada_cis", "poblacion"),
      labels = c("encuesta (bruta)", "encuesta (pesos CIS)", "población")
    )
  ) |>
  ggplot(aes(x = categoria, y = pct, fill = fuente)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  facet_wrap(~variable, scales = "free", ncol = 2) +
  scale_fill_manual(values = c("grey70", "steelblue", "firebrick")) +
  labs(
    x = NULL,
    y = "%",
    fill = NULL,
    title = "Distribución de cada variable: encuesta frente a población",
    subtitle = "Población: electorado actual (ECP, EPA y resultados del 23-J)"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

Lo que salta a la vista:

Estudios. Es la desviación más grande. En la encuesta, el 64,5 % tiene estudios altos. En la población, el 34,2 %. Las personas con estudios bajos son el 42,5 % de la población y solo el 14,0 % de la encuesta. Es un patrón conocido: quien tiene más estudios contesta más a las encuestas.
Recuerdo de voto. Ojo, aquí hablo de lo que la gente dice que votó en 2023, no de lo que votaría ahora. En la encuesta, el 29,5 % dice que votó al PSOE en 2023. En el censo, el PSOE sacó votos equivalentes al 21,1 %. Y solo el 12,4 % dice que se abstuvo, frente al 28,2 % real. SUMAR y OTROS también salen por encima.
Edad y sexo. Hay desviaciones, pero pequeñas: faltan jóvenes de 18 a 24 años y mayores de 65.

Los pesos del CIS arreglan bien los estudios, la edad y el sexo. Pero no arreglan el recuerdo de voto: con los pesos, el PSOE sigue en un 31,0 % y la abstención en un 15,6 %, muy lejos del 21,1 % y el 28,2 % de las urnas.

Por comunidad, las desviaciones vienen del diseño:

Mostrar código

marginales |>
  filter(variable == "ccaa") |>
  select(-variable) |>
  mutate(razon = encuesta / poblacion) |>
  arrange(desc(razon)) |>
  gt(rowname_col = "categoria") |>
  cols_label(
    encuesta = "Encuesta (bruta)",
    ponderada_cis = "Encuesta (pesos CIS)",
    poblacion = "Población",
    razon = "Encuesta / población"
  ) |>
  tab_spanner("%", columns = c(encuesta, ponderada_cis, poblacion)) |>
  fmt_number(
    c(encuesta, ponderada_cis, poblacion),
    decimals = 1,
    locale = "es"
  ) |>
  fmt_number(razon, decimals = 1, locale = "es", pattern = "{x}×") |>
  data_color(razon, palette = c("white", "firebrick"), domain = c(0.5, 4))

	%			Encuesta / población
	Encuesta (bruta)	Encuesta (pesos CIS)	Población	Encuesta / población
La Rioja	2,5	0,7	0,7	3,7×
Melilla	0,5	0,1	0,2	3,1×
Ceuta	0,5	0,2	0,2	3,0×
Cantabria	2,4	1,4	1,3	1,8×
Navarra	2,4	1,5	1,4	1,7×
Baleares	2,5	2,5	2,3	1,1×
Extremadura	2,5	2,6	2,4	1,0×
Asturias	2,4	2,3	2,3	1,0×
Canarias	4,6	5,4	4,6	1,0×
Madrid	13,7	12,9	13,8	1,0×
Cataluña	15,2	15,6	15,5	1,0×
Castilla-La Mancha	4,3	4,7	4,4	1,0×
Castilla y León	5,1	5,5	5,4	0,9×
Comunitat Valenciana	9,7	11,2	10,3	0,9×
Aragón	2,6	2,8	2,8	0,9×
Murcia	2,8	3,2	3,1	0,9×
Andalucía	16,4	17,4	18,2	0,9×
Galicia	5,6	5,7	6,2	0,9×
País Vasco	4,2	4,3	4,8	0,9×

El CIS entrevista de más a las comunidades pequeñas (La Rioja, Ceuta, Melilla, Cantabria, Navarra). Lo hace a propósito, para tener un mínimo de entrevistas en cada una, y luego lo corrige con los pesos.

Entonces, ¿qué hago?

El modelo multinivel ya usa sexo, edad, estudios, comunidad y recuerdo. Así que, si la encuesta tiene demasiada gente con estudios altos o demasiados votantes del PSOE, el modelo no se equivoca por eso al predecir en cada celda. El problema está en cuánto pesa cada celda al sumar.

Ese peso tiene que salir de la población, no de la encuesta. Y es lo que hace la P de MRP: multiplicar la predicción de cada celda por cuánta gente hay en esa celda en la población real. Para eso hace falta saber cuánta gente hay en cada combinación de comunidad, sexo, edad, estudios y recuerdo. Esa tabla no la publica nadie. Solo tenemos las tablas por separado que acabamos de ver, y con ellas la vamos a construir mediante raking.

La P: postestratificar

La postestratificación es una media ponderada. Para cada celda $j$ (una combinación de comunidad, sexo, edad, estudios y recuerdo) tengo dos cosas:

$\theta_j$: la probabilidad de votar a cada partido que predice el modelo para esa celda.
$N_j$: cuánta gente hay en esa celda en la población.

La estimación para toda la población es:

\[ \hat\theta = \frac{\sum_j N_j \, \theta_j}{\sum_j N_j} \]

Es decir: cuántos votantes de cada partido hay en cada celda ($N_j \theta_j$), sumados para todas las celdas y divididos entre la población total. Si quiero la estimación de una comunidad, hago la misma cuenta solo con las celdas de esa comunidad.

Los $\theta_j$ salen del modelo. Los $N_j$ deberían salir de la población. Y aquí está el problema: necesito saber cuánta gente hay en cada combinación de comunidad × sexo × edad × estudios × recuerdo, y esa tabla no existe. El INE no la publica, y el recuerdo de voto no sale en ningún censo. Lo que sí tengo son las tablas de la sección anterior, por separado:

sexo × edad × comunidad (ECP),
estudios × sexo × comunidad (EPA),
recuerdo × comunidad (resultados del 23-J).

Necesito una tabla conjunta que cumpla esas tres a la vez. La construyo con raking, también llamado IPF (Iterative Proportional Fitting).

Raking con dos variables

Antes de ir a las cinco dimensiones, lo veo con dos: estudios y recuerdo. Parto de la tabla cruzada de la encuesta, que llamo semilla:

Mostrar código

library(mipfp)

semilla_2d <- xtabs(~ edu + rec, cis)
semilla_2d
#>        rec
#> edu      PP PSOE VOX SUMAR ERC JUNTS BILDU PNV BNG CCA OTROS ABST NO_PODIA
#>   baja  107  206  40    43   8     3     2   0   2   1    26  114       10
#>   media 168  233  81    97  14    11     9   7   5   3    41  107       89
#>   alta  608  745 228   346  54    33    36  17  24   8   205  279        8

Y tengo los dos marginales de la población, que es lo que la tabla tiene que cumplir:

Mostrar código

obj_edu <- xtabs(N ~ edu, tg$edu_sexo_ccaa)
obj_rec <- xtabs(N ~ rec, tg$rec_ccaa)

obj_edu
#> edu
#>     baja    media     alta 
#> 15638656  8598849 12592200
obj_rec
#> rec
#>       PP     PSOE      VOX    SUMAR      ERC    JUNTS    BILDU      PNV 
#>  8093469  7762316  3034117  3014755   462890   392609   333321   275859 
#>      BNG      CCA    OTROS     ABST NO_PODIA 
#>   152348   114707  1107831 10396976  1688507

El IPF hace algo muy sencillo:

Reescala las filas de la semilla para que sumen el marginal de estudios.
Reescala las columnas para que sumen el marginal de recuerdo. Al hacerlo, las filas dejan de cuadrar un poco.
Repite los pasos 1 y 2 hasta que filas y columnas cuadran a la vez.

Con mipfp::Ipfp le paso la semilla, qué dimensiones tiene cada objetivo y los objetivos. Trabajo en proporciones y luego multiplico por el total, porque Ipfp exige que todos los objetivos sumen exactamente lo mismo y los redondeos dan problemas.

Mostrar código

ipf_2d <- Ipfp(
  seed = semilla_2d / sum(semilla_2d),
  target.list = list(1, 2), # dimensión 1 = edu, dimensión 2 = rec
  target.data = list(obj_edu / sum(obj_edu), obj_rec / sum(obj_rec)),
  tol = 1e-10
)

tabla_2d <- ipf_2d$x.hat * sum(obj_rec)
round(tabla_2d / 1000) # en miles de personas
#>        rec
#> edu       PP PSOE  VOX SUMAR  ERC JUNTS BILDU  PNV  BNG  CCA OTROS ABST
#>   baja  3060 3707 1079   902  159    89    55    0   36   32   361 5769
#>   media 1718 1499  782   728   99   117    89  120   33   34   204 1936
#>   alta  3315 2556 1173  1384  205   187   189  156   83   49   543 2692
#>        rec
#> edu     NO_PODIA
#>   baja       389
#>   media     1240
#>   alta        59

La tabla ya suma por filas lo que dice la EPA y por columnas lo que dicen las urnas. ¿Qué se queda de la encuesta? La asociación entre las dos variables. Lo veo con la razón de odds (odds ratio) de votar al PP frente al PSOE entre estudios altos y bajos:

Mostrar código

odds_ratio <- function(t) {
  (t["alta", "PP"] * t["baja", "PSOE"]) / (t["alta", "PSOE"] * t["baja", "PP"])
}

c(semilla = odds_ratio(semilla_2d), raking = odds_ratio(tabla_2d))
#>  semilla   raking 
#> 1.571197 1.571197

Es exactamente la misma. El raking cambia los totales y respeta las asociaciones que había en la semilla. Esa es su hipótesis de fondo: la estructura de la encuesta es buena, lo que falla son las proporciones. Lo que no está en ningún marginal, como la relación entre estudios y recuerdo, lo pone la encuesta.

Raking con cinco variables

Con las cinco variables la idea es la misma, pero la semilla tiene 19 × 2 × 7 × 3 × 13 = 10.374 celdas, y la encuesta tiene 4.018 entrevistas:

Mostrar código

semilla_enc <- xtabs(~ ccaa + sexo + edad + edu + rec, cis)
dn <- dimnames(semilla_enc)

mean(semilla_enc == 0)
#> [1] 0.835936

El 84 % de las celdas de la semilla están vacías. Y el IPF solo multiplica: una celda que empieza a 0 se queda a 0. Si uso la encuesta tal cual, la tabla dirá que en la población no hay nadie en esas celdas, y casi todas tienen gente. Así que mezclo la encuesta con una tabla uniforme, a partes iguales:

Mostrar código

alpha <- 0.5
semilla <- alpha *
  array(1, dim(semilla_enc), dn) /
  length(semilla_enc) +
  (1 - alpha) * semilla_enc / sum(semilla_enc)

Hay celdas que sí tienen que valer 0: son los ceros estructurales. Hay dos tipos:

Por edad. Quien tiene hoy 18-20 años no pudo votar el 23-J, así que su recuerdo es siempre NO_PODIA. Y quien tiene más de 20 años no puede estar en NO_PODIA.
Por territorio. Nadie pudo votar a ERC fuera de Cataluña, ni al BNG fuera de Galicia.

Los pongo a mano y el IPF los respeta:

Mostrar código

i18 <- which(dn$edad == "18-20")
inp <- which(dn$rec == "NO_PODIA")

semilla[,, i18, , -inp] <- 0 # 18-20 con cualquier recuerdo que no sea NO_PODIA
semilla[,, -i18, , inp] <- 0 # NO_PODIA con más de 20 años

for (p in names(territorio)) {
  semilla[!dn$ccaa %in% territorio[[p]], , , , p] <- 0 # partido fuera de su comunidad
}

mean(semilla == 0)
#> [1] 0.5691151

Ahora los objetivos. En el ejemplo de dos variables eran marginales de una sola variable. Aquí no: las fuentes publican tablas cruzadas, y eso es oro para el IPF. El INE da la población por sexo y edad y comunidad, no cada variable por separado. Así que los objetivos son tablas conjuntas de tres y dos variables:

Mostrar código

t1 <- xtabs(N ~ ccaa + sexo + edad, tg$sexo_edad_ccaa) # ECP
t2 <- xtabs(N ~ ccaa + sexo + edu, tg$edu_sexo_ccaa) # EPA
t3 <- xtabs(N ~ ccaa + rec, tg$rec_ccaa) # 23-J

dim(t1)
#> [1] 19  2  7
dim(t2)
#> [1] 19  2  3
dim(t3)
#> [1] 19 13

t1 es un array de 19 × 2 × 7: para cada comunidad, una tabla de sexo por edad. Por ejemplo, Madrid, en miles de personas:

Mostrar código

round(t1["Madrid", , ] / 1000)
#>     edad
#> sexo 18-20 21-24 25-34 35-44 45-54 55-64 65+
#>    H   106   146   390   402   479   393 492
#>    M   102   143   396   419   507   437 683

t2 es lo mismo con estudios en lugar de edad:

Mostrar código

round(t2["Madrid", , ] / 1000)
#>     edu
#> sexo baja media alta
#>    H  805   574 1035
#>    M  922   671 1089

Y t3 es el recuerdo de voto de cada comunidad:

Mostrar código

round(t3["Madrid", ] / 1000)
#>       PP     PSOE      VOX    SUMAR      ERC    JUNTS    BILDU      PNV 
#>     1444      994      500      551        0        0        0        0 
#>      BNG      CCA    OTROS     ABST NO_PODIA 
#>        0        0      106     1253      247

Cuanta más información conjunta dan las fuentes, menos tiene que inventar el raking. Con cinco variables hay diez pares posibles. Para cada par, la relación entre las dos variables sale de una fuente real o, si ningún objetivo las cruza, de la semilla, es decir, de la encuesta:

Par de variables	De dónde sale la relación
ccaa × sexo	ECP y EPA
ccaa × edad	ECP
ccaa × edu	EPA
ccaa × rec	23-J (más los ceros de territorio)
sexo × edad	ECP
sexo × edu	EPA
sexo × rec	encuesta
edad × edu	encuesta
edad × rec	encuesta (más los ceros estructurales)
edu × rec	encuesta

Seis de los diez pares los fijan fuentes oficiales. Los otros cuatro los pone la encuesta. Si encontrara una fuente con, por ejemplo, estudios por edad, la añadiría como cuarto objetivo y la encuesta tendría que poner un par menos.

Al llamar a Ipfp le digo qué dimensiones de la semilla toca cada objetivo. El orden de las dimensiones de la semilla es ccaa (1), sexo (2), edad (3), edu (4) y rec (5):

Mostrar código

ipf <- Ipfp(
  seed = semilla / sum(semilla),
  target.list = list(
    c(1, 2, 3), # t1: ccaa, sexo, edad
    c(1, 2, 4), # t2: ccaa, sexo, edu
    c(1, 5) # t3: ccaa, rec
  ),
  target.data = list(t1 / sum(t1), t2 / sum(t2), t3 / sum(t3)),
  iter = 1000,
  tol = 1e-10
)

length(ipf$evol.stp.crit) # iteraciones
#> [1] 13

tabla <- ipf$x.hat * sum(tg$total_ccaa$total)
dim(tabla)
#> [1] 19  2  7  3 13

Warning

Ipfp empareja los objetivos con la semilla por posición, no por nombre. Si los niveles de un factor están en distinto orden en la semilla y en un objetivo, la tabla sale permutada y no da ningún aviso. Por eso en 03_targets.R guardo todas las variables como factores con el mismo orden de niveles.

Compruebo que la tabla cumple los marginales. Por ejemplo, el de recuerdo:

Mostrar código

tibble(
  rec = dn$rec,
  objetivo = as.numeric(apply(t3, 2, sum)),
  tabla = as.numeric(apply(tabla, 5, sum))
) |>
  gt(rowname_col = "rec") |>
  cols_label(objetivo = "Objetivo (23-J)", tabla = "Tabla sintética") |>
  tab_spanner("Personas", columns = c(objetivo, tabla)) |>
  fmt_number(decimals = 0, locale = "es")

	Personas
	Objetivo (23-J)	Tabla sintética
PP	8.093.469	8.093.469
PSOE	7.762.316	7.762.316
VOX	3.034.117	3.034.117
SUMAR	3.014.755	3.014.755
ERC	462.890	462.890
JUNTS	392.609	392.609
BILDU	333.321	333.321
PNV	275.859	275.859
BNG	152.348	152.348
CCA	114.707	114.707
OTROS	1.107.831	1.107.831
ABST	10.396.976	10.396.976
NO_PODIA	1.688.507	1.688.507

Cuadra. Esta es la tabla de postestratificación sintética: la población estimada de cada una de las 10.374 celdas. “Sintética” porque no la ha medido nadie, la he fabricado con tres tablas reales y la estructura de la encuesta.

Para usarla con el modelo la paso a data frame, con una fila por celda, y quito las celdas a 0, que son los ceros estructurales:

Mostrar código

post <- as.data.frame.table(tabla, responseName = "N") |>
  filter(N > 1e-6) |>
  mutate(n = 1)

nrow(post)
#> [1] 4470
head(post)
#>        ccaa sexo  edad  edu rec        N n
#> 1 Andalucía    H 21-24 baja  PP 7688.767 1
#> 2    Aragón    H 21-24 baja  PP 2220.426 1
#> 3  Asturias    H 21-24 baja  PP 2344.072 1
#> 4  Baleares    H 21-24 baja  PP 1334.045 1
#> 5  Canarias    H 21-24 baja  PP 1895.940 1
#> 6 Cantabria    H 21-24 baja  PP 1449.799 1

Quedan 4.470 celdas con población. Las 5.904 que faltan son los ceros estructurales. Cada fila es una celda y N es cuánta gente hay en ella. La columna n = 1 la necesita brms para predecir, porque el modelo tiene trials(n): con n = 1 la predicción es una probabilidad por partido.

Juntando la M y la P

El modelo completo

Ahora sí, el modelo lleva el recuerdo de voto y los once partidos. Es el mismo modelo multinivel de la primera parte con un efecto aleatorio más, (1 | rec), y ajustado sobre las celdas de las cinco variables.

Primero agrego la encuesta por celda:

Mostrar código

celdas_enc <- dat |>
  count(ccaa, sexo, edad, edu, rec, voto) |>
  pivot_wider(names_from = voto, values_from = n, values_fill = 0)

celdas_enc$y <- as.matrix(celdas_enc[, voto_lv])
celdas_enc$n <- rowSums(celdas_enc$y)

Aquí hay un detalle que conviene ver despacio. y no es una columna normal: es una matriz metida dentro de una columna del data frame. Tiene una fila por celda y una columna por partido:

Mostrar código

class(celdas_enc$y)
#> [1] "matrix" "array"
dim(celdas_enc$y)
#> [1] 1442   11
head(celdas_enc$y)
#>      PP PSOE VOX SUMAR ERC JUNTS BILDU PNV BNG CCA OTROS
#> [1,]  0    0   1     0   0     0     0   0   0   0     0
#> [2,]  1    0   1     0   0     0     0   0   0   0     0
#> [3,]  0    0   0     0   0     0     0   0   0   0     1
#> [4,]  0    0   1     0   0     0     0   0   0   0     0
#> [5,]  2    0   0     0   0     0     0   0   0   0     1
#> [6,]  0    0   0     0   0     0     0   0   0   0     1

Si imprimo el data frame, y aparece como un bloque de columnas y[,"PP"], [,"PSOE"], etc.:

Mostrar código

celdas_enc |>
  select(ccaa, sexo, edad, edu, rec, y, n) |>
  arrange(desc(n))
#> # A tibble: 1,442 × 7
#>    ccaa      sexo  edad  edu   rec   y[,"PP"] [,"PSOE"] [,"VOX"]     n
#>    <fct>     <fct> <fct> <fct> <fct>    <int>     <int>    <int> <dbl>
#>  1 Andalucía M     45-54 alta  PP          16         0        4    20
#>  2 Andalucía M     55-64 alta  PP          17         0        0    17
#>  3 Madrid    M     45-54 alta  PP          15         0        2    17
#>  4 Madrid    M     65+   baja  PSOE         2        15        0    17
#>  5 Andalucía M     55-64 alta  PSOE         1        13        0    15
#>  6 Andalucía M     35-44 alta  PP           9         0        3    14
#>  7 Cataluña  M     65+   alta  PSOE         0        12        0    14
#>  8 Madrid    H     35-44 alta  PP          11         0        2    14
#>  9 Madrid    H     55-64 alta  PP          12         1        1    14
#> 10 Andalucía M     65+   alta  PSOE         0        11        1    13
#> # ℹ 1,432 more rows
#> # ℹ 1 more variable: y[4:11] <int>

Lo veo mejor en una tabla interactiva. Las columnas en azul son la matriz y: cuántas personas de la celda votarían a cada partido. n es su suma, el total de entrevistas de la celda:

Mostrar código

celdas_enc_dt <- celdas_enc |>
  select(ccaa, sexo, edad, edu, rec) |>
  bind_cols(
    as_tibble(celdas_enc$y) |> rename_with(~ paste0('y[, "', .x, '"]')),
    n = celdas_enc$n
  ) |>
  arrange(desc(n))

DT::datatable(
  celdas_enc_dt,
  rownames = FALSE,
  filter = "top",
  options = list(pageLength = 10, scrollX = TRUE)
) |>
  DT::formatStyle(
    columns = paste0('y[, "', voto_lv, '"]'),
    backgroundColor = "#dbe9f6"
  )

Esa es la forma que pide brms para un modelo multinomial con datos agregados. En la fórmula, y | trials(n) significa: en cada celda hay n personas, y y dice cuántas eligen cada categoría. Es la versión multinomial de la binomial con éxitos | trials(n).

Queda un detalle: los partidos regionales. Si ajusto el modelo tal cual, el modelo ve que en 18 comunidades nadie vota a ERC y lo trata como información: el efecto de comunidad de ERC tiene una media muy negativa, y el partial pooling acerca Cataluña a esa media. Es decir, subestimaría a ERC justo donde se presenta.

Pero esos ceros no son información sobre ERC, son ceros estructurales. Los meto en el modelo con un offset: para cada partido regional, una columna que vale 0 en las comunidades donde se presenta y −20 en las demás. Se suma al predictor lineal de ese partido, así que fuera de su territorio su probabilidad queda multiplicada por $e^{-20} \approx 2 \cdot 10^{-9}$. Es decir, prácticamente cero. Luego veremos que no siempre, y cómo se arregla.

Mostrar código

con_offsets <- function(d) {
  for (p in names(territorio)) {
    d[[paste0("off_", p)]] <- ifelse(d$ccaa %in% territorio[[p]], 0, -20)
  }
  d
}

celdas_enc <- con_offsets(celdas_enc)

celdas_enc |>
  filter(ccaa %in% c("Cataluña", "Andalucía")) |>
  distinct(ccaa, across(starts_with("off_")))
#> # A tibble: 2 × 7
#>   ccaa      off_ERC off_JUNTS off_BILDU off_PNV off_BNG off_CCA
#>   <fct>       <dbl>     <dbl>     <dbl>   <dbl>   <dbl>   <dbl>
#> 1 Andalucía     -20       -20       -20     -20     -20     -20
#> 2 Cataluña        0         0       -20     -20     -20     -20

brms permite una fórmula distinta para cada partido (cada mu). Todos llevan los mismos efectos, y los regionales además su offset:

Mostrar código

f_comun <- "sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa)"

f_partidos <- lapply(voto_lv[-1], function(p) {
  offset <- if (p %in% names(territorio)) paste0(" + offset(off_", p, ")") else ""
  as.formula(paste0("mu", p, " ~ ", f_comun, offset))
})

formula_mrp <- do.call(
  bf,
  c(list(as.formula(paste("y | trials(n) ~", f_comun))), f_partidos)
)
formula_mrp
#> y | trials(n) ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) 
#> muPSOE ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa)
#> muVOX ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa)
#> muSUMAR ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa)
#> muERC ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_ERC)
#> muJUNTS ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_JUNTS)
#> muBILDU ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_BILDU)
#> muPNV ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_PNV)
#> muBNG ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_BNG)
#> muCCA ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa) + offset(off_CCA)
#> muOTROS ~ sexo + (1 | edad) + (1 | edu) + (1 | rec) + (1 | ccaa)

Mostrar código

mod <- brm(
  formula_mrp,
  data = celdas_enc,
  family = multinomial(),
  prior = priors_multinom(voto_lv),
  chains = 4,
  cores = 4,
  iter = 2000,
  refresh = 0,
  silent = 2,
  control = list(adapt_delta = 0.95),
  backend = "cmdstanr",
  seed = 2026,
  file = file.path(dir_datos, "mod_mrp")
)

Una celda de la tabla

Tomo la primera fila de la tabla de postestratificación:

Mostrar código

post[1, ]
#>        ccaa sexo  edad  edu rec        N n
#> 1 Andalucía    H 21-24 baja  PP 7688.767 1

Es una celda: hombres de Andalucía, de 21 a 24 años, con estudios bajos, que votaron al PP en el 23-J. Según la tabla sintética hay unas 7.689 personas así. Eso es su $N_j$.

El modelo toma de esa fila las variables ccaa, sexo, edad, edu y rec, y predice la probabilidad de votar a cada partido. Pero no da una sola predicción. Es un modelo bayesiano, así que tengo 4.000 muestras de la posterior (4 cadenas × 1.000 iteraciones). Me quedo con 500, que bastan:

Mostrar código

post <- con_offsets(post) # la tabla también necesita los offsets

set.seed(2026)
ep <- posterior_epred(mod, newdata = post, ndraws = 500)

dim(ep)
#> [1]  500 4470   11

ep es un array de tres dimensiones:

500 draws de la posterior.
4.470 celdas, una por fila de post.
11 partidos.

Antes de seguir, los ceros estructurales. El offset deja la probabilidad de cada partido regional fuera de su territorio en casi 0, pero no siempre en 0 exacto. Fuera de Cataluña ningún dato informa del efecto de comunidad de ERC, así que el modelo lo saca de la distribución común. En algún draw ese efecto sale tan grande que compensa el −20 del offset:

Mostrar código

mascara <- sapply(voto_lv, function(p) {
  if (p %in% names(territorio)) {
    as.numeric(post$ccaa %in% territorio[[p]])
  } else {
    rep(1, nrow(post))
  }
})

# Probabilidad máxima de cada partido regional fuera de su territorio
fuga <- apply(ep, c(2, 3), max) * (1 - mascara)
round(apply(fuga[, names(territorio)], 2, max), 4)
#>    ERC  JUNTS  BILDU    PNV    BNG    CCA 
#> 0.0069 0.0053 0.0355 0.0010 0.0062 0.0006

Así que en cada draw pongo esas probabilidades a 0 exacto y renormalizo cada celda para que siga sumando 1:

Mostrar código

for (d in seq_len(dim(ep)[1])) {
  ep[d, , ] <- ep[d, , ] * mascara
  ep[d, , ] <- ep[d, , ] / rowSums(ep[d, , ])
}

max(apply(ep, c(2, 3), max) * (1 - mascara))
#> [1] 0

Las cinco primeras muestras para la primera celda:

Mostrar código

round(ep[1:5, 1, ], 3)
#>         PP  PSOE   VOX SUMAR ERC JUNTS BILDU PNV BNG CCA OTROS
#> [1,] 0.769 0.001 0.116 0.001   0     0     0   0   0   0 0.113
#> [2,] 0.802 0.001 0.145 0.000   0     0     0   0   0   0 0.051
#> [3,] 0.665 0.001 0.249 0.001   0     0     0   0   0   0 0.084
#> [4,] 0.780 0.003 0.147 0.002   0     0     0   0   0   0 0.068
#> [5,] 0.765 0.002 0.162 0.000   0     0     0   0   0   0 0.070

Cada fila es una muestra de la posterior, y suma 1. Todas dicen lo mismo con pequeñas diferencias: esta gente vota mayoritariamente al PP, con una parte que se va a VOX. Esa variación entre filas es la incertidumbre del modelo sobre esta celda. La media de las 500 es la $\theta_j$ de esta celda:

Mostrar código

theta_1 <- colMeans(ep[, 1, ])
round(theta_1, 3)
#>    PP  PSOE   VOX SUMAR   ERC JUNTS BILDU   PNV   BNG   CCA OTROS 
#> 0.704 0.002 0.219 0.001 0.000 0.000 0.000 0.000 0.000 0.000 0.075

La P a mano

Ahora multiplico por la $N_j$ de la celda. Así obtengo cuántas personas de esa celda votarían a cada partido:

Mostrar código

round(theta_1 * post$N[1])
#>    PP  PSOE   VOX SUMAR   ERC JUNTS BILDU   PNV   BNG   CCA OTROS 
#>  5411    15  1682     6     0     0     0     0     0     0   575

De los 7.689 hombres jóvenes andaluces con estudios bajos que votaron al PP, unos 5.400 lo seguirían votando y unos 1.700 se irían a VOX. Los partidos regionales salen con 0 exacto: en Andalucía no se presentan.

La postestratificación es hacer esto en las 4.218 celdas, sumar y dividir entre la población total. Para un solo draw, el primero:

Mostrar código

votantes_draw1 <- colSums(ep[1, , ] * post$N) # suma de N_j * theta_j
round(votantes_draw1)
#>       PP     PSOE      VOX    SUMAR      ERC    JUNTS    BILDU      PNV 
#> 10092595  8677805  5946717  3055844   753923   277465   608407   363808 
#>      BNG      CCA    OTROS 
#>   215899    76209  6761034

round(100 * votantes_draw1 / sum(post$N), 1)
#>    PP  PSOE   VOX SUMAR   ERC JUNTS BILDU   PNV   BNG   CCA OTROS 
#>  27.4  23.6  16.1   8.3   2.0   0.8   1.7   1.0   0.6   0.2  18.4

ep[1, , ] es una matriz de 4.470 celdas × 11 partidos. Al multiplicarla por post$N, cada fila se multiplica por la población de su celda. colSums suma por partido. Esa es exactamente la fórmula $\sum_j N_j \theta_j / \sum_j N_j$.

Y para tener la incertidumbre, repito lo mismo en cada uno de los 500 draws:

Mostrar código

mrp_draws <- apply(ep, 1, function(p) colSums(p * post$N) / sum(post$N)) |>
  t()

dim(mrp_draws)
#> [1] 500  11

est_mrp <- tibble(
  partido = voto_lv,
  mrp = 100 * colMeans(mrp_draws),
  q05 = 100 * apply(mrp_draws, 2, quantile, 0.05),
  q95 = 100 * apply(mrp_draws, 2, quantile, 0.95)
)

est_mrp |> mutate(across(-partido, ~ round(.x, 1)))
#> # A tibble: 11 × 4
#>    partido   mrp   q05   q95
#>    <chr>   <dbl> <dbl> <dbl>
#>  1 PP       26.8  25.2  28.5
#>  2 PSOE     24.2  22.9  25.7
#>  3 VOX      16.9  15.3  18.5
#>  4 SUMAR     8.6   7.7   9.6
#>  5 ERC       1.8   1.4   2.3
#>  6 JUNTS     1     0.7   1.4
#>  7 BILDU     1.4   1.2   1.8
#>  8 PNV       0.9   0.7   1.3
#>  9 BNG       0.6   0.5   0.9
#> 10 CCA       0.2   0.1   0.4
#> 11 OTROS    17.5  15.7  19.3

Ya está. Eso es MRP: 500 estimaciones de la intención de voto en toda la población, de las que saco la media y un intervalo al 90 %.

¿Cuánto cambia?

Lo comparo con la intención de voto bruta de la encuesta, con la ponderada con los pesos del CIS y con la estimación que publicó el CIS para este barómetro:

Mostrar código

# Estimación del CIS (página 3 de 3577_Estimacion.pdf), agrupada como aquí:
# SUMAR = SUMAR + Podemos; OTROS = SALF + UPN + otros partidos + en blanco
cis_publicado <- c(
  PP = 25.5,
  PSOE = 31.0,
  VOX = 16.6,
  SUMAR = 5.7 + 3.7,
  ERC = 2.5,
  JUNTS = 0.7,
  BILDU = 1.3,
  PNV = 0.7,
  BNG = 0.7,
  CCA = 0.2,
  OTROS = 1.8 + 0.1 + 8.6 + 1.0
)

comparacion <- est_mrp |>
  mutate(
    bruto = as.numeric(100 * prop.table(table(dat$voto))[partido]),
    pesos_cis = as.numeric(
      100 * tapply(dat$peso, dat$voto, sum)[partido] / sum(dat$peso)
    ),
    cis_publicado = cis_publicado[partido]
  ) |>
  select(partido, bruto, pesos_cis, mrp, q05, q95, cis_publicado)

comparacion |>
  mutate(partido = if_else(partido == "SUMAR", "SUMAR y aliados", partido)) |>
  gt(rowname_col = "partido") |>
  cols_merge(c(q05, q95), pattern = "{1} – {2}") |>
  cols_label(
    bruto = "Bruto",
    pesos_cis = "Pesos CIS",
    mrp = "MRP",
    q05 = "Intervalo 90 %",
    cis_publicado = "CIS publicado"
  ) |>
  tab_spanner("% sobre voto válido", columns = everything()) |>
  fmt_number(decimals = 1, locale = "es") |>
  tab_style(
    style = list(cell_text(weight = "bold"), cell_fill(color = "#fde0dd")),
    locations = cells_body(columns = mrp)
  )

	% sobre voto válido
	Bruto	Pesos CIS	MRP	Intervalo 90 %	CIS publicado
PP	25,7	24,2	26,8	25,2 – 28,5	25,5
PSOE	28,4	30,6	24,2	22,9 – 25,7	31,0
VOX	14,5	15,7	16,9	15,3 – 18,5	16,6
SUMAR y aliados	10,6	10,0	8,6	7,7 – 9,6	9,4
ERC	2,0	2,0	1,8	1,4 – 2,3	2,5
JUNTS	0,9	0,7	1,0	0,7 – 1,4	0,7
BILDU	1,6	1,3	1,4	1,2 – 1,8	1,3
PNV	0,6	0,6	0,9	0,7 – 1,3	0,7
BNG	1,0	0,8	0,6	0,5 – 0,9	0,7
CCA	0,2	0,1	0,2	0,1 – 0,4	0,2
OTROS	14,5	13,9	17,5	15,7 – 19,3	11,5

Mostrar código

comparacion |>
  pivot_longer(
    c(bruto, pesos_cis, mrp, cis_publicado),
    names_to = "metodo",
    values_to = "pct"
  ) |>
  mutate(
    metodo = factor(
      metodo,
      levels = c("bruto", "pesos_cis", "mrp", "cis_publicado"),
      labels = c("bruto", "pesos CIS", "MRP", "CIS publicado")
    ),
    partido = factor(
      partido,
      levels = rev(voto_lv),
      labels = rev(c(
        "PP", "PSOE", "VOX", "SUMAR y aliados",
        "ERC", "Junts", "EH Bildu", "PNV", "BNG", "CC", "OTROS"
      ))
    ),
    q05 = if_else(metodo == "MRP", q05, NA_real_),
    q95 = if_else(metodo == "MRP", q95, NA_real_)
  ) |>
  ggplot(aes(x = pct, y = partido, color = metodo)) +
  geom_pointrange(
    aes(xmin = q05, xmax = q95),
    position = position_dodge(width = 0.6),
    na.rm = TRUE
  ) +
  scale_color_manual(values = c("grey60", "steelblue", "firebrick", "black")) +
  labs(
    x = "% sobre voto válido",
    y = NULL,
    color = NULL,
    title = "Intención de voto según el método",
    subtitle = "MRP con intervalo al 90 %"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

Donde más se nota es en el PSOE:

	PSOE
Bruto	28,4 %
Pesos del CIS	30,6 %
CIS publicado	31,0 %
MRP	24,2 % (22,9 a 25,7)

El MRP le quita unos 7 puntos respecto a lo que publicó el CIS. Mi interpretación es que viene del recuerdo de voto. En la encuesta, el 29,5 % dice que votó al PSOE en 2023, y en las urnas fue el 21,1 % del electorado. Los pesos del CIS no lo corrigen. Al postestratificar con la tabla sintética, cada celda de votantes del PSOE en 2023 pesa lo que pesa en la población, no lo que pesa en la encuesta. El PP y VOX suben un poco por el mismo motivo, y SUMAR baja.

La otra diferencia grande está en OTROS: 14,5 % en bruto, 11,5 % en la estimación del CIS y 17,5 % con MRP. Son dos efectos distintos:

Bruto frente al CIS. El 4,2 % del voto válido de la encuesta es voto en blanco, y el CIS lo deja en un 1,0 % en su estimación (página 3 del PDF de la estimación). Eso explica casi toda la diferencia. Aquí no hay “cocina”, así que el blanco se queda como lo dice la gente.
MRP frente al bruto. Quien se abstuvo en 2023 y ahora dice que votaría, se va mucho a OTROS: un 33 % según el modelo. En la encuesta esos abstencionistas son el 7,5 % de los que declaran voto válido. En la tabla de postestratificación son el 28,2 % del electorado. Al darles su peso real, OTROS sube.

El punto 2 tiene trampa, y la cuento en los avisos del final.

Por comunidad

El CIS no publica estimaciones por comunidad, y con razón: con 61 entrevistas en Navarra o 133 en el País Vasco, el voto bruto no dice gran cosa. El MRP sí puede: cada comunidad toma prestada información del resto para todo lo que no es propio de ella, y la tabla de postestratificación pone el peso de cada celda. Aquí es donde se nota tener a los partidos regionales separados.

Hago la misma cuenta $\sum_j N_j \theta_j / \sum_j N_j$, pero solo con las celdas de cada comunidad. Al lado pongo el resultado del 23-J como referencia, no como la verdad: es de hace tres años.

Mostrar código

ccaa_reg <- c("Cataluña", "País Vasco", "Navarra", "Galicia", "Canarias")

# La misma cuenta que para España, pero solo con las celdas de cada comunidad
mrp_ccaa_draws <- map(ccaa_reg, function(cc) {
  i <- post$ccaa == cc
  apply(ep[, i, ], 1, function(p) colSums(p * post$N[i]) / sum(post$N[i])) |>
    t()
}) |>
  set_names(ccaa_reg)

# Referencia: el 23-J en cada comunidad, sobre votos emitidos
res_23j <- tg$rec_ccaa |>
  filter(ccaa %in% ccaa_reg, !rec %in% c("ABST", "NO_PODIA")) |>
  group_by(ccaa) |>
  mutate(pct_23j = 100 * N / sum(N)) |>
  ungroup() |>
  transmute(ccaa = as.character(ccaa), partido = as.character(rec), pct_23j)

est_ccaa <- imap_dfr(mrp_ccaa_draws, function(d, cc) {
  tibble(
    ccaa = cc,
    partido = voto_lv,
    mrp = 100 * colMeans(d),
    q05 = 100 * apply(d, 2, quantile, 0.05),
    q95 = 100 * apply(d, 2, quantile, 0.95)
  )
}) |>
  left_join(res_23j, by = c("ccaa", "partido")) |>
  mutate(pct_23j = replace_na(pct_23j, 0)) |>
  filter(mrp >= 1 | pct_23j >= 1) |>
  left_join(
    dat |> count(ccaa, name = "entrevistas") |> mutate(ccaa = as.character(ccaa)),
    by = "ccaa"
  ) |>
  mutate(
    partido = if_else(partido == "SUMAR", "SUMAR y aliados", partido),
    grupo = paste0(ccaa, " (", entrevistas, " entrevistas)")
  )

est_ccaa |>
  select(grupo, partido, mrp, q05, q95, pct_23j) |>
  gt(groupname_col = "grupo", rowname_col = "partido") |>
  cols_merge(c(q05, q95), pattern = "{1} – {2}") |>
  cols_label(
    mrp = "MRP",
    q05 = "Intervalo 90 %",
    pct_23j = "23-J"
  ) |>
  tab_spanner("% de voto", columns = everything()) |>
  fmt_number(decimals = 1, locale = "es") |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = mrp)
  ) |>
  tab_footnote("23-J: porcentaje sobre votos emitidos (OTROS incluye blanco y nulo).")

	% de voto
	MRP	Intervalo 90 %	23-J
Cataluña (458 entrevistas)
PP	15,9	14,0 – 17,8	13,2
PSOE	25,3	23,0 – 27,9	34,2
VOX	12,4	10,4 – 14,7	7,7
SUMAR y aliados	9,3	7,7 – 11,4	13,9
ERC	11,4	8,8 – 14,5	13,0
JUNTS	6,4	4,7 – 8,7	11,1
OTROS	19,3	16,5 – 22,4	6,9
País Vasco (133 entrevistas)
PP	12,7	10,2 – 15,3	11,4
PSOE	18,3	15,7 – 20,8	25,1
VOX	7,3	5,2 – 9,8	2,6
SUMAR y aliados	7,2	5,5 – 9,4	11,0
BILDU	24,1	19,4 – 29,6	23,8
PNV	19,3	13,8 – 26,5	23,9
OTROS	11,2	8,1 – 14,5	2,3
Navarra (61 entrevistas)
PP	17,8	15,6 – 20,4	16,5
PSOE	21,2	18,6 – 24,0	27,1
VOX	10,0	7,8 – 12,7	5,6
SUMAR y aliados	8,4	6,4 – 10,6	12,7
BILDU	20,1	13,2 – 27,2	17,1
OTROS	22,4	17,9 – 27,5	20,9
Galicia (176 entrevistas)
PP	32,2	30,0 – 34,5	43,2
PSOE	22,5	20,0 – 25,1	29,6
VOX	12,7	10,2 – 14,9	4,7
SUMAR y aliados	6,4	4,6 – 8,1	10,8
BNG	9,9	7,3 – 13,7	9,4
OTROS	16,3	13,4 – 19,9	2,3
Canarias (148 entrevistas)
PP	26,9	23,8 – 30,0	30,1
PSOE	25,2	22,4 – 27,9	33,0
VOX	16,7	13,6 – 20,6	7,5
SUMAR y aliados	8,1	6,2 – 10,2	10,4
CCA	4,8	2,2 – 8,5	11,2
OTROS	18,3	14,2 – 22,1	7,8
23-J: porcentaje sobre votos emitidos (OTROS incluye blanco y nulo).

Algunas cosas que se ven:

Los partidos regionales salen cerca del 23-J donde hay datos. EH Bildu en el País Vasco (24,1 % frente a 23,8 %), BNG en Galicia (9,9 % frente a 9,4 %) y ERC en Cataluña (11,4 % frente a 13,0 %).
Otros salen bastante por debajo: PNV (19,3 % frente a 23,9 %), Junts (6,4 % frente a 11,1 %) y Coalición Canaria (4,8 % frente a 11,2 %). No sé si es un cambio real o un límite del modelo. Con 19 votantes del PNV o 6 de CC en la encuesta, el modelo tiene poco con lo que trabajar.
VOX y OTROS salen por encima del 23-J en todas. Es la tendencia nacional de la encuesta. El modelo solo tiene un efecto de comunidad por partido, así que el cambio respecto a 2023 se reparte de forma parecida en todas las comunidades.
Los intervalos reflejan los datos. En Navarra, con 61 entrevistas, EH Bildu va del 13 % al 27 %. En Cataluña, con 458, los intervalos son mucho más estrechos.

Jugando con la posterior

Hasta aquí he resumido cada posterior con su media y un intervalo. Pero tengo mucho más: 500 escenarios completos y coherentes entre sí. En cada draw, los once partidos suman 100 %. Así que cualquier pregunta sobre combinaciones de partidos se responde contando en cuántos draws se cumple.

Por ejemplo, el bloque PP + VOX. Sumo los dos partidos draw a draw y miro su distribución:

Mostrar código

pp_vox <- 100 * (mrp_draws[, "PP"] + mrp_draws[, "VOX"])

c(
  "PP + VOX > 45 %" = mean(pp_vox > 45),
  "PP + VOX > 50 %" = mean(pp_vox > 50),
  "PSOE > PP" = mean(mrp_draws[, "PSOE"] > mrp_draws[, "PP"]),
  "VOX > 20 %" = mean(mrp_draws[, "VOX"] > 0.20)
)
#> PP + VOX > 45 % PP + VOX > 50 %       PSOE > PP      VOX > 20 % 
#>           0.104           0.000           0.026           0.002

Mostrar código

tibble(pp_vox = pp_vox) |>
  ggplot(aes(x = pp_vox)) +
  geom_histogram(bins = 40, fill = "steelblue", color = "white") +
  geom_vline(xintercept = c(45, 50), linetype = "dashed", color = "firebrick") +
  labs(
    x = "PP + VOX (% sobre voto válido)",
    y = "Draws",
    title = "Distribución posterior del bloque PP + VOX",
    subtitle = "500 draws del MRP. Líneas en el 45 % y el 50 %"
  ) +
  theme_minimal()

El bloque PP + VOX se mueve entre el 41 % y el 47 % en los 500 draws, con media en el 43,6 %. La probabilidad de que pase del 45 % es del 10 %. La de que pase del 50 % es 0: ningún draw llega. La de que el PSOE quede por delante del PP es del 3 %.

Fíjate en que esto no se puede hacer con los intervalos de cada partido por separado. El intervalo de PP + VOX no es la suma de los dos intervalos, porque cuando el PP sube, VOX suele bajar: la correlación entre sus draws es de −0,37. En cada draw esa relación ya está dentro.

Lo mismo sirve por comunidad. Con los draws de cada comunidad puedo preguntar quién gana cada duelo:

Mostrar código

duelo <- function(cc, a, b) {
  mean(mrp_ccaa_draws[[cc]][, a] > mrp_ccaa_draws[[cc]][, b])
}

tribble(
  ~comunidad, ~a, ~b,
  "Cataluña", "ERC", "JUNTS",
  "Cataluña", "PSOE", "ERC",
  "País Vasco", "PNV", "BILDU",
  "Navarra", "BILDU", "PP",
  "Galicia", "BNG", "PSOE",
  "Canarias", "CCA", "VOX"
) |>
  mutate(probabilidad = pmap_dbl(list(comunidad, a, b), duelo)) |>
  gt(groupname_col = "comunidad") |>
  cols_merge(c(a, b), pattern = "{1} > {2}") |>
  cols_label(a = "Pregunta", probabilidad = "Probabilidad") |>
  fmt_percent(probabilidad, decimals = 0, locale = "es")

Pregunta	Probabilidad
Cataluña
ERC > JUNTS	99%
PSOE > ERC	100%
País Vasco
PNV > BILDU	19%
Navarra
BILDU > PP	66%
Galicia
BNG > PSOE	0%
Canarias
CCA > VOX	0%

En Cataluña, ERC queda por delante de Junts en el 99 % de los draws, y el PSOE por delante de ERC en todos. En el País Vasco, el MRP ve a EH Bildu por delante del PNV: el PNV solo gana en el 19 % de los draws. En Navarra, Bildu supera al PP en el 66 %, que es casi tirar una moneda.

Antes de sacar conclusiones

Cuatro avisos:

La estimación publicada por el CIS no es solo ponderar. Lleva su propia “cocina”: cómo reparte a los indecisos, qué hace con quien no contesta, etc. Aquí no hay nada de eso, porque todo se calcula sobre voto válido declarado.
No hay elecciones con las que comparar, así que no puedo decir qué estimación acierta más. Lo que sí puedo decir es que el MRP usa más información de la población, el recuerdo de voto, y que esa información corrige una desviación grande de la muestra.
Todo el electorado vota. El modelo estima la intención entre quienes declaran un voto válido, y la postestratificación la aplica a todo el electorado, abstencionistas incluidos. Es decir, supone que en cada celda vota todo el mundo, y con la misma intención que quien contesta. Por eso los abstencionistas de 2023 pesan tanto y suben OTROS. Una estimación electoral en serio necesitaría además un modelo de participación: la probabilidad de que cada celda vaya a votar.
Las probabilidades de la posterior son optimistas. Solo recogen la incertidumbre del modelo: tratan la tabla del raking como si fuera exacta, y no incluyen el sesgo de la propia encuesta ni lo que harán los indecisos. Además son porcentajes de voto, no escaños: pasar del 50 % de los votos no es lo mismo que tener mayoría absoluta. Y con 500 draws, una probabilidad de 0 % quiere decir “menos de 1 entre 500”, no “imposible”.

Pues hasta otra, y que ustedes voten libres, cuando se pueda.