10.9 Diagnósticos generales

Cuando generamos una muestra de la distribución posterior usando MCMC, sin importar el método (Metrópolis, Gibbs, HMC), buscamos que:

Los valores simulados sean representativos de la distribución posterior, esto implica que no deben estar influenciados por el valor inicial (arbitrario) y deben explorar todo el rango de la posterior.
Debemos tener suficientes simulaciones de tal manera que las estimaciones sean precisas y estables.
Queremos tener un método eficiente para generar las simulaciones.

En la práctica intentamos cumplir lo más posible estos objetivos, pues aunque en principio los métodos MCMC garantizan que cadena infinitamente largas lograran una representación perfecta, siempre debemos tener un criterio para cortar la cadena y evaluar la calidad de las simulaciones.

Representatividad

Para determinar la convergencia de la cadena es conveniente realizar más de una cadena, buscamos ver si realmente se ha olvidado el estado inicial y revisar que algunas cadenas no hayan quedado atoradas en regiones inusuales del espacio de parámetros.

set.seed(337687)

norm_fit_1 <- sampling(object = stan_norm_cpp, data = data_list,
  chains = 3 , iter = 40 , warmup = 20)
#> Warning: The largest R-hat is 1.14, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#tail-ess

norm_posterior_inc_warmup <- extract(norm_fit_1, inc_warmup = TRUE, 
  permuted = FALSE)

mcmc_trace(norm_posterior_inc_warmup,  pars = c("mu", "sigma2"), n_warmup = 20,
  facet_args = list(nrow = 2, labeller = label_parsed)) + 
  facet_text(size = 15)
#> Error in mcmc_trace(norm_posterior_inc_warmup, pars = c("mu", "sigma2"), : could not find function "mcmc_trace"

La gráfica anterior nos puede ayudar a determinar si elegimos un periodo de calentamiento adecuado o si alguna cadena está alejada del resto.

Además de realizar gráficas podemos usar la medida de convergencia \(\hat{R}\) que la función rstan() calcula por omisión:

norm_fit_1
#> Inference for Stan model: modelo_normal.
#> 3 chains, each with iter=40; warmup=20; thin=1; 
#> post-warmup draws per chain=20, total post-warmup draws=60.
#> 
#>          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
#> mu       1.84    0.05 0.31   1.18   1.73   1.88   2.00   2.41    40 1.04
#> sigma2   4.76    0.11 0.89   3.24   4.06   4.71   5.23   6.56    61 0.97
#> lp__   -71.06    0.20 1.24 -74.81 -71.56 -70.47 -70.30 -70.06    38 1.07
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Dec 10 23:25:35 2019.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

La medida \(\hat{R}\) se conoce como el factor de reducción potencial de escala o diagnóstico de convergencia de Gelman-Rubin, esta es una estimación de la posible reducción en la longitud de un intervalo de confianza si las simulaciones continuran infinitamente. \(\hat{R}\) es aproximadamente la raíz cuadrada de la varianza de todas las cadenas juntas dividida entre la varianza dentro de cada cadena. Si \(\hat{R}\) es mucho mayor a 1 esto indica que las cadenas no se han mezclado bien. Una regla usual es iterar hasta alcanzar un valor \(\hat{R} \leq 1.1\) para todos los parámetros.

\[\hat{R} = \sqrt{\frac{\hat{d}+3}{\hat{d}+1}\frac{\hat{V}}{W}}\]

donde \(B\) es la varianza entre las cadenas, \(W\) es la varianza dentro de las cadenas

\[B = \frac{N}{M-1}\sum_m (\hat{\theta_m} - \hat{\theta})^2\] \[W = \frac{1}{M}\sum_m \hat{\sigma_m}^2\] Y \(\hat{V}\) es una estimación del varianza de posteriro de \(\theta\):

\[\hat{V} = \frac{N-1}{N}W + \frac{M+1}{MN}B\]

Precisión

Una vez que tenemos una muestra representativa de la distribución posterior, nuestro objetivo es asegurarnos de que la muestra es lo suficientemente grande para producir estimaciones estables y precisas de la distribución.

Para ello usaremos otra salida de la función jags: n.eff que es el tamaño efectivo de muestra, si las simulaciones fueran independientes n.eff sería el número total de simulaciones; sin embargo, las simulaciones de MCMC suelen estar correlacionadas, el tamaño efectivo nos dice que tamaño de muestra de observaciones independientes nos daría la misma información que las simulaciones de la cadena.

\[NEM = \frac{N}{1+2\sum_{k=1}^\infty ACF(k)} \]

Usualmente nos gustaría obtener un tamaño efectivo de al menos \(100\).

norm_fit
#> Inference for Stan model: modelo_normal.
#> 3 chains, each with iter=1000; warmup=500; thin=1; 
#> post-warmup draws per chain=500, total post-warmup draws=1500.
#> 
#>          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
#> mu       1.91    0.01 0.31   1.30   1.70   1.92   2.12   2.53  1071 1.00
#> sigma2   4.72    0.03 0.92   3.26   4.07   4.62   5.26   6.81   981 1.00
#> lp__   -71.04    0.05 0.99 -73.64 -71.39 -70.76 -70.33 -70.07   455 1.01
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Dec 10 23:25:34 2019.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

Eficiencia

Hay varias maneras para mejorar la eficiencia de un proceso MCMC:

Paralelizar, no disminuimos el número de pasos en las simulaciones pero podemos disminuir el tiempo que tarda en correr.
Cambiar la parametrización del modelo o transformar los datos.
Adelgazar la muestra: esto nos puede ayudar a disminuir el uso de memoria, consiste en guardar únicamente los \(k\)-ésimos pasos de la cadena y resulta en cadenas con menos autocorrelación (en el caso de la función sampling() el adelgazamiento está controlado por el parámetro thin).

Recomendaciones generales

Gelman and Hill (2007) recomienda los siguientes pasos cuando uno esta simulando de la posterior:

Cuando definimos un modelo por primera vez establecemos un valor bajo para iter. La razón es que la mayor parte de las veces los modelos no funcionan a la primera por lo que sería pérdida de tiempo dejarlo correr mucho tiempo antes de descubrir el problema.
Si las simulaciones no han alcanzado convergencia aumentamos las iteraciones a \(500\) ó \(1000\) de tal forma que las corridas tarden segundos o unos cuantos minutos.
Si tarda más que unos cuantos minutos (para problemas del tamaño que veremos en la clase) y aún así no alcanza convergencia entonces juega un poco con el modelo (por ejemplo intenta transformaciones lineales), para JAGS Gelman sugiere más técnicas para acelerar la convergencia en el capitulo \(19\) del libro Data Analysis Using Regression and Multilevel/Hierarchical models. En el caso de Stan veremos ejemplos de reparametrización, y se puede leer más en la guía.
Otra técnica conveniente cuando se trabaja con bases de datos grandes (sobretodo para la parte exploratoria) es trabajar con un subconjunto de los datos, quizá la mitad o una quinta parte.

Referencias

Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Vol. Analytical methods for social research. New York: Cambridge University Press.