Twitter: @nelsonamayad
Publicación: 17-04-2018
Última actualización: 20-05-2018.

“Far better an approximate answer to the right question, which is often vague, than the exact answer to the wrong question, which can always be made precise.”
- John Tukey


Este es la primera de varias recetas. Todas tienen el mismo propósito: ver qué se puede inferir de las encuestas de intención de voto para las elecciones presidenciales de Colombia en 2018. Cada receta combina tres cosas chéveres: estadística bayesiana, reproducibilidad y computación open-source. Las tendencias en las encuestas se ven mejor acá.

La primera receta es simple y tiene varias limitaciones. Dos en especial: no modela la determinación simultánea de la votación para todos los candidatos, ni el sesgo de respuesta de las encuestas. De cualquier forma, acá se pone todo sobre la mesa. Las críticas constructivas, bienvenidas.

Ya que los resultados se pueden interpretar como un pronóstico electoral, esta página va en reversa. Primero el plato (los resultados), después los ingredientes (datos), luego la receta (modelo), y al final la preparación en la cocina (código). Entre más se baje la página, más técnico se pone este asunto.

Utiliza RStan y lo demás se hace con el tidyverse de R.

Buen provecho.

Plato


Estos son los resultados del modelo. Son promedios de la distribución posterior estimada para cada candidato, así como los intervalos HPD (higher posterior density) de 90% sobre esos parámetros. Para tener una referencia, el resultado se compara con un promedio simples sobre las 18 encuestas que se han hecho después de las elecciones legislativas del 11 de marzo.

Como escribí en la primera publicación, esperaba que conforme salieran más encuestas el intervalo se redujera. Estos intervalos se han apretado con cada estimación porque la varianza de la intención de voto para todos los candidatos ha caído.

Ingredientes


El modelo utiliza las encuestas que han salido hasta la última fecha de actualización.

Antes de las consultas del 11 de marzo, y de la adhesión de Juan Carlos Pinzón a la campaña de German Vargas Lleras (el 16 de marzo de 2016), las encuestas estaban identificando un conjunto ruidoso de candidatos. Por esa razón, este modelo solo tiene en cuenta las encuestas realizadas después de las elecciones legislativas.

Encuestas

Los datos básicos de las encuestas se pueden importar desde GitHub con RCurl.

library(RCurl)

# Importar encuestas desde GitHub:
encuestas <- read.csv(text=getURL("https://raw.githubusercontent.com/nelsonamayad/Elecciones-presidenciales-2018/master/Elecciones%202018/encuestas2018.csv"))

Priors

El modelo utiliza como priors para la estimación de un parámetro la proporción de votos promedio \(\mu_{candidato}\) y la desviación estándar \(\sigma_{candidato}\) que han registrado las encuestas para cada candidato. Los demás parámetros tienen priors poco informativos.

library(tidyverse)
library(kableExtra)

# Preparar data frame para calcular priors
priors <- encuestas %>% select(n,fecha,ivan_duque,gustavo_petro,sergio_fajardo,german_vargas_lleras,humberto_delacalle) %>% filter(as.Date(fecha)>as.Date("2018-03-11"))
priors <- priors %>% select(-n,-fecha) %>% gather(candidato, int_voto)

# Calcular priors
priors <- priors %>% group_by(candidato) %>% mutate(prior_mu=mean(int_voto),prior_sigma=sd(int_voto))
priors <- priors %>% distinct(prior_mu,prior_sigma) 

# Tabla de priors
priors %>% kable("html", digits=0,caption = "Priors por candidato") %>% kable_styling(full_width = F)
Priors por candidato
candidato prior_mu prior_sigma
ivan_duque 38 3
gustavo_petro 27 3
sergio_fajardo 13 3
german_vargas_lleras 8 2
humberto_delacalle 3 1

Receta


El modelo parte del supuesto de que la proporción de votos \(\pi\) que obtiene un candidato en las elecciones en el momento t es un reflejo de las preferencias que tiene la sociedad por ese candidato antes de las elecciones:

\[\pi_{candidato,t} \sim Normal(\pi_{candidato,t-1}, \sigma_{candidato,t-1})\] Como nadie es adivino para saber esas preferencias, solo se observan mediciones ruidosas de esa relación: las encuestas de intención de voto. Aunque no se puede conocer la proporción de votos que recibirá cada candidato antes del día de las elecciones, esa proporción es una función de la proporción de intención de voto \(\lambda\) que hayan capturado las encuestas que se hayan realizado antes de esa fecha.

\[\pi_{candidato,t-1} \sim Normal(\lambda_{candidato,t-1}, \sigma_{candidato,t-1})\] La proporción de votos para cada candidato se aproxima mediante un modelo lineal sobre las siguientes características de las encuestas: 1) el tamaño de la muestra de cada encuesta (m), 2) el márgen de error de la encuesta (e), 3) los días que pasaron entre la encuesta y la estimación (d), 4) una dummy para el tipo de encuesta (telefónica o presencial) (tipo). Además, se incluyen efectos aleatorios por encuestadora que permiten incorporar la variación a ese nivel.

Preparación


Este es el modelo completo, con los priors para cada parámetro. Los únicos priors informados son los que determinan el parámetro que captura el promedio y desviación estándar de cada candidato, y estos se actualizan con cada estimación del modelo cuando sale una nueva encuesta.

\[\small \lambda_{candidato,t} \sim Normal(\mu_{candidato,t},\sigma_{candidato,t})\] \[\small \mu_{candidato,t} = \alpha_{t}+\alpha_{encuestadora[i]}+\beta_1*m+\beta_2*e+\beta_3*d+\beta_4*tipo \] \[\small\alpha_{t} \sim Normal(\mu_{candidato},\sigma_{candidato}) \] \[\small\beta_1,\beta_2,\beta_3,\beta_4 \sim Normal(0,10) \] \[\small\alpha_{encuestadora[i]} \sim Normal(\mu, \sigma) \] \[\small\mu \sim Normal(0,10) \] \[\small\sigma \sim HalfCauchy(0,5) \] \[\small\sigma_{candidato} \sim HalfCauchy(0,5) \]

Alistamiento de los datos

Hay que hacer unos cuantos ajustes a los datos de las encuestas antes de estimar el modelo:

library(lubridate)

#1. Depurar encuestas:
df <- encuestas %>% select(n,fecha, encuestadora,ivan_duque,gustavo_petro,sergio_fajardo,german_vargas_lleras, humberto_delacalle, margen_error,tipo, muestra_int_voto)

#2. Solo las encuestas post 2018-03-11:
df <- df %>% filter(as.Date(fecha, tz="GMT") >= as.Date('2018-03-11', tz="GMT"))

#3. Crear variable duracion:
df <- df %>% mutate(dd = as.Date(as.character(today()), format="%Y-%m-%d") - as.Date(as.character(fecha), format="%Y-%m-%d"))
df <- df %>% mutate(dd = as.numeric(dd))
df <- df %>% mutate(dd = 100*(dd/sum(dd)))

#4. Codificar encuestadoras:
df <- df %>% mutate(enc=encuestadora)
df$encuestadora <- match(df$encuestadora, unique(df$encuestadora))

#5. Dummy tipo de encuesta (=1 si presencial):
df <- df %>% mutate(tipo=ifelse(tipo=="Presencial",1,0))

#6. Otros ajustes:
df <- df %>% rename(m_error = margen_error)

#7. Pasar a formato largo:
df <- df %>% select(-fecha,-n) %>% gather(candidato, int_voto, ivan_duque,gustavo_petro,sergio_fajardo,german_vargas_lleras,humberto_delacalle)

#8. Crear data frames por candidato:
id <- df %>% filter(candidato=="ivan_duque") 
gp <- df %>% filter(candidato=="gustavo_petro") 
sf <- df %>% filter(candidato=="sergio_fajardo") 
gvl <- df %>% filter(candidato=="german_vargas_lleras") 
hdlc <- df %>% filter(candidato=="humberto_delacalle")

Estimación

Este es el código para estimar el modelo para cada candidato. Solo se necesitan los datos cargados en R y tener el paquete RStan instalado (ver instrucciones acá)

El muestreo del modelo se hace en Stan, que para cada candidato utiliza su respectivo data frame y priors.

Por ejemplo, para el candidato Sergio Fajardo se utiliza el data frame sf y los priors \(\small\mu_{candidato}=12\) y \(\small\sigma_{candidato}=3\) en un objeto Stan de nombre fajardo.stan:

data{
    int<lower=1> N;
    int<lower=1> N_encuestadora;
    real int_voto[N];
    int encuestadora[N];
    real muestra_int_voto[N];
    real m_error[N];
    real dd[N];
    real tipo[N];
}
parameters{
    real a1;
    vector[N_encuestadora] a_;
    real a_enc;
    real<lower=0> s_enc;
    real a2;
    real a3;
    real a4;
    real a5;
    real<lower=0> s;
}
model{
    vector[N] m;
    s ~ cauchy( 0 , 5 );
    a5 ~ normal( 0 , 10 );
    a4 ~ normal( 0 , 10 );
    a3 ~ normal( 0 , 10 );
    a2 ~ normal( 0 , 10 );
    s_enc ~ cauchy( 0 , 5 );
    a_enc ~ normal( 0 , 10 );
    a_ ~ normal( a_enc , s_enc );
    a1 ~ normal( 12 , 3 );  //Priors Fajardo: mu=12, sd=3
    for ( i in 1:N ) {
        m[i] = a1+a_[encuestadora[i]]+a2*muestra_int_voto[i]+a3*m_error[i]+a4*dd[i]+a5*tipo[i];
    }
    int_voto ~ normal( m , s );
}
generated quantities{
    vector[N] m;
    real dev;
    dev = 0;
    for ( i in 1:N ) {
        m[i] = a1+a_[encuestadora[i]]+a2 * muestra_int_voto[i]+a3*m_error[i]+a4*dd[i]+a5*tipo[i];
    }
    dev = dev + (-2)*normal_lpdf( int_voto | m , s );
}


Ahora a RStan: 

library(rstan)
options(mc.cores = parallel::detectCores())
fajardo_fit <- stan(file='fajardo.stan',data=list(N=17,N_encuestadora=6,int_voto=sf$int_voto,encuestadora=sf$encuestadora, muestra_int_voto=sf$muestra_int_voto,m_error=sf$m_error,dd=sf$dd,tipo=sf$tipo),control=list(adapt_delta=0.95),iter=4000,chains=4)
## Warning: There were 11 divergent transitions after warmup. Increasing adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems


A pesar de las divergencias iniciales, el modelo converge rápido para todos los candidatos. Al fin y al cabo es muy simple y tiene pocas observaciones.

Por sugerencia de @infrahumano, incluyo dos gráficas antes de ir a shinystan: trace plot con bayesplot para ver cómo se comportaron las 4 cadenas en los parámetros clave, y una comparación a la carrera entre el promedio de cada parámetro y su valor observado.

Ahora un resumen de todos los parámetros del modelo:

library(bayesplot)
posterior <- as.matrix(fajardo_fit)
color_scheme_set("green")
mcmc_intervals(posterior,prob=0.9,prob_outer = 0.99,point_est="mean")



Para cerrar, se pueden inspeccionar otras características del modelo, diagnosticar otras partes del muestreo, ver los coeficientes y visualizar contrafactuales simulados con shinystan. 

library(shinystan)
launch_shinystan(fajardo_fit)

Algunas referencias en desorden


McElreath, R. (2015). Statistical Rethinking. Texts in Statistical Science.

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis. Boca Raton, FL: CRC press.

Stan Development Team (2016) Stan Modeling Language: User’s Guide and Reference Manual. Version 2.14.0.

Gelman, A., & Azari, J. (2017). 19 things we learned from the 2016 election. Statistics and Public Policy, 4(1), 1-10.

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper). Bayesian analysis, 1(3), pp.515-534.

Linzer, D. A. (2013). Dynamic Bayesian forecasting of presidential elections in the states. Journal of the American Statistical Association, 108(501), 124-134.

Bunker, K., & Bauchowitz, S. (2016). Electoral Forecasting and Public Opinion Tracking in Latin America: An Application to Chile. Politica, 54(2).

Shirani-Mehr, H., Rothschild, D., Goel, S., & Gelman, A. (2018). Disentangling bias and variance in election polls. Journal of the American Statistical Association, (just-accepted), 1-23.

Wickham, H., & Grolemund, G. (2016). R for data science: import, tidy, transform, visualize, and model data. O’Reilly Media, Inc.

Este proyecto de Pierre-Antoine Kremp para las presidenciales en EEUU de 2016

Esta entrada del blog de Jim Savage

Esta visualización del muestreo que hacen algoritmos como Marvok Chain Monte Carlo, Metropolis-Hastings y el de Stan: No-U-Turn sampler