Plato, cocina, receta, ingredientes #2

Receta para un plato mixto

Twitter: @nelsonamayad
Publicación: 25-04-2018
Última actualización: 19-05-2018

“Linear regression is the geocentric model of applied statistics.”
- Richard McElreath

La primera receta para hacerle seguimiento a las encuestas logra mucho con muy poco. Pero desaprovecha la poca información que hay.

Esta segunda receta le saca más jugo a lo poco que hay. Requiere más tiempo en la cocina y usa nuevos ingredientes, pero el resultado lo justifica.

El truco del plato mixto es explotar relaciones entre las características de las encuestas; por ejemplo, que las encuestas más pequeñas son, a la vez, las que se hacen por teléfono. En términos estadísticos, esta receta es un modelo mixto que extiende los efectos aleatorios por encuestadora que tenía la primera receta (\(\alpha_{encuestadora[i]}\)) a los coeficientes de las otras variables explicativas (\(\beta_{encuestadora[i]}\)).

¿Qué se logra con eso? Un mejor plato. ¿Por qué? No tengo la menor idea sobre cómo decir eso en español, o en metáfora culinaria, pero la idea es que cuando hay clusters que no tienen un órden determinado en los datos, como en este caso son las encuestadoras, se puede modelar simultáneamente lo que está dentro esos clusters y usar esas inferencias para agrupar (pool) información entre parámetros.

Siéntase libre de leer por encima, o ignorar, lo que le parezca demasiado técnico. No se pierde de mucho. De cualquier forma, la apuesta de este ejercicio no cambia: la transparencia es algo que se hace, no algo que se dice. Aquí no hay salsa secreta.

Buen provecho

---

Plato

Primero veamos cómo se comparan los resultados del segundo plato mixto con los del primer plato simple. En ambos casos el resultado es el promedio de la estimación de los parámetros y el HPDI de 90% de cada uno.

Los valores cambian con cada estimación. Eso refleja la incorporación de nueva información y la forma como la receta los procesa.

---

Ingredientes

Esta receta utiliza la misma información que la anterior: las encuestas que han salido después de las elecciones legislativas del 11 de marzo.

Encuestas

Las encuestas se importan 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

Este modelo comienza con los mismos priors del plato simple:

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

Pero es no es todo. Como este plato modela la relación entre todas las caraterísticas de las encuestas, necesita priors para las varianzas y correlaciones entre todos los parámetros.

Advertencia: lo que viene es técnico pero importante

El plato mixto utiliza un prior específico con la distribución Lewandowski-Kurowicka-Joe (LKJ para sus amigos), que permite definir con un sólo parámetro \(\eta\) la diagnonal de la matriz de correlación.

Utilizando el prior \(\eta=2\) las correlaciones están débilmente informadas y son escépticas de niveles extremos como -1 (perfecta correlación negativa) o 1 (perfecta correlación positiva). Ese es el prior que recomiendan el manual de Stan y texto McElreath para situaciones en las que esas relaciones no se conocen muy bien.

No voy a decir más sobre esto excepto que este truco permite lo que el plato mixto busca hacer: aprovechar mejor la poca información que hay.

---

Receta

Esta receta tiene la misma base que la anterior. La proporción de votos para cada candidato sale de modelo lineal sobre características observables 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).

La diferencia entre esta receta y la anterior es añadir (\(\beta_{encuestadora[i]}\)) a los coeficientes de las otras variables explicativas. Suena simple pero no lo es: implica incluir una matriz de varianza/covarianza que describe cómo se relaciona la distribución posterior de cada parámetro con los otros (ver capítulo 13 de Statistical Rethinking).

---

Preparación

Este es el modelo con los respectivos priors para cada parámetro. Nuevamente, los únicos priors informados son los que determinan el parámetro que captura el promedio y desviación estándar de cada candidato. El otro prior clave es el que va en la matriz de correlación LKJ.

Nunca fui muy hábil en álgebra lineal, así que esto va con uno que otro detalle innecesario y de pronto un horror en notación:

\[\small \lambda_{candidato,t} \sim Normal(\mu_{candidato,t},\sigma_{candidato,t})\] \[\small \mu_{candidato,t} = \alpha_{encuestadora[i]}+\beta_{1,encuestadora[i]}*m+\beta_{2,encuestadora[i]}*e+\beta_{3,encuestadora[i]}*d+\beta_{4,encuestadora[i]}*tipo \] \[\left[\begin{array} {rrr} \alpha_{encuestadora[i]}\\ \beta_{1,encuestadora[i]}\\ \beta_{2,encuestadora[i]}\\ \beta_{3,encuestadora[i]}\\ \beta_{4,encuestadora[i]}\\ \end{array}\right] = MVNormal \Bigg( \left[\begin{array} {rrr} \alpha_{encuestadora[i]}\\ \beta_{1,encuestadora[i]}\\ \beta_{2,encuestadora[i]}\\ \beta_{3,encuestadora[i]}\\ \beta_{4,encuestadora[i]}\\ \end{array}\right],S \Bigg)\]

\[ \mathbf{S} = \mathbf{x} \cdot \mathbf{R} \cdot \mathbf{x} \] \[ \mathbf{x} = \left(\begin{array}{rrr} \sigma_{\alpha} & 0 & 0 & 0 & 0\\ 0 & \sigma_{\beta_1} & 0 & 0 & 0\\ 0 & 0 & \sigma_{\beta_2} & 0 & 0\\ 0 & 0 & 0 & \sigma_{\beta_3} & 0\\ 0 & 0 & 0 & 0 & \sigma_{\beta_4}\\ \end{array}\right) \] \[ \mathbf{R} \sim LKJcorr(2) \]

\[\small\alpha_{encuestadora} \sim Normal(\mu_{candidato},\sigma_{candidato}) \]

\[\small\beta_1,\beta_2,\beta_3,\beta_4 \sim Normal(0,10) \]

\[\small\beta_{encuestadora} \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

Este paso es idéntico al del plato simple:

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

El plato mixto también se prepara con RStan, y en este caso voy a usar al candidato Iván Duque como ejemplo.

El código es miedoso, y tuve que verificarlo varias veces con el paquete rethinking de McElreath. Pero logra su cometido.

El modelo utiliza el data frame id, así como los priors \(\small\mu_{candidato}=39\), \(\small\sigma_{candidato}=4\) y para la matriz de correlaciones LKJ \(\small\eta=2\). Toda la receta va en el siguiente objeto duque.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{
    vector[N_encuestadora] b4_encuestadora;
    vector[N_encuestadora] b3_encuestadora;
    vector[N_encuestadora] b2_encuestadora;
    vector[N_encuestadora] b1_encuestadora;
    vector[N_encuestadora] a_encuestadora;
    real a;
    real b1;
    real b2;
    real b3;
    real b4;
    vector<lower=0>[5] s_encuestadora;
    real<lower=0> s;
    corr_matrix[5] Rho;
}
transformed parameters{
    vector[5] v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[N_encuestadora];
    vector[5] Mu_ab1b2b3b4;
    cov_matrix[5] SRS_s_encuestadoraRho;
    for ( j in 1:N_encuestadora ) {
        v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[j,1] = a_encuestadora[j];
        v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[j,2] = b1_encuestadora[j];
        v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[j,3] = b2_encuestadora[j];
        v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[j,4] = b3_encuestadora[j];
        v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora[j,5] = b4_encuestadora[j];
    }
    for ( j in 1:5 ) {
        Mu_ab1b2b3b4[1] = a;
        Mu_ab1b2b3b4[2] = b1;
        Mu_ab1b2b3b4[3] = b2;
        Mu_ab1b2b3b4[4] = b3;
        Mu_ab1b2b3b4[5] = b4;
    }
    SRS_s_encuestadoraRho = quad_form_diag(Rho,s_encuestadora);
}
model{
    vector[N] m;
    Rho ~ lkj_corr( 2 );
    s ~ cauchy( 0 , 5 );
    s_encuestadora ~ cauchy( 0 , 5 );
    b4 ~ normal( 0 , 10 );
    b3 ~ normal( 0 , 10 );
    b2 ~ normal( 0 , 10 );
    b1 ~ normal( 0 , 10 );
    a ~ normal( 39 , 4 );
    v_a_encuestadorab1_encuestadorab2_encuestadorab3_encuestadorab4_encuestadora ~ multi_normal( Mu_ab1b2b3b4 , SRS_s_encuestadoraRho );
    for ( i in 1:N ) {
        m[i] = a_encuestadora[encuestadora[i]] + b1_encuestadora[encuestadora[i]] *      muestra_int_voto[i] + b2_encuestadora[encuestadora[i]] * m_error[i] +      b3_encuestadora[encuestadora[i]] * dd[i] + b4_encuestadora[encuestadora[i]] *      tipo[i];
    }
    int_voto ~ normal( m , s );
}
generated quantities{
    vector[N] m;
    real dev;
    dev = 0;
    for ( i in 1:N ) {
        m[i] = a_encuestadora[encuestadora[i]] + b1_encuestadora[encuestadora[i]] *      muestra_int_voto[i] + b2_encuestadora[encuestadora[i]] * m_error[i] +      b3_encuestadora[encuestadora[i]] * dd[i] + b4_encuestadora[encuestadora[i]] *      tipo[i];
    }
    dev = dev + (-2)*normal_lpdf( int_voto | m , s );
}

Ahora nos vamos a RStan para prepara el plato:

library(rstan)
options(mc.cores = parallel::detectCores())
duque_fit <- stan(file='duque.stan',data=list(N=17,N_encuestadora=6,int_voto=id$int_voto,encuestadora=id$encuestadora, muestra_int_voto=id$muestra_int_voto,m_error=id$m_error,dd=id$dd,tipo=id$tipo),control=list(adapt_delta=0.95),iter=4000,chains=4)

## Warning: There were 256 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 2234 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: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low

## Warning: Examine the pairs() plot to diagnose sampling problems

El plato mixto necesita mucho más tiempo para cocinarse que el simple, pero al final sale. Aunque duque.stan tiene muchas divergencias iniciales en el muestreo, el plato sale del horno sin quemarse.

El resultado es un plato mixto que sabe mejor que el plato simple. Puede que el sabor sea demasiado bueno (i.e. overfitting).

Para los que tuvieron la paciencia de bajar hasta acá, va un premio: los valores observados para Iván Duque en las encuestas (en naranja) vs. los valores estimados por el plato mixto (en rojo).

Este plato tiene muchísimos parámetros, así que abajo sólo presento los más relevantes:

library(bayesplot)
posterior <- as.matrix(duque_fit)
color_scheme_set("orange")
mcmc_areas(posterior,prob=0.9,prob_outer = 0.99,point_est="mean",pars=c("m[1]","m[2]","m[3]","m[4]","m[5]","m[6]","m[7]","m[8]","m[9]","m[10]","m[11]","m[12]","m[13]","m[14]","m[15]","m[16]","m[17]"))

Para todo lo demás shinystan.

library(shinystan)

## Loading required package: shiny

## 
## This is shinystan version 2.5.0

launch_shinystan(duque_fit)

## 
## Launching ShinyStan interface... for large models this  may take some time.

## 
## Listening on http://127.0.0.1:6075

---

Referencias

McElreath, R. (2015). Statistical Rethinking. Texts in Statistical Science. Capítulo 13.

McElreath, R. “Multilevel Regression as Default”

Gelman, A. “Multilevel (Hierarchical) Modeling: What It Can and Cannot Do”.

Entrada del blog stlaPblog sobre LKJ en Stan

Lewandowski, D., Kurowicka, D., & Joe, H. (2009). Generating random correlation matrices based on vines and extended onion method. Journal of multivariate analysis, 100(9), 1989-2001.

Stan Development Team. 2017. RStan: the R interface to Stan. R package version 2.16.2. http://mc-stan.org

Stan Development Team. 2017. ShinyStan: Interactive Visual and Numerical Diagnostics and Posterior Analysis for Bayesian Models. R package version 2.4.0. http://mc-stan.org

Silver, N. (2017) “The Media has a probability problem”. 538.