14  🍲 Plato Mixto

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

📅 Publicado 25 Abril 2018, Actualizado 10 Febrero 2022

El Plato Simple logra mucho con muy poco. Pero desaprovecha la poca información que hay. El Plato Mixto le saca más jugo a lo 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 y no presencialmente. 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 de esos clusters y usar esas inferencias para agrupar (pool) información entre los 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.


Platos servidos

Primero veamos cómo se comparan los resultados de Plato Mixto con los del 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. Y los mismos priors que el Plato Simple.

Priors por candidato
Candidato \(\mu_{prior}\) \(\sigma_{prior}\)
Ivan Duque 38.1 3.3
Gustavo Petro 26.8 3.4
Sergio Fajardo 13.1 2.9
German Vargas Lleras 7.5 2.1
Humberto de la Calle 3.0 1.1
1 Encuestas disponibles: 18

Priors Plato Mixto


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.

Ojo, lo que sigue es algo técnico

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 el 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{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) \]


Datos

Este paso es idéntico al Plato Simple.

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 los priors \(\small\mu_{candidato}=38\), \(\small\sigma_{candidato}=3\) y para la matriz de correlaciones LKJ \(\small\eta=2\). Toda la receta va en el siguiente objeto duque.stan:

Code
data{
    vector[18] int_voto;
    int tipo[18];
    vector[18] dd;
    vector[18] m_error;
    int muestra_int_voto[18];
    int encuestadora[18];
}
parameters{
    vector[7] b4_encuestadora;
    vector[7] b3_encuestadora;
    vector[7] b2_encuestadora;
    vector[7] b1_encuestadora;
    vector[7] a_encuestadora;
    real a;
    real b1;
    real b2;
    real b3;
    real b4;
    vector<lower=0>[5] s_encuestadora;
    real s;
    corr_matrix[5] Rho;
}
model{
    vector[18] 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( 38 , 4 );
    {
    vector[5] YY[7];
    vector[5] MU;
    MU = [ a , b1 , b2 , b3 , b4 ]';
    for ( j in 1:7 ) YY[j] = [ a_encuestadora[j] , b1_encuestadora[j] , b2_encuestadora[j] , b3_encuestadora[j] , b4_encuestadora[j] ]';
    YY ~ multi_normal( MU , quad_form_diag(Rho , s_encuestadora) );
    }
    for ( i in 1:18 ) {
        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[18] m;
    real dev;
    dev = 0;
    for ( i in 1:18 ) {
        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:

Code
library(rstan)

options(mc.cores = parallel::detectCores())

duque_fit <- stan(file='duque.stan',
                  data=list(N=18,
                            N_encuestadora=7,
                            int_voto=id_2018_simple$int_voto,
                            encuestadora=id_2018_simple$encuestadora,
                            muestra_int_voto=id_2018_simple$muestra_int_voto,
                            m_error=id_2018_simple$m_error,
                            dd=id_2018_simple$dd,
                            tipo=id_2018_simple$tipo),
                  control=list(adapt_delta=0.95),
                  iter=4000,
                  chains=4,
                  cores = 4)


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

Code
library(bayesplot)
color_scheme_set("orange")
  
# Crear matriz de muestras de la distribucion posterior
posterior_mixto <- as.matrix(duque_fit)
  
#Trace plot para los parametros de interes
mcmc_trace(posterior_mixto,pars=vars(starts_with("m[")))


El resultado es un plato que sabe mejor que el 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:

Code
library(bayesplot)

posterior_mixto <- as.matrix(duque_fit)

color_scheme_set("orange")

mcmc_areas(posterior_mixto,
           prob=0.9,
           prob_outer = 0.99,
           point_est="mean",
           pars=vars(starts_with("m[")))


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.