Twitter: @nelsonamayad
Publicación: 10-06-2018

“This self-deceit, this fatal weakness of mankind, is the source of half the disorders of human life”
- Adam Smith


Llegó la hora de ver si la presidencia de Colombia en 2018 queda en cabeza de @IvanDuque o @GustavoPetro. Ya se hizo la tarea para la primera vuelta y los resultados fueron muy buenos, a pesar de que no capturaron la sorpresa fajardista que le sacó varios puntos a Petro y a De la Calle en la última semana.

Este postre busca lograr lo que los anteriores no pudieron: estimar qué tan probable es que gane uno u otro candidato. Para hacer eso el postre le aplica una transformación sencillita al plato Simple: en vez de servir dos recetas por candidato, se sirve una sola logística para ambos. Esto permite generar un pronóstico probabilístico para la segunda vuelta.

¿Por qué hacerlo así? Porque según las encuestas la victoria es para Duque: todas coinciden en que el candidato de la coalición de derecha obtiene una mayor proporción de votos que su contrincante. Es más, el resultado del plato mixto, al que mejor le fue en la primera vuelta, aplicado a las encuestas que salieron después del 27 de mayo da como resultado que Duque obtendrá 51% de la votación, con un rango intervalo HPD de 95% de entre 46% y 57%, y Petro obtendría 37%, con un rango HPD de 95% entre 34% y 41%. Así que en vez de decir lo obvio, nos vamos por lo menos evidente, aunque sea más difícil de masticar.

Buen provecho.

Postre Logístico


Este postre recoge la probabilidad de que @IvanDuque sea presidente, estimada por el modelo que se describe abajo.

Estas probabilidades salen de sacar muchas simulaciones de cada uno de los parámetros del modelo estimado. El resultado es una densidad de la probabilidad de que Duque obtenga más de 50% de la votación. La distribución parece bimodal porque unas encuestas le dan menos de 50% de la votación a Duque, pero el grueso dice lo contrario: el 20% de las simulaciones estima que Duque gana con más o menos 53% de probabilidad.

@IvanDuque tiene más o menos 65% de probabilidad de obtener más del 50% de la votación.


library(ggplot2)
library(RCurl)

#Importar simulaciones
id.4cast <- read.csv(text=getURL("https://raw.githubusercontent.com/nelsonamayad/nelsonamayad.github.io/master/logis/id4cast.csv"))

#postre logis
ggplot(id.4cast)+
  geom_density(aes(x=id4cast,fill="orangered3",color="grey60"))+
  geom_vline(xintercept=0.5,linetype="dotted")+
labs(x="Probabilidad de que Duque obtenga más de 50% de la votación",y="Densidad Probabilidades estimadas",
       title="Postre logístico", 
       subtitle="10.000 simulaciones del modelo logístico para la 2da vuelta", 
       caption = "Fuente: Cálculos @nelsonamayad")+
  theme(legend.position = "none",
        panel.background = element_blank())+
  xlim(0.35,0.65)+
  scale_fill_brewer(palette = "Spectral")
## Warning: Removed 6935 rows containing non-finite values (stat_density).

Ingredientes


Los únicos ingredientes de esta receta son las encuestas que han salido desde la primera vuelta. Como priors se toman los promedios y desviaciones estándar de cada candidato. 

library(tidyverse)
library(RCurl)
library(lubridate)

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

# 2. Alistamiento de los datos
d <- d %>% filter(n>45) %>% select(encuestadora,ivan_duque,gustavo_petro,fecha,muestra_int_voto,margen_error,tipo,municipios)
d <- d %>% mutate(id = ivan_duque*muestra_int_voto/100, gp = gustavo_petro*muestra_int_voto/100,n=muestra_int_voto,m_error=margen_error,tipo = ifelse(d$tipo=="Presencial",1,0))
d <- d %>% mutate(id = round(id,digits=0),gp = round(gp,digits=0))
f <- today()    #FECHA ESPECIFICA
d <- d %>% mutate(dd = as.Date(as.character(f), format="%Y-%m-%d") - 
                    as.Date(as.character(fecha), format="%Y-%m-%d"))
d <- d %>% mutate(dd = as.numeric(dd))
d$encuestadora <- match(d$encuestadora, unique(d$encuestadora))

Receta


Este postre modela directamente las encuestas como una distribución binomial donde cada encuesta es un ensayo:

\[ N^{Duque}_i \sim \textrm{Binomial}(N_i, \pi^{Duque}_i) \] donde \(N_i\) es la muestra de cada encuesta y \(\pi^{Duque}_i\) es la proporción de la intención de voto para Iván Duque.

Ahora, la proporción de intención de voto para Duque \(\pi^{Duque}_i\) se determina a través de una función logística (log-link function) y, como es costumbre en este recetario, efectos aleatorios por encuestadora así:

\[\textrm{logit} (\pi^{Duque}_i) = \alpha_{encuestadora}\]

Preparación


Este es el postre completo, con los priors del caso:
\[ N^{Duque}_i \sim \textrm{Binomial}(N_i, \pi^{Duque}_i) \] \[\textrm{logit} (\pi^{Duque}_i) = \alpha_{encuestadora}\] \[\small\alpha_{encuestadora} \sim Normal(\alpha,\sigma) \] \[\small\alpha \sim Normal(50, 5) \] \[\small\sigma \sim HalfCauchy(0,5) \]

Estimación en RStan


Este es el modelo en el siempre espeluznante código de RStan: 

data{
    int<lower=1> N;
    int<lower=1> N_encuestadora;
    int id[N];
    int n[N];
    int encuestadora[N];
}
parameters{
    vector[N_encuestadora] a;
    real<lower=0> s;
}
model{
    vector[N] p;
    s ~ cauchy( 0 , 5 );
    a ~ normal( 0 , 10 );
    a ~ normal( a , s );
    for ( i in 1:N ) {
        p[i] = a[encuestadora[i]];
    }
    id ~ binomial_logit( n , p );
}
generated quantities{
    vector[N] p;
    real dev;
    dev = 0;
    for ( i in 1:N ) {
        p[i] = a[encuestadora[i]];
    }
    dev = dev + (-2)*binomial_logit_lpmf( id | n , p );
}


Ahora meter el postre al horno de RStan: 

library(rstan)
options(mc.cores = parallel::detectCores())
logis_fit <- stan(file='logis.stan',data=list(N=7,N_encuestadora=6,n=d$n,encuestadora=d$encuestadora,id=d$id),control=list(adapt_delta=0.95),iter=4000,chains=4)
## Warning: There were 7198 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: There were 339 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems


Con algo de trabajo, y a pesar de muchas divergencias iniciales en el muestreo, el postre no se quema y sale del horno.

Veámos cómo se comparan los parametros estimados con los datos observados: 

library(tidyverse)

#Muestreo de la distribucion posterior
posterior <- as.matrix(logis_fit)

# Parametros vs Observados
ps <- posterior %>% as.tibble() %>% select(7:13) %>% summarize_all(funs(mean)) %>% t() %>% as.data.frame()

# Convertir parametros con funcion inversa logistica
inv_log <- function(x) {
  1/(1+exp(-x))
}
ps <- ps %>% mutate_all(funs(inv_log))

#Valores observados
obs <- d %>% mutate(id = id/n,t=seq(from=1,to=7))

#Bind
ps_obs <- cbind(ps,obs)

#Grafica
ggplot(ps_obs, aes(x=t))+
    geom_point(aes(y=id),size=4,color="orangered")+
    geom_point(aes(y=ps),shape=10,size=5,color="red2")+
    theme_classic()+labs(y="intencion de voto",x="")+
  ylim(0.3,0.7)


Para terminar, vamos las densidades de los parámetros estimados usando bayesplot y lanzar shinystan: 

#Grafica
library(bayesplot)
posterior <- as.matrix(logis_fit)
color_scheme_set("orange")
mcmc_areas(posterior,prob=0.95,prob_outer = 0.99,point_est="mean",pars=c("p[1]","p[2]","p[3]","p[4]","p[5]","p[6]","p[7]"))

#Shinystan
library(shinystan)
launch_shinystan(logis_fit)


Bonus: Estimación con map2stan


Además de la siempre abstrusa representación del código en RStan, en este postre incluyo una alternativa que de paso le hace bombo al paquete rethinking que preparó uno de los mejores cocineros: Richard McElreath.

La especificación del modelo en el código del paquete rethinking, que traduce a un lenguaje más amable el código de RStan, se ve así: 

library(rethinking)

# Modelo logístico en map2stan:
logis <- map2stan(
  alist(
    #modelo
    id ~ dbinom(n,p),
    logit(p) <- a[encuestadora],
    #priors
    a[encuestadora] ~ dnorm(a,s),
    a ~ dnorm(50,5),
    s ~ dcauchy(0,5)
     ),
  data=d,
  control=list(adapt_delta=0.96),
  iter=4000, warmup=1000, chains=4, cores=2)
## Warning: There were 11225 divergent transitions after warmup. Increasing adapt_delta above 0.96 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning in map2stan(alist(id ~ dbinom(n, p), logit(p) <- a[encuestadora], : There were 11225 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.


Referencias


Este proyecto de Pierre-Antoine Kremp, para las presidenciales en EEUU de 2016, fue la inspiración para estas recetas.

McElreath, R. (2015). Statistical Rethinking. Texts in Statistical Science. Bendito sea Richard McElreath por este texto.

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

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.

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