sandboxodelamahttp://www.odelama.com/sandbox/odelamaikiwiki2020-12-31T20:12:34ZIntroducción a Regresiones Bayesianas Multinivelhttp://www.odelama.com/sandbox/test/Oscar de Lama2020-08-19T12:51:00Z2020-08-19T12:51:00Z
<div class="options" data-disqus="no" data-eocline="yes" data-mathjax="yes" data-title-icon="lab" data-uuid="e7d3f6e6-2aa8-4a71-a256-d9432b1119fd" style="display:none"></div>
<p>Frecuentemente los datos a analizar pueden ser agrupados por uno o más conceptos y es adecuado entender la data desde la perspectiva de dichos grupos. Ejemplos de grupos son el tipo de tratamiento efectuado al paciente, la región donde trabaja el empleado, o el colegio donde estudia un alumno.</p>
<p>Cuando se presenta este caso, una alternativa de análisis es hacer regresiones separadas para cada grupo, pero el problema entonces es que el modelo en cada grupo ignora por completo a las observaciones en los otros. Es decir, se ignora el patrón general de los datos, y hay un sobreajuste (<em>overfitting</em>) cuando hay pocas observaciones en un grupo.</p>
<p>El otro extremo es ignorar los grupos y hacer una sola regresión para todos los datos. Si se incluye el grupo, como un predictor de tipo categórico en la regresión, el resultado es similar al caso mencionado en el párrafo anterior. Si no se incluye el grupo, se ignora la tendencia al interior de cada uno de ellos y entonces hay un sub-ajuste (<em>underfitting</em>), sobre todo cuando hay marcadas diferencias en las observaciones entre grupos.</p>
<p>Las regresiones multinivel incorporan en un mismo modelo tanto la tendencia al interior de cada grupo como el patrón global de los datos. </p>
<p>El resultado de una regresión multinivel incluye la desviación estándar de los coeficientes <em>entre</em> los grupos, permitiendo entender cuando estos grupos son <em>imprescindibles</em> en el modelo. Sin embargo, en el caso en el que las observaciones son homogéneas entre los grupos, y éstos no son <em>imprescindibles</em>, la regresión multinivel <em>revierte automáticamente</em> a una regresión sin grupos.</p>
<p>Deberíamos usar regresiones multinivel como forma de análisis predeterminado.</p>
<p>Algunos estadísticos se refieren a los modelos multinivel con nombres diferentes, algunos refiriéndose a variantes especializadas, y otros usándolos indistintamente. Los sinónimos más comunes de “multinivel” son jerárquicos y de efectos mixtos. Algunos llaman a los coeficientes que varían por grupo "efectos aleatorios" a los globales "efectos fijos" y a la combinación de ellos "efectos mixtos".</p>
<div class="ContentTable">
<header class="tocHeader">Contents</header>
</div>
<h1 id="introduccin">Introducción</h1>
<p>En una regresión lineal clásica los coeficientes de la regresión son los mismos para todas las observaciones.</p>
<p>\begin{align}
y^{(i)} &= \beta_0 + \beta_1 x_1^{(i)} + \beta_2 x_2^{(i)} + \dots + \beta_p x_p^{(i)} +\epsilon^{(i)}, ~para~i=1, ..., N
\label{eq:clasLinReg} \\
\epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2) \
\label{eq:linRegEps}
\end{align}</p>
<p>En la ecuación \eqref{eq:clasLinReg} $y$ es un vector con los valores de la variable dependiente (llamada también variable de resultado), y la notación para referirse a su valor en la posición indicada por $i$ es $y^{(i)}$. Los vectores tienen forma vertical, son columnas. Las variables ${x_1^{(i)}, x_2^{(i)}, \dots}$ son las predictoras de $y^{(i)}$. En notación más compacta usamos</p>
<p>\begin{aligned}
y^{(i)} = {X}^{(i)} {\beta} + \epsilon^{(i)}
\label{eq:clasLinRegMat}
\end{aligned}</p>
<p>donde ${X}$ es una matriz de dimensiones $N \times p$, es decir con $N$ filas conteniendo cada una los $p$ predictores de cada observación, incluyendo una columna de $1$s para la intercepción, de manera que $X^{(i)}$ son los predictores de la observación en la posición $i$. A su vez $\beta$ es un vector con los coeficientes de la regresión, entonces ${X}^{(i)} {\beta}$ es un vector donde cada fila es el producto interno de predictores por coeficientes. Finalmente, el prerrequisito estándar $\eqref{eq:linRegEps}$ es que los errores $\epsilon^{(i)}$ sean independientes y que estén idénticamente distribuidos con distribución Normal y varianza $\sigma_{\epsilon}^2$.</p>
<p>Las regresiones multinivel (RMN) consideran la posibilidad de que la data puede separarse en grupos donde la regresión entre $y$ y sus predictores varía con el grupo. Por ejemplo, considere el rendimiento de los empleados de una compañía que varía entre oficinas y regiones. Otro ejemplo típico es la nota de alumnos agrupada por salón, curso, o colegio; más adelante estudiaremos un conjunto de datos reales de este tema.</p>
<p>De la misma forma como las regresiones clásicas son extendidas con el modelo linear generalizado (GLM) o aun más con modelo aditivo generalizado (GAM) las regresiones multinivel también son extendidas en esos contextos.</p>
<p>Lo interesante de las RMN es que este agrupamiento no es un prerrequisito. Aunque el diagnóstico de la regresión nos va a decir si el agrupamiento incluido en el modelo es <em>necesario</em> o no, el hecho de que no lo sea, no perjudica a la regresión resultante. Podemos decir que sin agrupamiento significativo la RMN revierte automáticamente a una sin grupos.</p>
<p><em>Conjunto de Datos "Toy"</em><span id="toyModel"></span></p>
<p>Para entender la mecánica de las RMN, practicaremos inicialmente con un conjunto de datos simulados denominado <em>"Toy"</em>, posteriormente estudiaremos datos reales. Los datos simulados constan de 504 observaciones de las variables <code>(x, y, g)</code>, donde -como ya supondrá- usaremos <code>x</code> como predictora; <code>y</code> como dependiente, y <code>g</code> es una variable categórica con el número del grupo al que pertenece cada observación (<a href="http://www.odelama.com/sandbox/xyg.csv">puede descargar la data aquí</a>, puede ver <a href="http://www.odelama.com/sandbox/#genToyData">el código con el fue generada aquí</a>, y el <a href="http://www.odelama.com/sandbox/#toyModel">modelo estadístico detrás aquí</a>). Tenemos 49 grupos, así el valor de <code>g</code> varía entre 1 y 49. El número de observaciones por grupo varía entre 1 y 21 ($n_j$), y el promedio es de 10 observaciones por grupo.</p>
<table class="standard center">
<thead>
<tr>
<th>x </th>
<th>y </th>
<th>g</th>
</tr>
</thead>
<tbody>
<tr>
<td>0.273 </td>
<td>8.9 </td>
<td>1</td>
</tr>
<tr>
<td>3.308 </td>
<td>15.6 </td>
<td>1</td>
</tr>
<tr>
<td>2.720 </td>
<td>15.2 </td>
<td>1</td>
</tr>
<tr>
<td>0.617 </td>
<td>11.6 </td>
<td>1</td>
</tr>
<tr>
<td>0.252 </td>
<td>8.0 </td>
<td>1</td>
</tr>
<tr>
<td>5.337 </td>
<td>21.3 </td>
<td>1</td>
</tr>
<tr>
<td>1.735 </td>
<td>12.4 </td>
<td>2</td>
</tr>
<tr>
<td>5.413 </td>
<td>27.1 </td>
<td>2</td>
</tr>
<tr>
<td>0.147 </td>
<td>8.8 </td>
<td>2</td>
</tr>
</tbody>
</table>
<h1 id="toyvariandosolointercepcinengrupos">Toy | Variando solo Intercepción en Grupos</h1>
<p>Para introducir la notación de las RMN, digamos que la data puede separarse en $J$ grupos. Inicialmente estudiaremos el caso en el que solo la intercepción de la regresión varía con el grupo:</p>
<p>\begin{align}
y^{(i)} = \alpha_0^{(j)} + {X}^{(i)} {\beta} +\epsilon^{(i)} ,~para~i=1, \dots ,N \
\label{eq:rmnInterc}
\end{align}</p>
<p>Como en el lado izquierdo de la ecuación \eqref{eq:rmnInterc} tenemos $y^{(i)}$, está implícito que el número de grupo $j$ en $\alpha_0^{(j)}$ se obtiene a partir de $i$, que es el dato de entrada. Es decir, en la matriz de datos hay una columna con el número de grupo de cada observación. En nuestro ejemplo, ésa es la variable <code>g</code>, que contiene número de grupo, por lo que $j=g[i]$. Para denotar esto podríamos usar $\alpha_0^{(g[i])}$ o tal vez $\alpha_0^{(j_i)}$, pero se hace más difícil la lectura de la ecuación y simplemente usamos $\alpha_0^{(j)}$ con la salvedad que acabamos de explicar.</p>
<p>Para ajustar el modelo a los datos podríamos usar el paquete R <code>lme4</code>, que emplea estimaciones puntuales de máxima verosimilitud, y para cálculos con inferencia bayesiana completa podemos usar <code>brms</code> o <code>rstanarm</code> o incluso <code>rstan</code>. La ventaja de <code>lme4</code> es que es muy rápida y se puede usarse para explorar los datos, pero para el modelo final es recomendable <code>brms</code> o <code>rstanarm</code>, quienes <em>por detrás</em> <a href="http://www.odelama.com/sandbox/#stan">llaman a <code>stan</code></a>. En nuestro caso <a href="http://www.odelama.com/sandbox/#Burkner-1">vamos a usar <code>brms</code></a>. Afortunadamente <code>lme4</code>, <code>brms</code> y <code>rstanarm</code> emplean la misma notación que emplearemos acá.</p>
<p>Para nuestro conjunto de datos, el modelo en la notación de \eqref{eq:rmnInterc} simplemente es</p>
<p>\begin{align}
y^{(i)} = \alpha_0^{(j)} + \beta_0 + {\beta_1}{x}^{(i)} +\epsilon^{(i)} ,~para~i=1, \dots ,N
\label{eq:rmnIntercEx1}
\end{align}</p>
<p>mientras que, el mismo, en la notación de fórmula de <code>brms</code> es</p>
<div class="highlight-r"><pre class="hl">y <span class="hl opt">~</span> <span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">|</span> g<span class="hl opt">)</span>
</pre></div>
<p>donde <code>1 + x</code> indica que queremos una regresión con intercepción <code>1</code> y con <code>x</code> como predictor, comunes a toda la población. Mientras que <code>(1 | g)</code> denota que queremos la intercepción <code>1</code> variando por cada grupo en <code>g</code>. El <code>1</code> global puede ser omitido, pues por defecto estas fórmulas (en R) asumen que se desea la intercepción, aunque es mejor ser explícito. Mas bien, cuando no se desea intercepción se debe usar <code>-1</code> o <code>0</code> como en <code>y ~ -1 + x + (1 | g)</code> o <code>y ~ 0 + x + (1 | g)</code>.</p>
<p>Antes de calcular la regresión de nuestro ejemplo debemos terminar de presentar las ecuaciones que explican cómo es calculado el modelo multinivel. Para ello vamos a darle más generalidad a la ecuación \eqref{eq:rmnInterc}: en vez de $\alpha_0^{(j)}$ vamos a usar $\pmb{\alpha}^{(j)}$ y luego agregar que $\alpha_0^{(j)}$ es su único componente (en este ejemplo):</p>
<p>\begin{align}
y^{(i)} &= \pmb{\alpha}^{(j)} + {X}^{(i)} {\beta} +\epsilon^{(i)} ,~para~i=1, \dots ,N \label{eq:rmnIntercF1} \\
\pmb{\alpha}^{(j)} &= \alpha_0^{(j)} + \eta^{(j)}, ~para~j=1, \dots ,J \label{eq:rmnIntercF2} \\
\eta^{(j)} &\sim \mathcal{N}(0, \sigma_\alpha^2) \label{eq:rmnIntercF3} \\
\epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2) \label{eq:rmnIntercF4}
\end{align}</p>
<p>El término $\eta^{(j)}$ en las ecuaciones \eqref{eq:rmnIntercF2} y \eqref{eq:rmnIntercF3} nos dice que todas las $J$ intercepciones $\alpha_0^{(j)}$ vienen de una misma distribución (normal) con varianza $\sigma_\alpha^2$. Otra forma de expresar estas dos ecuaciones, condensadas en una sola, es:</p>
<p>\begin{align}
\pmb{\alpha}^{(j)} \sim \mathcal{N}(\mu_{\alpha}, \sigma_\alpha^2) \label{eq:muRmnIntc}
\end{align}</p>
<p>Note que $\mu_{\alpha}$ y $\sigma_\alpha$ no dependen de $j$, lo cual indica que el origen de todas las intercepciones en los grupos es una misma distribución. Esto traerá beneficios que más adelante comentaremos.</p>
<p>Por otro lado, $\beta$ es el vector con los coeficientes que no varían con los grupos. A menudo se denomina <em>"efectos aleatorios"</em> a la intercepción y coeficientes que varían por grupo. En nuestra notación estos están en $\pmb{\alpha}^{(j)}$, aunque por el momento solo tenemos la intercepción. En esa misma terminología, a la intercepción y variables en ${X}^{(i)}$ se les denomina <em>"efectos fijos"</em>, y cuando un modelo tiene ambos efectos dicen que es de <em>"efectos mixtos"</em>. Sin embargo, nos referiremos a ellas como a nivel del grupo o de la población (globales), respectivamente.</p>
<p><em>Con Intercepción a Nivel de Grupo y Población.</em> </p>
<p>Note que el modelo \eqref{eq:rmnIntercEx1} incluye intercepción tanto a nivel población como por grupo ($\beta_0$ y $\alpha_0^{(j)}$). Esto suena a algo redundante, y de hecho lo es. Para resolver esto, el modelo multinivel va a estimar como intercepción a nivel de población a la media de las intercepciones $(\beta_0 = \mu_\alpha)$ y como intercepción a nivel de grupo $\alpha_0^{(j)}$ a la variación por grupo con respecto a esta media. Es decir, en este modelo $\mu_{\alpha}=0$ en la ecuación \eqref{eq:muRmnIntc} .</p>
<h2 id="ajustedelmodeloaladataspanidtoyfitintcspan">Ajuste del Modelo a la Data<span id="toyFitIntc"></span></h2>
<p>Ajustaremos el modelo con la función <code>brm()</code> en <a href="http://www.odelama.com/sandbox/#Burkner-1">el paquete <code>brms</code> (del Inglés <strong>b</strong>ayesian <strong>r</strong>egression <strong>m</strong>odels using <strong>s</strong>tan)</a>. Vamos a usar las distribuciones <em>anteriores</em> (<em>prior distributions</em>) predeterminadas de <code>brms</code> (más adelante hablaremos de estas distribuciones). Luego de ajustar el modelo imprimimos el resumen del ajuste.</p>
<div class="highlight-r"><pre class="hl">brm2_intc <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>y <span class="hl opt">~</span> <span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">|</span> g<span class="hl opt">),</span> data<span class="hl opt">,</span> iter<span class="hl opt">=</span><span class="hl num">2500</span><span class="hl opt">)</span>
<span class="hl kwd">summary</span><span class="hl opt">(</span>brm2_intc<span class="hl opt">)</span>
</pre></div>
<p><span id="toyIntcOutput"></span></p>
<pre><code>Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ 1 + x + (1 | g)
Data: data (Number of observations: 504)
Samples: 4 chains, each with iter = 2500; warmup = 1250; thin = 1;
total post-warmup samples = 5000
Group-Level Effects:
~g (Number of levels: 49)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 2.96 0.34 2.36 3.69 1.01 535 987
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 9.36 0.47 8.45 10.30 1.03 410 810
x 3.70 0.06 3.58 3.81 1.00 4237 3256
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 2.07 0.07 1.94 2.22 1.00 4317 3699
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
</code></pre>
<p>Note que en el resumen hay tres grandes bloques: </p>
<ol>
<li><code>Group-Level Effects</code>: Coeficientes a nivel de grupo, donde nuestro único grupo por la variable <code>g</code> ocupa la sección <code>~g</code>. </li>
<li><code>Population-Level Effects</code>: Coeficientes a nivel de población (globales). </li>
<li><code>Family Specific Parameters</code>; Parámetros auxiliares del modelo. </li>
</ol>
<p>Con los datos en esas secciones podemos ver que la realización de las ecuaciones \eqref{eq:rmnIntercF1}-\eqref{eq:rmnIntercF4}, para nuestro modelo es:</p>
<p>\begin{align*}
y^{(i)} &= \pmb{\alpha}^{(j)} + 9.36 + 3.70 x +\epsilon^{(i)} ,~para~i=1, \dots ,504 \\
\pmb{\alpha}^{(j)} &= \alpha_0^{(j)} + \eta^{(j)}, ~para~j=1, \dots ,49 \\
\eta^{(j)} &\sim \mathcal{N}(0, 2.96^2) \\
\epsilon^{(i)} &\sim \mathcal{N}(0, 2.07^2)
\end{align*}</p>
<p>Antes de ver más detalles de qué hemos obtenido al ajustar el modelo con inferencia Bayesiana completa, necesitamos entender mejor en qué consiste ésta y cómo se analizan sus resultados.</p>
<h1 id="inferenciabayesianacompleta">Inferencia Bayesiana Completa</h1>
<p>Como <code>brms</code> ajusta modelos usando inferencia Bayesiana completa, varios de los pasos siguientes en este artículo refieren a términos y análisis propios de ese enfoque. Vamos a hacer un alto al modelado y en esta sección vamos a revisar este tema, que puede saltar si ya está familiarizado con esto. ¡Si no es así, bienvenido! Está en un buen lugar para aprender esto desde cero. Al principio el tema es un poco denso. Por favor tenga paciencia y en breve verá como todas las piezas encajan.</p>
<p><em>Modelo de Datos eCompo</em><span id="eCompoModel"></span></p>
<p>En esta sección usaremos un <a href="http://www.odelama.com/sandbox/simple.rda">conjunto de datos denominado <em>"eCompo"</em></a>. Se trata de una serie de lecturas "$x$" de un componente electrónico, que en una prueba ha sido sometido a diferentes valores de un estímulo "$y$". Queremos conocer los dos parámetros $\beta$ y $\sigma_\epsilon$ en el modelo que se muestra a continuación:</p>
<p>\begin{align}
y &= \beta \cdot x + \epsilon \\
\epsilon &\sim \mathcal{N}(0, \sigma_{\epsilon}^2) \label{eq:demoPostLik} \\
\implies y - (\beta \cdot x) &\sim \mathcal{N}(0, \sigma_{\epsilon}^2) \label{eq:demoPostEps}
\end{align}</p>
<p>La distribución esperada de los residuos \eqref{eq:demoPostLik} tiene la varianza de éstos, pero está centrada en cero. Es decir, se espera que <em>el centro</em> de los residuos sea cero. Vamos a estimar la varianza $\sigma_{\epsilon}^2$ con la de los residuos en la data. </p>
<p><em>Ejercicio para la Intuir la Probabilidad Conjunta</em></p>
<p>En el gráfico <em>"Beta vs Residuos"</em> (abajo), se muestra con color verde la distribución esperada de los residuos, y con morado su distribución real $\epsilon = y - (\beta \cdot x)$ según la data. Con $\beta=9.3$ en (1) los residuos se aproximan a la distribución esperada. Conforme $\beta$ se incrementa, como hasta $\beta=15.3$ en (2), los residuos se esparcen y su centro se mueve hacia la izquierda, alejándose de lo esperado. Como su varianza aumenta, su distribución esperada se aplana <em>tratando de seguirlos</em>, pero el calce es cada vez peor: la distribución real es cada vez más diferente a la esperada. En este último caso decimos que $\beta=15$ "hace poco probable <em>observar</em> los datos que tenemos".</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/moving-beta.png">
<img alt="Beta vs los Residuos" class="pic-border" data-width="750" src="http://www.odelama.com/sandbox/moving-beta.png" title="Beta vs los Residuos" />
</a>
<p class="caption-text"><em>Beta vs Residuos</em>. En verde se muestra la distribución esperada de los residuos según \eqref{eq:demoPostLik}: con centro en cero y la varianza de éstos. En morado su distribución real. En el momento (1) el valor de $\beta$ produce residuos cuya distribución es relativamente acorde a la esperada. En (2) valores crecientes de $\beta$ causan residuos cada vez más negativos con una varianza creciente. Su distribución esperada se aplana <em>tratando de seguir a los residuos</em>, pero el calce es cada vez menor.</p>
</div>
<p>La idea de este ejercicio es motivar la intuición de que los parámetros de un modelo tienen <em>la distribución conjunta</em> de la probabilidad de observar la data según el modelo. Cuando movemos uno o más parámetros, usamos la data que tenemos para ver como varía <em>"la probabilidad de verla"</em> según el modelo. Esto es, analizamos la probabilidad conjunta de los parámetros <em>"condicionados a la data"</em>. </p>
<p>Aun cuando en el sumario de <code>brm()</code> hemos visto estimaciones puntuales de los parámetros de nuestro <a href="http://www.odelama.com/sandbox/#eCompoModel">modelo eCompo</a>, en realidad los ha calculado a partir de muestras de la distribución conjunta de esos parámetros <em>condicionados a la data</em>. Veamos a continuación como obtener esa distribución.</p>
<p>La expresión formal de la probabilidad conjunta de los parámetros de nuestro modelo multinivel ( $\alpha_0^{(j)}, \beta, \sigma_\alpha, \sigma_{\epsilon}$) condicionada a la data ($y, X$) es:</p>
<p>$$\begin{aligned}
\pi(\alpha_0^{(j)}, \beta, \sigma_\alpha, \sigma_{\epsilon} \mid y, X)
\end{aligned}$$</p>
<p><em>Notación Simplificada</em> <span id="notaSimpli"></span></p>
<p>Por simplicidad, llamaremos $\theta$ al vector con todos los parámetros del modelo. A su vez, con $y$ representaremos a toda la data, así que $X$ estará también implícita. Finalmente, en estas expresiones, $\pi(\cdot)$ es la <em>función de densidad</em> de $(\cdot)$, es decir es una función diferente para cada $(\cdot)$ diferente, que corresponde a una <em>función de masa</em> en el caso discreto. Usaremos el término <em>distribución</em> con el mismo sentido que <em>función de densidad</em>.</p>
<p>Con esta notación simplificada, lo que buscamos es:</p>
<p>$$\begin{aligned}
\pi(\theta \mid y)
\end{aligned}$$</p>
<p>Por la regla de Bayes sabemos que:</p>
<p>\begin{align}
\pi(\theta \mid y) &= \frac{\pi(y, \theta)}{\pi(y)} = \frac{\pi(y \mid \theta) \cdot \pi(\theta)}{\pi(y)} \label{eq:posterior} \\
\pi(\theta \mid y) &\propto \pi(y \mid \theta) \cdot \pi(\theta) \label{eq:propPost}
\end{align}</p>
<p>En \eqref{eq:posterior}, el numerador es la probabilidad conjunta de los parámetros y la data, a la cual nos hemos referido en el ejercicio anterior. El denominador $\pi(y)$ es una constante que no depende de $\theta$ sino solo de la data, por lo que $\pi(\theta \mid y)$ en \eqref{eq:propPost} es proporcional al numerador de \eqref{eq:posterior}. Podemos considerar el lado derecho \eqref{eq:propPost} como la distribución no normalizada de $\pi(\theta \mid y)$, es decir cuya integral no es igual a $1$.</p>
<p>La ecuación \eqref{eq:propPost} es clave en la inferencia Bayesiana, y cada término tiene un nombre al que nos referiremos frecuentemente. </p>
<ul>
<li>$\pi(\theta)$ es la <em>distribución "<strong>anterior</strong>"</em> (<em>prior distribution</em>): la que asumimos tienen los parámetros antes de ver la data. Nos referiremos a esta distribución como <em>anterior</em> o <em>prior</em>.</li>
<li>$\pi(y \mid \theta)$ es la distribución de muestreo de la data <em>(sampling distribution)</em>. Aquí miramos este término como función de $\theta$, pues la data $y$ se mantiene constante, y en ese sentido este término es llamado la <em>"<strong>verosimilitud</strong>" (likelihood)</em> de nuestro modelo. </li>
<li>$\pi(\theta \mid y)$ es la <em>distribución "<strong>posterior</strong>"</em>: Es la distribución de los parámetros condicionados a la data. A la que nos referimos como la probabilidad de ver la data que tenemos, con el modelo y parámetros que tenemos.</li>
</ul>
<p>Como $\theta$ es un vector, estos términos corresponden a distribuciones <em>conjuntas</em> de probabilidad de todos ellos. Por no haber aun introducido la distribución <em>"anterior"</em>, implícitamente nos hemos referido a la <em>verosimilitud</em> cuando hemos hablado de la <em>"posibilidad de observar la data"</em> en el ejercicio anterior, que corresponde al caso cuando la distribución <em>anterior</em> es uniforme.</p>
<h2 id="ejemplodeinferenciabayesianasinstan">Ejemplo de Inferencia Bayesiana sin Stan</h2>
<p>En esta sección veremos un ejemplo práctico de los conceptos introducidos en la sección anterior, analizando numéricamente (sin usar directa ni indirectamente <code>stan</code>) <a href="http://www.odelama.com/sandbox/#eCompoModel">el modelo eCompo</a>.</p>
<p><em>Distribución Anterior</em></p>
<p>Los expertos que hemos consultado nos dicen que $\beta$ debe de estar —a lo mucho— entre $1$ y $25$, con una distribución aproximadamente Normal, lo cual implica su centro en $(25+1)/2=13$, con una desviación estándar (al 95%) de $(25-1)/4=6$, por lo que asumimos como <em>anterior</em> de $\beta$ una distribución Normal con esos parámetros. Los expertos también nos dicen, que la distribución Exponencial es adecuada como <em>prior</em> de una desviación estándar usando como parámetro $\lambda$ de la misma al recíproco de la desviación estándar de $y$, por lo que adoptamos como <em>anterior</em> de $\sigma_\epsilon$ a tal distribución con $\lambda=1/sd(y)$.</p>
<p>El <em>prior</em> $\pi(\theta)$ de este modelo es:</p>
<div class="highlight-r"><pre class="hl">prior_beta_mean <span class="hl opt"><-</span> <span class="hl num">13</span>
prior_beta_sd <span class="hl opt"><-</span> <span class="hl num">6</span>
prior_sigma_lambda <span class="hl opt"><-</span> <span class="hl num">1</span> <span class="hl opt">/</span> <span class="hl kwd">sd</span><span class="hl opt">(</span>y<span class="hl opt">)</span>
prior_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) {</span>
log_prior_beta <span class="hl opt"><-</span> <span class="hl kwd">dnorm</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> mean <span class="hl opt">=</span> prior_beta_mean<span class="hl opt">,</span> sd <span class="hl opt">=</span> prior_beta_sd<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span>
log_prior_sigma <span class="hl opt"><-</span> <span class="hl kwd">dexp</span><span class="hl opt">(</span>sigma<span class="hl opt">,</span> rate <span class="hl opt">=</span> prior_sigma_lambda<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span>
result <span class="hl opt"><-</span> log_prior_beta <span class="hl opt">+</span> log_prior_sigma
<span class="hl kwa">if</span> <span class="hl opt">(</span>log <span class="hl opt">!=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span> result <span class="hl opt"><-</span> <span class="hl kwd">exp</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl kwd">return</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl opt">}</span>
</pre></div>
<p>Donde <code>beta</code> y <code>sigma</code> en el código corresponden a $\beta$ y $\sigma_\epsilon$ en nuestro modelo. Como estos parámetros son independientes, su <em>prior</em> conjunto es el producto de los <em>priors</em> de cada parámetro. Tal como vimos, esta distribución no depende de los datos, sino solo de $\theta$, y de los <em>hiper-parametros</em> con los que modelamos los <em>priors</em> individuales (como <code>prior_beta_sd</code> o <code>prior_sigma_lambda</code>), que forman parte de nuestro modelo. Como un tema meramente técnico (estabilidad y precisión numérica), los cálculos son hechos en el espacio logarítmico.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/prior_density.png">
<img alt="Distribución Anterior" class="pic-border" data-width="850" src="http://www.odelama.com/sandbox/prior_density.png" title="Distribución Anterior" />
</a>
<p class="caption-text"><em>Distribución Anterior</em>. La distribución <em>anterior</em> cubre el área de valores plausibles de los parámetros del modelo.</p>
</div>
<p>Una distribución <em>anterior</em> debe cubrir con holgura el dominio de los posibles valores de los parámetros, que en el el gráfico <em>"Distribución Anterior"</em> (arriba) vemos como ocurre en nuestro ejemplo. </p>
<p><em>Distribución de Verosimilitud</em></p>
<p>La distribución de la <em>verosimilitud</em> $\pi(y \mid \theta)$ es:</p>
<div class="highlight-r"><pre class="hl">likelihood_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> y<span class="hl opt">,</span> x<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) {</span>
result <span class="hl opt"><-</span> <span class="hl num">0</span>
<span class="hl kwa">for</span> <span class="hl opt">(</span>i <span class="hl kwa">in</span> <span class="hl kwd">seq_along</span><span class="hl opt">(</span>y<span class="hl opt">)) {</span>
mean <span class="hl opt"><-</span> beta <span class="hl opt">*</span> x<span class="hl opt">[</span>i<span class="hl opt">]</span>
result <span class="hl opt"><-</span> result <span class="hl opt">+</span> <span class="hl kwd">dnorm</span><span class="hl opt">(</span>y<span class="hl opt">[</span>i<span class="hl opt">],</span> mean <span class="hl opt">=</span> mean<span class="hl opt">,</span> sd <span class="hl opt">=</span> sigma<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span>
<span class="hl opt">}</span>
<span class="hl kwa">if</span> <span class="hl opt">(</span>log <span class="hl opt">!=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span> result <span class="hl opt"><-</span> <span class="hl kwd">exp</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl kwd">return</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl opt">}</span>
</pre></div>
<p>Como esperábamos, la <em>verosimilitud</em> es función tanto de la data como de $\theta$. En nuestro modelo, la ecuación \eqref{eq:demoPostLik} define la probabilidad de ver cada observación en la data. Como éstas son independientes su probabilidad conjunta es el producto de la <em>verosimilitud</em> de todas las observaciones. </p>
<p>Note que, estrictamente hablando, <code>likelihood_density()</code> y <code>prior_density()</code> son funciones de densidad no normalizadas. Para que sean normalizadas se requiere que la integral, con respecto $\theta$, de los valores que produce cada una de ellas sea $1$. Por ejemplo, en el caso de la distribución <em>anterior</em>, esto es $\int_{-\infty}^{+\infty} \pi(\theta)~\textrm{d}\theta = 1$, lo cual no es cierto. Solo falta dividir el resultado de cada función por el valor de su integral. Más adelante nos ocuparemos de esto.</p>
<p><em>Distribución Conjunta</em></p>
<p>Con lo calculado anteriormente, la distribución conjunta $\pi(y \mid \theta) \cdot \pi(\theta)$ es trivial:</p>
<div class="highlight-r"><pre class="hl">joint_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) {</span>
result <span class="hl opt"><-</span>
<span class="hl kwd">prior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) +</span>
<span class="hl kwd">likelihood_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data$y<span class="hl opt">,</span> data$x<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span>
<span class="hl kwa">if</span> <span class="hl opt">(</span>log <span class="hl opt">!=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span> result <span class="hl opt"><-</span> <span class="hl kwd">exp</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl kwd">return</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl opt">}</span>
</pre></div>
<p><em>Distribución Posterior</em></p>
<p>La distribución conjunta $\pi(y \mid \theta) \cdot \pi(\theta)$ calculada por la función <code>joint_density()</code> es todo lo que necesitamos en \eqref{eq:propPost}. Aun así vamos a calcular $\pi(y)$, mostrando que no es más que la integral de <code>joint_density()</code> sobre $\theta$, una constante con la que la dividiremos para normalizarla y obtener así la distribución de densidad <em>posterior</em> $\pi(\theta \mid y)$ según \eqref{eq:posterior}. </p>
<p>$$\begin{aligned}
\pi(y) &= \int_{-\infty}^{+\infty} \pi(y \mid \theta) \cdot \pi(\theta)~\textrm{d}\theta \\
\pi(y) &= \int_0^{+\infty} \int_{-\infty}^{+\infty} \textrm{joint_density}(\beta, \sigma_\epsilon \mid data)~\textrm{d}\beta~\textrm{d}\sigma_\epsilon
\end{aligned}$$</p>
<p>La distribución <em>posterior</em> $\pi(\theta \mid y)$ tiene la forma:</p>
<div class="highlight-r"><pre class="hl">posterior_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) {</span>
result <span class="hl opt"><-</span> <span class="hl kwd">joint_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">) -</span> <span class="hl kwd">log</span><span class="hl opt">(</span>pi_y<span class="hl opt">)</span>
<span class="hl kwa">if</span> <span class="hl opt">(</span>log <span class="hl opt">!=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span> result <span class="hl opt"><-</span> <span class="hl kwd">exp</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl kwd">return</span><span class="hl opt">(</span>result<span class="hl opt">)</span>
<span class="hl opt">}</span>
</pre></div>
<p>Donde <code>pi_y</code> es el valor de la integral que debemos calcular, es decir corresponde a $\pi(y)$. Tenemos la dificultad de que cuando integramos <code>pi_y = 1</code> la función <code>posterior_density()</code> retorna valores menores que $3.4 \times 10^{-47}$ y la integración es numéricamente inestable, por lo que vamos a asignarle a <code>pi_y</code> un valor temporal mientras integramos y luego ajustaremos su valor con el resultado de la integral.</p>
<div class="highlight-r"><pre class="hl">pi_y <span class="hl opt"><-</span> <span class="hl num">3.4e-47</span>
int_pi_y <span class="hl opt"><-</span> <span class="hl kwd">integrate</span><span class="hl opt">(</span><span class="hl kwa">function</span><span class="hl opt">(</span>sigmas<span class="hl opt">) {</span>
<span class="hl kwd">sapply</span><span class="hl opt">(</span>sigmas<span class="hl opt">,</span> <span class="hl kwa">function</span><span class="hl opt">(</span>sigma<span class="hl opt">) {</span>
<span class="hl kwd">integrate</span><span class="hl opt">(</span>
<span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">)</span> <span class="hl kwd">posterior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">F</span><span class="hl opt">),</span>
<span class="hl opt">-</span><span class="hl kwb">Inf</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span>
<span class="hl opt">)</span>$value
<span class="hl opt">})</span>
<span class="hl opt">},</span> <span class="hl num">0</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span><span class="hl opt">)</span>$value
pi_y <span class="hl opt"><-</span> pi_y <span class="hl opt">*</span> int_pi_y
</pre></div>
<p>El valor final es $\pi(y)=2.798 \times 10^{-46}$, con lo que <code>posterior_density()</code> es ahora una función de densidad normalizada. Ahora podemos regresar al ejercicio inicial y calcular la probabilidad de (1) $\pi(\beta=9.3, \sigma_\epsilon=36.6 \mid y)=0.071$ y de (2) $\pi(\beta=15, \sigma_\epsilon=55.1 \mid y)=1.6 \times 10^{-25}$. Donde vemos que es practicamente imposible <em>observar</em> la data que tenemos con $\beta=15$.</p>
<p><em>Gráfico de la distribución Posterior</em></p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/posterior_density.png">
<img alt="Distribución Posterior" class="pic-border" data-width="900" src="http://www.odelama.com/sandbox/posterior_density.png" title="Distribución Posterior" />
</a>
<p class="caption-text"><em>Distribución Posterior</em>. Distribución <em>posterior</em> calculada numéricamente, con algunas líneas de contorno para notar mejor su forma. Puede verse que $\sigma_\epsilon$ tiene una cola más larga a la derecha.</p>
</div>
<p>En el gráfico <em>"Proporcional al Posterior"</em> (arriba), vemos que en nuestro caso $\pi(\theta)$ tiene un punto modal prominente. Algunos algoritmos de estimación puntual de tipo <em>máxima verosimilitud</em> estiman la posición de esta moda como solución a los parámetros del modelo. En una estimación Bayesiana eso corresponde a la moda de la <em>posterior</em> cuando los <em>priors</em> son uniformes, pero nosotros no usamos solo las estimaciones puntuales de $\theta$, sino toda su distribución, incorporando en nuestros cálculos la incertidumbre que el modelo tiene acerca de su valor. </p>
<p>Es posible obtener la posición de la moda de la <em>posterior</em> usando la función <code>optimizing()</code> de <code>rstan</code>. En nuestro caso podemos hallarla con un optimizador numérico en R.</p>
<div class="highlight-r"><pre class="hl">res <span class="hl opt"><-</span> <span class="hl kwd">nlm</span><span class="hl opt">(</span><span class="hl kwa">function</span><span class="hl opt">(</span>theta<span class="hl opt">) -</span><span class="hl kwd">posterior_density</span><span class="hl opt">(</span>theta<span class="hl opt">[</span><span class="hl num">1</span><span class="hl opt">],</span> theta<span class="hl opt">[</span><span class="hl num">2</span><span class="hl opt">],</span> data<span class="hl opt">),</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">7</span><span class="hl opt">,</span> <span class="hl num">11</span><span class="hl opt">))</span>
res
</pre></div>
<pre><code>$minimum
[1] 2.61698
$estimate
[1] 9.303106 35.262023
$gradient
[1] -1.038726e-07 -6.448118e-09
</code></pre>
<p>Como <code>nlm</code> (Non-Linear Minimization) calcula mínimos, invertimos el signo de <code>posterior_density</code>. Asimismo, como <code>log</code> es una función monotónica optimizamos en el espacio logarítmico por conveniencia numérica. Obtenemos $\beta=9.3031$ y $\sigma_\epsilon=35.2620$.</p>
<h2 id="comparacinconajusteobtenidoconstan">Comparación con Ajuste Obtenido con Stan</h2>
<p>En este artículo, obtendremos muestras de <em>la posterior</em> con el apoyo de <code>stan</code>, mediante <code>brms</code>. A continuación, vamos a ajustar <a href="http://www.odelama.com/sandbox/#eCompoModel">el modelo eComp</a> directamente con <code>stan</code> (y no vía <code>brms</code>) para mejor control de todos los detalles.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">library</span><span class="hl opt">(</span>rstan<span class="hl opt">)</span>
<span class="hl kwd">options</span><span class="hl opt">(</span>mc<span class="hl opt">.</span>cores <span class="hl opt">=</span> parallel<span class="hl opt">::</span><span class="hl kwd">detectCores</span><span class="hl opt">())</span>
rstan<span class="hl opt">::</span><span class="hl kwd">rstan_options</span><span class="hl opt">(</span>auto_write <span class="hl opt">=</span> <span class="hl kwb">TRUE</span><span class="hl opt">)</span>
model<span class="hl opt">.</span>stan <span class="hl opt"><-</span> <span class="hl str">"</span>
<span class="hl str">data {</span>
<span class="hl str"> int<lower=2> N;</span>
<span class="hl str"> vector[N] x;</span>
<span class="hl str"> vector[N] y;</span>
<span class="hl str"> real prior_beta_sd;</span>
<span class="hl str"> real prior_beta_mean;</span>
<span class="hl str"> real<lower=0> prior_sigma_lambda;</span>
<span class="hl str">}</span>
<span class="hl str"></span>
<span class="hl str">parameters {</span>
<span class="hl str"> real beta;</span>
<span class="hl str"> real<lower=0> sigma;</span>
<span class="hl str">}</span>
<span class="hl str"></span>
<span class="hl str">model {</span>
<span class="hl str"> vector[N] mu;</span>
<span class="hl str"> beta ~ normal(prior_beta_mean, prior_beta_sd);</span>
<span class="hl str"> sigma ~ exponential(prior_sigma_lambda);</span>
<span class="hl str"> mu = beta * x;</span>
<span class="hl str"> y ~ normal(mu, sigma);</span>
<span class="hl str">}</span>
<span class="hl str"></span>
<span class="hl str">"</span>
data<span class="hl opt">.</span>stan <span class="hl opt"><-</span> <span class="hl kwd">list</span><span class="hl opt">(</span>
N <span class="hl opt">=</span> <span class="hl kwd">length</span><span class="hl opt">(</span>y<span class="hl opt">),</span>
x <span class="hl opt">=</span> x<span class="hl opt">,</span>
y <span class="hl opt">=</span> y<span class="hl opt">,</span>
prior_beta_sd <span class="hl opt">=</span> prior_beta_sd<span class="hl opt">,</span>
prior_beta_mean <span class="hl opt">=</span> prior_beta_mean<span class="hl opt">,</span>
prior_sigma_lambda <span class="hl opt">=</span> prior_sigma_lambda
<span class="hl opt">)</span>
fit<span class="hl opt">.</span>stan <span class="hl opt"><-</span> <span class="hl kwd">stan</span><span class="hl opt">(</span>model_code <span class="hl opt">=</span> model<span class="hl opt">.</span>stan<span class="hl opt">,</span> data <span class="hl opt">=</span> data<span class="hl opt">.</span>stan<span class="hl opt">,</span> chains<span class="hl opt">=</span><span class="hl num">4</span><span class="hl opt">,</span> iter<span class="hl opt">=</span><span class="hl num">2250</span><span class="hl opt">,</span> warmup<span class="hl opt">=</span><span class="hl num">1000</span> <span class="hl opt">,</span>refresh <span class="hl opt">=</span> <span class="hl num">0</span><span class="hl opt">)</span>
stan_samples <span class="hl opt"><-</span> as<span class="hl opt">.</span>data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>fit<span class="hl opt">.</span>stan<span class="hl opt">)</span>
<span class="hl kwd">summary</span><span class="hl opt">(</span>fit<span class="hl opt">.</span>stan<span class="hl opt">,</span> probs<span class="hl opt">=</span><span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">.25</span><span class="hl opt">,</span><span class="hl num">.50</span><span class="hl opt">,</span><span class="hl num">.75</span><span class="hl opt">))</span>$summary <span class="hl opt">%>%</span> <span class="hl kwd">print</span><span class="hl opt">(</span>digits<span class="hl opt">=</span><span class="hl num">4</span><span class="hl opt">)</span>
</pre></div>
<p><span id="stanEComp"></span></p>
<pre><code> mean se_mean sd 25% 50% 75% n_eff Rhat
beta 9.305 0.006899 0.3983 9.035 9.312 9.568 3332 1
sigma 38.647 0.121053 6.7265 33.899 37.826 42.289 3088 1
lp__ -104.999 0.025306 1.0639 -105.401 -104.684 -104.250 1768 1
</code></pre>
<p>En nuestro código arriba, la variable <code>stan_samples</code> contiene las $5000$ muestras recibidas de <code>stan</code>, que está compuesta de los valores de las variables <code>beta</code>, <code>sigma</code> y <code>lp__</code>, tal como puede verse en el listado de salida. En <code>stan_samples</code>, la variable <code>lp__</code> contiene el logaritmo de <em>la posterior</em> de nuestro modelo más los ajuste por <a href="https://mc-stan.org/docs/2_24/stan-users-guide/changes-of-variables.html">los Jacobianos requeridos para pasar parámetros de una escala restringida a una no restringida</a>. En nuestro caso $\sigma_\epsilon$ es la única variable restringida (a ser positiva), lo cual compensamos restando su logaritmo.</p>
<p>Si lo anterior no quedó muy claro, digamos que <code>lp__</code> no es exactamente el logaritmo de la posterior como la calculamos nosotros, sino que contiene además ciertas adiciones de <code>stan</code> que debemos deshacer antes de poder compararla con la que calculamos nosotros.</p>
<div class="highlight-r"><pre class="hl">lpd_num <span class="hl opt"><-</span> purrr<span class="hl opt">::</span><span class="hl kwd">map2_dbl</span><span class="hl opt">(</span>
stan_samples$beta<span class="hl opt">,</span> stan_samples$sigma<span class="hl opt">,</span>
<span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">)</span> <span class="hl kwd">posterior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">)</span>
<span class="hl opt">)</span>
plot_data <span class="hl opt"><-</span> data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>
lpd_num <span class="hl opt">=</span> lpd_num<span class="hl opt">,</span>
lpd_stan <span class="hl opt">=</span> stan_samples$lp__<span class="hl opt">,</span>
beta <span class="hl opt">=</span> stan_samples$beta<span class="hl opt">,</span>
sigma <span class="hl opt">=</span> stan_samples$sigma
<span class="hl opt">) %>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>
diff <span class="hl opt">=</span> lpd_num <span class="hl opt">- (</span>lpd_stan <span class="hl opt">-</span> <span class="hl kwd">log</span><span class="hl opt">(</span>sigma<span class="hl opt">))</span>
<span class="hl opt">)</span>
<span class="hl kwd">sd</span><span class="hl opt">(</span>plot_data$diff<span class="hl opt">)</span>
</pre></div>
<pre><code>[1] 3.235096e-14
</code></pre>
<p>En el código arriba, en <code>lpd_num</code> calculamos el logaritmo de <em>la posterior</em> para los pares de $\beta$ y $\sigma_\epsilon$ en las muestras devueltas por <code>stan</code>, a fin de comparar el resultado de nuestros cálculos con los de <code>stan</code> en <code>lp__</code>. Al calcular en <code>diff</code> la diferencia entre los logaritmos de <em>las posteriores</em>, a la de <code>stan</code> le restamos el logaritmo de <code>sigma</code> por la razón comentada antes. Luego de eso, ambos logaritmos deberían diferir solo por una constate, lo cual verificamos es lo que ocurre, pues la desviación estándar de tal diferencia es prácticamente cero.</p>
<p>En otras palabras, nuestros cálculos de <em>la posterior</em>, <em>verosimilitud</em> y <em>posterior</em> coinciden perfectamente con los que hace <code>stan</code> ✅.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/prop_to_posterior.png">
<img alt="Niveles de ICC" class="pic-border" data-width="900" src="http://www.odelama.com/sandbox/prop_to_posterior.png" title="Niveles de ICC" />
</a>
<p class="caption-text"><em>Proporcional al Posterior</em>. Muestras recibidas de <code>stan</code>. La distribución calculada por <code>stan</code> es igual a la que calculamos nosotros y mostramos en el gráfico anterior. La escala vertical aquí es diferente a la nuestra debido a que ésta es calculada con \eqref{eq:propPost} y la nuestra con \eqref{eq:posterior}.</p>
</div>
<p>La distribución hallada por <code>stan</code> es igual a la que calculamos nosotros y mostramos en el gráfico <em>"Distribución Posterior"</em> (más arriba). La escala vertical aquí es diferente a la nuestra, pues aquí corresponde a \eqref{eq:propPost} y la nuestra a \eqref{eq:posterior}. Aunque eso es irrelevante, pues lo que importa es la posición de las muestras en el espacio de $\theta$.</p>
<p><em>Distribución Marginal</em></p>
<p>Si bien los parámetros tienen una distribución conjunta, es decir el valor de cada uno condiciona el de los demás, podemos obtener su distribución marginal integrando la distribución conjunta <em>posterior</em> sobre los demás parámetros.</p>
<p>\begin{align*}
\pi(\beta) &= \int_0^{+\infty} \textrm{posterior_density}(\beta, \sigma_\epsilon \mid data) ~\textrm{d}\sigma_\epsilon \\
\pi(\sigma_\epsilon) &= \int_{-\infty}^{+\infty} \textrm{posterior_density}(\beta, \sigma_\epsilon \mid data)~\textrm{d}\beta
\end{align*}</p>
<p>Integramos numéricamente la distribución <em>posterior</em> para obtener las marginales, luego las graficamos, comparándolas con los valores obtenidos con <code>stan</code>:</p>
<div class="highlight-r"><pre class="hl">sigma_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>sigmas<span class="hl opt">) {</span>
<span class="hl kwd">sapply</span><span class="hl opt">(</span>sigmas<span class="hl opt">,</span> <span class="hl kwa">function</span><span class="hl opt">(</span>sigma<span class="hl opt">) {</span>
<span class="hl kwd">integrate</span><span class="hl opt">(</span><span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">)</span> <span class="hl kwd">posterior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">F</span><span class="hl opt">), -</span><span class="hl kwb">Inf</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span><span class="hl opt">)</span>$value
<span class="hl opt">})</span>
<span class="hl opt">}</span>
<span class="hl slc"># Asegurar el area bajo la curva integrada sea igual a 1</span>
sigmas <span class="hl opt"><-</span> <span class="hl kwd">seq</span><span class="hl opt">(</span><span class="hl num">22.2</span><span class="hl opt">,</span> <span class="hl num">75</span><span class="hl opt">,</span> length<span class="hl opt">.</span>out <span class="hl opt">=</span> <span class="hl num">200</span><span class="hl opt">)</span>
plot_sigma <span class="hl opt"><-</span> data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>sigma<span class="hl opt">=</span>sigmas<span class="hl opt">,</span> dens<span class="hl opt">=</span><span class="hl kwd">sigma_density</span><span class="hl opt">(</span>sigmas<span class="hl opt">))</span>
diffs <span class="hl opt"><-</span> <span class="hl kwd">diff</span><span class="hl opt">(</span>sigmas<span class="hl opt">)</span>
diffs <span class="hl opt"><-</span> <span class="hl kwd">c</span><span class="hl opt">(</span>diffs<span class="hl opt">[</span><span class="hl num">1</span><span class="hl opt">],</span> diffs<span class="hl opt">)</span>
area <span class="hl opt"><-</span> <span class="hl kwd">sum</span><span class="hl opt">(</span>diffs <span class="hl opt">*</span> plot_sigma$dens<span class="hl opt">)</span>
plot_sigma$dens <span class="hl opt"><-</span> plot_sigma$dens<span class="hl opt">/</span> area
plot_sigma <span class="hl opt">%>%</span>
<span class="hl kwd">ggplot</span><span class="hl opt">() +</span>
<span class="hl kwd">geom_line</span><span class="hl opt">(</span><span class="hl kwd">aes</span><span class="hl opt">(</span>sigma<span class="hl opt">,</span> dens<span class="hl opt">),</span> color <span class="hl opt">=</span> <span class="hl str">"darkblue"</span><span class="hl opt">) +</span>
<span class="hl kwd">geom_ribbon</span><span class="hl opt">(</span>
<span class="hl kwd">aes</span><span class="hl opt">(</span>sigma<span class="hl opt">,</span> ymax <span class="hl opt">=</span> dens<span class="hl opt">,</span> fill <span class="hl opt">=</span> <span class="hl str">"numeric"</span><span class="hl opt">),</span>
ymin <span class="hl opt">=</span> <span class="hl num">0</span><span class="hl opt">,</span> alpha <span class="hl opt">=</span> <span class="hl num">0.4</span>
<span class="hl opt">) +</span>
<span class="hl kwd">geom_histogram</span><span class="hl opt">(</span>
data <span class="hl opt">=</span> stan_samples<span class="hl opt">,</span>
<span class="hl kwd">aes</span><span class="hl opt">(</span>x<span class="hl opt">=</span> sigma<span class="hl opt">,</span> y <span class="hl opt">= ..</span>density<span class="hl opt">..,</span> fill <span class="hl opt">=</span> <span class="hl str">"stan"</span><span class="hl opt">),</span> alpha <span class="hl opt">=</span> <span class="hl num">0.4</span><span class="hl opt">,</span> bins<span class="hl opt">=</span><span class="hl num">45</span>
<span class="hl opt">) +</span>
<span class="hl kwd">scale_fill_manual</span><span class="hl opt">(</span>values <span class="hl opt">=</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl str">"numeric"</span> <span class="hl opt">=</span> <span class="hl str">"blue"</span><span class="hl opt">,</span> <span class="hl str">"stan"</span> <span class="hl opt">=</span> <span class="hl str">"darkgreen"</span><span class="hl opt">)) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>
title <span class="hl opt">=</span> <span class="hl str">"Distribución Marginal de Sigma"</span><span class="hl opt">,</span>
x <span class="hl opt">=</span> <span class="hl str">"Sigma"</span><span class="hl opt">,</span>
y <span class="hl opt">=</span> <span class="hl str">"Densidad"</span>
<span class="hl opt">) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>fill<span class="hl opt">=</span><span class="hl str">"Método"</span><span class="hl opt">) +</span>
<span class="hl kwd">theme</span><span class="hl opt">(</span>legend<span class="hl opt">.</span>position <span class="hl opt">=</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">0.7</span><span class="hl opt">,</span> <span class="hl num">0.6</span><span class="hl opt">))</span>
</pre></div>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/beta-sigma-marginal.png">
<img alt="Distribución Marginal" class="pic-border" data-width="900" src="http://www.odelama.com/sandbox/beta-sigma-marginal.png" title="Distribución Marginal" />
</a>
<p class="caption-text"><em>Distribuciones Marginales</em>. .</p>
</div>
<p>Estas curvas muestran que la densidad de las muestras que devuelve <code>stan</code> es proporcional a la de la distribución de <em>la posterior</em>. Esta es una propiedad importante que aprovecharemos en la siguiente sección. </p>
<p><em>Estimación del Valor Esperado de los Parámetros</em></p>
<p>Calcularemos numéricamente los valores esperados de $\beta$ y $\sigma_\epsilon$ ($\mu_\beta$ y $\mu_{\sigma_\epsilon}$ respectivamente) para compararlos con los valores estimados con las muestras de <code>stan</code>.</p>
<p>$$
\mu_\beta = \int_0^{+\infty} \int_{-\infty}^{+\infty} \beta \cdot \pi(\beta, \sigma_\epsilon \mid y) ~\textrm{d}\beta~\textrm{d}\sigma_\epsilon
$$</p>
<div class="highlight-r"><pre class="hl">stan_beta_mu <span class="hl opt"><-</span> stan_summary<span class="hl opt">[</span><span class="hl str">"beta"</span><span class="hl opt">,</span> <span class="hl str">"mean"</span><span class="hl opt">]</span>
stan_sigma_mu <span class="hl opt"><-</span> stan_summary<span class="hl opt">[</span><span class="hl str">"sigma"</span><span class="hl opt">,</span> <span class="hl str">"mean"</span><span class="hl opt">]</span>
stan_beta_mu_se <span class="hl opt"><-</span> stan_summary<span class="hl opt">[</span><span class="hl str">"beta"</span><span class="hl opt">,</span> <span class="hl str">"se_mean"</span><span class="hl opt">]</span>
stan_sigma_mu_se <span class="hl opt"><-</span> stan_summary<span class="hl opt">[</span><span class="hl str">"sigma"</span><span class="hl opt">,</span> <span class="hl str">"se_mean"</span><span class="hl opt">]</span>
beta_mu <span class="hl opt"><-</span> <span class="hl kwd">integrate</span><span class="hl opt">(</span><span class="hl kwa">function</span><span class="hl opt">(</span>sigmas<span class="hl opt">) {</span>
<span class="hl kwd">sapply</span><span class="hl opt">(</span>sigmas<span class="hl opt">,</span> <span class="hl kwa">function</span><span class="hl opt">(</span>sigma<span class="hl opt">) {</span>
<span class="hl kwd">integrate</span><span class="hl opt">(</span>
<span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">)</span> beta <span class="hl opt">*</span> <span class="hl kwd">posterior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">F</span><span class="hl opt">),</span>
<span class="hl opt">-</span><span class="hl kwb">Inf</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span>
<span class="hl opt">)</span>$value
<span class="hl opt">})</span>
<span class="hl opt">},</span> <span class="hl num">0</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span><span class="hl opt">)</span>$value
sigma_mu <span class="hl opt"><-</span> <span class="hl kwd">integrate</span><span class="hl opt">(</span><span class="hl kwa">function</span><span class="hl opt">(</span>sigmas<span class="hl opt">) {</span>
<span class="hl kwd">sapply</span><span class="hl opt">(</span>sigmas<span class="hl opt">,</span> <span class="hl kwa">function</span><span class="hl opt">(</span>sigma<span class="hl opt">) {</span>
<span class="hl kwd">integrate</span><span class="hl opt">(</span>
<span class="hl kwa">function</span><span class="hl opt">(</span>beta<span class="hl opt">)</span> sigma <span class="hl opt">*</span> <span class="hl kwd">posterior_density</span><span class="hl opt">(</span>beta<span class="hl opt">,</span> sigma<span class="hl opt">,</span> data<span class="hl opt">,</span> log <span class="hl opt">=</span> <span class="hl kwb">F</span><span class="hl opt">),</span>
<span class="hl opt">-</span><span class="hl kwb">Inf</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span>
<span class="hl opt">)</span>$value
<span class="hl opt">})</span>
<span class="hl opt">},</span> <span class="hl num">0</span><span class="hl opt">, +</span><span class="hl kwb">Inf</span><span class="hl opt">)</span>$value
<span class="hl opt">(</span>stan_beta_mu <span class="hl opt">-</span> beta_mu<span class="hl opt">) /</span> stan_beta_mu_se
<span class="hl opt">(</span>stan_sigma_mu <span class="hl opt">-</span> sigma_mu<span class="hl opt">) /</span> stan_sigma_mu_se
</pre></div>
<pre><code>[1] -0.1928841
[1] 0.4831648
</code></pre>
<p>La herramienta de integración numérica que usamos nos dice que esos estimados tienen al menos 4 dígitos decimales de precisión. Asimismo, como referencia, hemos ajustado el modelo con una regresión lineal clásica (columna <code>lm</code>). Los parámetros del modelo que usamos para simular esta data se muestran en la columna <code>source_value</code> y la posición de la moda está en <code>post_mode</code>.</p>
<pre><code>parameter stan_mean stan_se numeric_mean lm source_value post_mode
beta 9.3048 0.0069 9.3062 9.2898 8.9 9.3031
sigma 38.6470 0.1211 38.5885 36.6511 30.4 35.2620
</code></pre>
<p>Si no tuviéramos un estimado más preciso del valor de los parámetros, asumiríamos que el error absoluto debe de ser menor que dos errores estándar <code>stan_se</code> (al $97.7\%$), siguiendo una distribución normal. Comparados con nuestros estimados, los errores en los de <code>stan</code> son apenas $-0.19$ y $0.48$ errores estándar. Es decir, los errores <code>stan_se</code> son estimados conservadores, tal como lo afirma su documentación.</p>
<p>La media de empírica de nuestros parámetros, calculada con los valores de las muestras de <code>stan</code>, se aproxima a la esperanza de sus valores a medida que las muestras tienen una densidad proporcional a la de <em>la posterior</em> y corresponden a muestras independientes tomadas de la misma. El número efectivo de muestras independientes que hemos obtenido de stan está representado por <a href="http://www.odelama.com/sandbox/#stanEComp">la variable <code>n_eff</code> que aparece en sus estimados</a>.</p>
<p><em>Aplicación de la Distribución Posterior</em></p>
<p>La gran ventaja de tener la distribución posterior de los parámetros es que podemos derivar cualquier expresión que dependa de ellos. Por ejemplo, supongamos que queremos saber en nuestro modelo si es posible que $y$ para $x=24$ puede llegar a ser mayor que su valor para $x=28$ considerando tanto la distribución de muestreo de $y$ como la incertidumbre que tenemos de los parámetros.</p>
<div class="highlight-r"><pre class="hl">y_density <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">(</span>x<span class="hl opt">) {</span>
<span class="hl kwd">map2_dbl</span><span class="hl opt">(</span>
stan_samples$beta <span class="hl opt">*</span> x<span class="hl opt">,</span> stan_samples$sigma<span class="hl opt">,</span>
<span class="hl kwa">function</span><span class="hl opt">(</span>mu<span class="hl opt">,</span> sd<span class="hl opt">)</span> <span class="hl kwd">rnorm</span><span class="hl opt">(</span><span class="hl num">1</span><span class="hl opt">,</span> mu<span class="hl opt">,</span> sd<span class="hl opt">)</span>
<span class="hl opt">)</span>
<span class="hl opt">}</span>
y_22 <span class="hl opt"><-</span> <span class="hl kwd">y_density</span><span class="hl opt">(</span>x<span class="hl opt">=</span><span class="hl num">22</span><span class="hl opt">)</span>
y_28 <span class="hl opt"><-</span> <span class="hl kwd">y_density</span><span class="hl opt">(</span>x<span class="hl opt">=</span><span class="hl num">28</span><span class="hl opt">)</span>
<span class="hl kwd">mean</span><span class="hl opt">(</span>y_22 <span class="hl opt">></span> y_28<span class="hl opt">)</span>
</pre></div>
<pre><code>[1] 0.1532
</code></pre>
<p>La respuesta es, sí es posible, y puede ocurrir en el $15\%$ de los casos. ¿Cómo calcularía algo así con una estimación puntual clásica?</p>
<p><em>Conclusión</em></p>
<p>Hemos necesitado usar un ejemplo trivial, solo para que nos permita mostrar en un solo gráfico la distribución <em>posterior</em> y podamos tener la noción de que lo que hace <code>stan</code> es exactamente lo que hemos descrito aquí. Sin embargo, nosotros hemos "hallado" la distribución con una cuadrícula de $100 \times 100$ puntos, pero no hemos muestreado los parámetros de allí, los hemos calculado mediante integrales evaluadas numéricamente. Más adelante resolveremos un modelo con $227$ parámetros, y un análisis con cuadrícula sería imposible, requería $100^{227}=10^{454}$ puntos y tendríamos que calcular integrales múltiples con $227$ niveles (orden) y aun si llegáramos a poder plantear algo así ¿cuánto tiempo de ejecución se requería?.</p>
<blockquote>
<p><em>Stan</em> es una plataforma de alto rendimiento para el modelado estadístico, que obtiene muestras de los valores de los parámetros, tomadas de la distribución <em>posterior</em> del modelo.</p>
</blockquote>
<h1 id="toyvariandosolointercepcinengruposcontinuacin">Toy | Variando solo Intercepción en Grupos (continuación)</h1>
<p>En las primeras secciones ajustamos el modelo "variando solo intercepción por grupo" a la data del <a href="http://www.odelama.com/sandbox/#toyModel">conjunto de datos <em>Toy</em></a>. Aquí revisamos lo que hicimos y continuamos el análisis.</p>
<h2 id="distribucinposteriorobtenidaconstan">Distribución Posterior obtenida con Stan</h2>
<p>En <a href="http://www.odelama.com/sandbox/#toyFitIntc">esta sección anterior ajustamos el modelo a la data</a> usando la función <code>brm()</code> y obtuvimos <a href="http://www.odelama.com/sandbox/#toyIntcOutput">este listado</a>. La función <code>brm()</code> transformó nuestra formula <code>y ~ 1 + x + (1 | g)</code> en un modelo en <a href="https://mc-stan.org/">código <em>stan</em></a>, incluyendo en él las distribuciones <em>anteriores</em> por defecto de <code>brms</code>. Con ése código y nuestra data, <code>brm()</code> llamó a <code>stan</code> para que "ajuste" el modelo a la data. La función <code>brm()</code> recibió de <code>stan</code> $5000$ muestras de la distribución <em>posterior</em> de los parámetros de nuestro modelo \eqref{eq:rmnIntercEx1} y en base a estas calculó los parámetros reportados en el listado mencionado. </p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/brm2_intc_dens_pars.png">
<img alt="Distribución posterior de $\beta$" class="pic-border" data-width="700" src="http://www.odelama.com/sandbox/brm2_intc_dens_pars.png" title="Distribución posterior de $\beta$" />
</a>
<p class="caption-text">Distribución posterior de los coeficientes $\beta_0$ y $\beta_1$ del modelo en \eqref{eq:rmnIntercEx1}.</p>
</div>
<p>La cantidad de muestras que se obtienen depende de los parámetros <code>chains</code>, <code>iter</code> y <code>warmup</code> en el llamado a <code>brm()</code>: $nsamples=chains \times (iter-warmup)$, en nuestro caso <code>=4*(2500-1250)=5000</code>. </p>
<p>Siempre existe la posibilidad de que <em>stan</em> no pueda recorrer satisfactoriamente todo el espacio muestral de los parámetros, y que por ello las muestras que retorne no sean representativas de la verdadera distribución <em>posterior</em>. Para mantenernos al tanto de eso <code>stan</code> emite estadísticas como <code>Rhat</code> y <code>Bulk_ESS</code> que aparecen en el sumario (para más detalles vea <a href="https://mc-stan.org/misc/warnings.html#bulk-ess">Brief Guide to Stan’s Warnings</a>). Lo deseable es que <code>Rhat</code> sea menor que <code>1.1</code>, siendo <code>1</code> su mejor valor posible. <code>Tail_ESS</code> y <code>Tail_ESS</code> deben ser mayor que $100 \times chains$; en nuestro ajuste usamos <code>iter = 2500</code> pues con su valor por defecto <code>2000</code> no pasábamos estos chequeos, ahora todo está bien ✅.</p>
<h2 id="estimacionesdelosparmetros">Estimaciones de los Parámetros</h2>
<p>El paquete R <code>brms</code> nos permite obtener la distribución de cada parámetro, pero también permite obtener directamente una estimación puntual de cada parámetro del modelo, que por simplicidad es lo que vamos a usar en este ejemplo.</p>
<p>Los coeficientes a nivel de población <a href="http://www.odelama.com/sandbox/#toyIntcOutput">están en el sumario</a>, pero no aquellos a nivel de grupo, como $\alpha_0^{(j)}$ en nuestro modelo. Con la función <code>ranef()</code> (<strong>ran</strong>dom <strong>ef</strong>fects: <em>efectos aleatorios</em>) se puede obtener el coeficiente para cada grupo, y con <code>fixef()</code> (<strong>fix</strong>ed <strong>ef</strong>fects: <em>efectos fijos</em> ) los coeficientes a nivel población, mientras que los estimados consolidando ambos orígenes se obtienen con <code>coef()</code>. </p>
<p>Por ejemplo, en nuestro ajuste la intercepción a nivel población es $9.36$ (<code>fixef()</code>: $\beta_0$) y la intercepción de los grupos <code>1</code> y <code>2</code> es $-3.34$ y $-2.31$ (<code>ranef()</code>: $\alpha_0^{(1)}$ y $\alpha_0^{(2)}$), entonces la intercepción <em>consolidada</em> de estos dos grupos es $6.01$ y $7.04$ (<code>coef()</code>: $\alpha_0^{(j)} + \beta_0$), pues $6.01= 9.36-3.34$ y $7.04= 9.36-2.31$.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">coef</span><span class="hl opt">(</span>brm2_intc<span class="hl opt">)</span>
</pre></div>
<pre><code>$g
, , Intercept
Estimate Est.Error Q2.5 Q97.5
1 6.009167 0.8352279 4.409063 7.671721
2 7.046681 0.4624261 6.146496 7.931515
3 10.660634 0.5802230 9.518355 11.818091
:
, , x
Estimate Est.Error Q2.5 Q97.5
1 3.696783 0.0568952 3.585757 3.807682
2 3.696783 0.0568952 3.585757 3.807682
3 3.696783 0.0568952 3.585757 3.8076
:
[truncado por razones de espacio]
</code></pre>
<h1 id="correlacinintra-clases">Correlación Intra-Clases</h1>
<p>Una pregunta relevante en RMN es ¿cuanta variación hay entre los grupos comparada con la variación entre individuos de la población? Esta proporción es medida por la <em>correlación intraclases</em> <em>(ICC: intraclass correlation):</em></p>
<p>\begin{equation}
ICC = \frac{\sigma_\alpha^2}{\sigma_\alpha^2 + \sigma_{\epsilon}^2}
\label{eq:icc}
\end{equation}</p>
<p>ICC es la proporción de la varianza en la variable de resultado $y$ que puede ser explicada por la varianza entre grupos. En nuestro modelo $ICC=2.96^2/(2.96^2+2.07^2) = 67\%$. A mayor ICC más relevante es el uso de grupos (RMN) en el modelo.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/icc_levels.png">
<img alt="Niveles de ICC" class="pic-border" data-width="1116" src="http://www.odelama.com/sandbox/icc_levels.png" title="Niveles de ICC" />
</a>
<p class="caption-text">Niveles de ICC en la posición de puntos agrupados por colores.</p>
</div>
<p>En la imagen arriba hay nubes de puntos de 6 colores diferentes representando la pertenencia de cada punto a un grupo. La posición de los puntos y la de los centros de las nubes tienen varianza $\sigma_{\epsilon}^2$ y $\sigma_\alpha^2$ respectivamente. La ICC aumenta de izquierda a derecha (20%, 50% y 90%), ilustrando que a mayor ICC mayor es la varianza entre grupos que entre puntos.</p>
<h1 id="regresinmultinivelrequerida">Regresión Multinivel Requerida</h1>
<p>Intuitivamente, un mayor valor de $\sigma_\alpha$ respecto a $\sigma_{\epsilon}$ corresponde a una mayor relevancia en el uso de grupos, pero ¿A partir de cuál valor son requeridas las RMN? Cuando el índice de efecto del diseño (<em>design effect index</em>) es mayor que $1.1$ se recomienda el uso de RMN [<a href="http://www.odelama.com/sandbox/#Lai-1">Lai et al. 2015</a>]. Este <em>índice de efecto del diseño (deff)</em> se define como:</p>
<p>\begin{equation}
de\textit{ff} = 1+(n_j-1) \cdot ICC
\label{eq:deff}
\end{equation}</p>
<p>Donde ${n}_j$ es el numero promedio de observaciones por grupo, e ICC la correlación intra-clases. En nuestro ejemplo $de\textit{ff} = 1+(10.3-1) \times 0.67= 7.2$ que <em>justifica</em> el uso de una regresión multinivel <em>con grupos de</em> <code>g</code>.</p>
<p>Si simplificamos la ecuación \eqref{eq:deff}, condicionada a que $de\textit{ff}>1.1$, se encuentra que lo que se requiere es:</p>
<p>$$ \frac{\sigma_{\alpha}}{\sigma_\epsilon} > \frac{1}{\sqrt{10 \cdot n_j-11}} $$</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/whenRmnRequired.png">
<img alt="¿Cuándo se requiere una RMN?" class="pic-border" data-width="800" src="http://www.odelama.com/sandbox/whenRmnRequired.png" title="¿Cuándo se requiere una RMN?" />
</a>
<p class="caption-text"><em>¿Cuándo se requiere una RMN?</em>. Cuando la proporción de $\sigma_\alpha$ entre $\sigma_\epsilon$ está arriba de la curva correspondiente a $\textrm{deff}>1.1$.</p>
</div>
<p>Que para $n_j$ entre 5 y 10 requiere que $\sigma_\alpha$ sea mayor que entre 0.16 y 0.11 veces el valor de $\sigma_{\epsilon}$, respectivamente. Lo cual no es una valla muy difícil de pasar.</p>
<h1 id="toyvariandointercepcin:bondaddelajuste">Toy | Variando Intercepción: Bondad del Ajuste</h1>
<p>Hay múltiples y complementarias formas de evaluar la calidad de un modelo ajustado con <em>stan</em>. Aquí veremos un par de ellas para tener idea de cómo progresan nuestros ejemplos.</p>
<h2 id="chequeopredictivodelaposteriorspanidpp-checkspan">Chequeo Predictivo de la <em>Posterior</em><span id="pp-check"></span></h2>
<p>Cada conjunto de parámetros correspondiente a una misma muestra de la <em>posterior</em> de nuestro modelo corresponde a una <em>"solución potencial"</em> del mismo. Con un conjunto de estos parámetros, junto con $X$, podemos replicar los valores de $y$, que denominaremos $\tilde{y}$, <em>prediciendo</em> sus valores según esta <em>solución potencial</em>. Si hacemos estas replicaciones con varias de estas <em>soluciones</em> tomadas al azar de la posterior, obtendremos réplicas $\tilde{y}$ que en conjunto describen tanto la incertidumbre que tenemos de los parámetros, así como la proveniente de la distribución de muestreo en el modelo. A la larga, lo que deberíamos esperar —si nuestro modelo es correcto— es que la distribución de estas réplicas $\tilde{y}$, denominada distribución <em>predictiva (de la) posterior</em>, sea similar a la de la $y$ observada en la data. El chequear cuánto se cumple o no esto se denomina <em>chequeo predictivo (de la) posterior</em>.</p>
<p>La distribución de estas réplicas $\pi(\tilde{y} \mid \tilde{x}, x, y)$ en la <a href="http://www.odelama.com/sandbox/#notaSimpli">notación simplificada</a> es $\pi(\tilde{y} \mid y)$, la cual obtenemos según:</p>
<p>\begin{align}
\pi(\tilde{y} \mid y) &= \int \pi(\tilde{y}, \theta \mid y) \, \textrm{d}\theta \nonumber \\
&= \int \pi(\tilde{y} \mid \theta) \cdot \pi(\theta \mid y)
\, \textrm{d}\theta \label{eq:posPredFactors} \\
\end{align}</p>
<p>El primer factor en \eqref{eq:posPredFactors} corresponde a la incertidumbre por la distribución de muestreo, y el segundo a la de los parámetros del modelo. En la práctica obtenemos $\tilde{y}$ según lo describimos inicialmente:</p>
<p>\begin{align}
\theta^{(m)} \sim p(\theta \mid y) \label{eq:sampleTheta} \\
\tilde{y}^{(m)} \sim p(y \mid \theta^{(m)}) \label{eq:sampleY}
\end{align}</p>
<p>En \eqref{eq:sampleTheta} tomamos muestras de $\theta^{(m)}$ de <em>la posterior</em> y en \eqref{eq:sampleY} usamos cada muestra para muestrear a su vez valores de $\tilde{y}^{(m)}$ según su distribución de muestreo. Note cómo estos pasos corresponden a la incertidumbre de los parámetros \eqref{eq:sampleTheta} y de la distribución de muestreo \eqref{eq:sampleY}.</p>
<p>El paquete R <code>bayesplot</code>, desarrollado por el mismo equipo de desarrolladores de <em>stan</em>, incluye decenas de tipos de gráficos para comparar las réplicas $\tilde{y}$ con la data observada $y$ (<a href="https://mc-stan.org/bayesplot/reference/PPC-overview.html">vea más detalles aquí</a>). Afortunadamente, <code>brms</code> hace fácil acceder a esos gráficos con su función <code>pp_check()</code>.</p>
<p>Con la función <code>pp_check()</code> (<a href="https://mc-stan.org/bayesplot/reference/PPC-overview.html"><em>posterior predictive checking</em></a>) se obtienen diversos gráficos variando su parámetro <code>type</code>. Por defecto esta función muestra en color celeste la DPP ($y_{rep}$: $y$ <em>replicada</em>) y con azul oscuro la distribución en la data observada $y$. </p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">pp_check</span><span class="hl opt">(</span>brm2_intc_slo<span class="hl opt">,</span> nsamples<span class="hl opt">=</span><span class="hl num">100</span><span class="hl opt">) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>title<span class="hl opt">=</span><span class="hl str">"Toy | Variando solo Intercepción y Pendiente por grupo"</span><span class="hl opt">) +</span>
<span class="hl kwd">theme</span><span class="hl opt">(</span>legend<span class="hl opt">.</span>position <span class="hl opt">=</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">0.8</span><span class="hl opt">,</span> <span class="hl num">0.5</span><span class="hl opt">))</span>
</pre></div>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/pp_check_only_intc.png">
<img alt="Bondad de modelo Solo Intercepción" class="pic-border" data-width="600" src="http://www.odelama.com/sandbox/pp_check_only_intc.png" title="Bondad de modelo Solo Intercepción" />
</a>
<p class="caption-text"><em>Distribución Predictiva Posterior del Modelo Toy Variando solo Intercepción por Grupo</em>. Nuestro modelo se aproxima, pero no lo suficiente a la data observada.</p>
</div>
<p>Si bien no queremos sobre ajustar el modelo a la data, esperamos que siga su patrón general. En el gráfico vemos que nuestro modelo es deficiente en varios tramos, como en las colas.</p>
<h2 id="criteriodeinformacinloo">Criterio de Información Loo</h2>
<p>A veces se requiere un índice que cuantifique la bondad del ajuste. En nuestros ejemplos <a href="http://mc-stan.org/loo/">usaremos <code>loo()</code></a>. </p>
<p><a href="https://mc-stan.org/loo/"><code>loo</code></a> y WAIC (<a href="https://en.wikipedia.org/wiki/Watanabe%E2%80%93Akaike_information_criterion"><em>widely applicable information criterion</em></a>) son métodos similares —con los mismos supuestos— para estimar la capacidad predictiva fuera-de-la-muestra de un modelo, esto en contraposición a estimaciones de que tan bien el modelo se ajustó a la data, como lo hacen índices del tipo R2. </p>
<p>Se recomienda preferir <code>loo</code> sobre WAIC porque tiene la ventaja de poder darse cuenta cuándo su cálculo ha encontrado dificultades que no hacen confiable su resultado. Esta basado en la <em>verosimilitud</em> de la capacidad predictiva del modelo, calculada usando validación cruzada del tipo deja-uno-afuera (LOO-CV: <em>Leave-one-out cross-validation</em>).</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">loo</span><span class="hl opt">(</span>brm2_intc<span class="hl opt">)</span>
</pre></div>
<pre><code> Computed from 5000 by 504 log-likelihood matrix
Estimate SE
elpd_loo -1111.3 19.2
p_loo 50.8 4.3
looic 2222.6 38.4
------
Monte Carlo SE of elpd_loo is 0.2.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 495 98.2% 533
(0.5, 0.7] (ok) 9 1.8% 103
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.
</code></pre>
<p>El índice <code>loo</code> es reportado de dos maneras (vea <a href="https://mc-stan.org/loo/reference/loo-glossary.html">más detalles aquí</a>): </p>
<ul>
<li><code>elpd_loo</code> es un valor que buscamos maximizar (mayor es mejor). </li>
<li><code>looic</code> es un valor que buscamos minimizar (menor es mejor).</li>
</ul>
<p><code>eplpd_loo</code> es el estimado del logaritmo de la verosimilitud predictiva de $y$ (<em>expected log predictive density</em>) y puede ser positivo o negativo. <code>looic</code> se calcula fácilmente a partir de <code>elpd_loo</code> pues $looic \equiv -2 \times elpd\_loo$.</p>
<p><code>p_loo</code> es un estimador de la probabilidad de que el modelo prediga valores correctos de $y$ en datos futuros. Asintóticamente —bajo ciertas condiciones— <code>p_loo</code> puede interpretarse como el número efectivo de parámetros. Cuando el modelo tiene un buen comportamiento se observa que $p\_loo < N$ y $p\_loo < p$, donde $p$ es el número total de parámetros en el modelo y $N$ el número de observaciones. Cuando ésto no se cumple, señala que el modelo tiene una capacidad predictiva muy débil y puede indicar un error de especificación grave en el modelo.</p>
<p><em>Valores de loo() para nuestro modelo</em></p>
<p>En nuestro modelo, tenemos $p=53$ parámetros (49 intercepciones de grupos + dos coeficientes a nivel población + dos varianzas) y $N=504$, entonces $p\_loo < 53$ y $p\_loo < 504$ por lo que <code>p_loo</code> no es señal de problemas. Los demás estimadores lucen bien respecto a su <code>SE</code> (error estándar) y los usaremos cuando comparemos este modelo con otro. El diagnóstico <code>Pareto</code> muestra que <code>loo</code> ha sido estimado de manera confiable.</p>
<h1 id="reversinaunaregresinsingrupos">Reversión a una Regresión sin Grupos</h1>
<p>Las ecuaciones \eqref{eq:rmnIntercF1}-\eqref{eq:rmnIntercF4}, en la notación de de nuestro modelo en \eqref{eq:rmnIntercEx1}, corresponden a</p>
<p>\begin{align}
y^{(i)} &= \alpha_0^{(j)} + \beta_0 + {\beta_1}{x}^{(i)} +\epsilon^{(i)}
,~para~j=1, \dots ,J \label{eq:simpToyIntc} \\
\alpha_0^{(j)} &\sim \mathcal{N}(0, \sigma_\alpha^2) \nonumber \\
\epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2) \nonumber
\end{align}</p>
<p>En la distribución de $\alpha_0^{(j)}$ podemos ver que conforme $\sigma_\alpha$ se aproxima a cero, entonces $\alpha_0^{(j)}$ se aproxima también a $0$:</p>
<p>$$ \lim_{\sigma_\alpha \to 0} \alpha_0^{(j)} \to 0$$</p>
<p>Es decir, cuando no hay varianza entre grupos, $\alpha_0^{(j)}$ en la ecuación \eqref{eq:simpToyIntc} se vuelve cero, y ésta revierte a</p>
<p>\begin{align*}
y^{(i)} &= \beta_0 + {\beta_1}{x}^{(i)} +\epsilon^{(i)} ,~para~j=1, \dots ,J \\
\epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2)
\end{align*}</p>
<p>una regresión sin grupos.</p>
<h1 id="modelovariandointercepcinypendiente">Modelo Variando Intercepción y Pendiente</h1>
<p>Cuando la intercepción y la pendiente varían por grupos, el modelo aparte de incluir el término $\alpha_1^{(j)} x$ en la parte que cambia por grupo, ahora incluye que $\alpha_0^{(j)}$ y $\alpha_1^{(j)}$ pertenecen a una distribución Normal bivariada con correlación $\rho$. Esta es una extensión natural al caso en el que solo la intercepción varía por grupos, donde $\rho$ no es una restricción adicional, pues claramente puede ser cero.</p>
<p>\begin{align*}
y^{(i)} &= (\alpha_0^{(j)} + \alpha_1^{(j)} x_1) + {X}^{(i)} {\beta} +\epsilon^{(i)} ,~para~i=1, \dots ,N \\
\\
\begin{pmatrix}
\alpha_0^{(j)} \\ \alpha_1^{(j)}
\end{pmatrix}
&\sim
\mathcal{N} \begin{pmatrix}
\begin{pmatrix}
0 \\
0
\end{pmatrix}
,
\begin{pmatrix}
\sigma_{\alpha_0}^2 & \rho~ \sigma_{\alpha_0} \sigma_{\alpha_1} \\
\rho~ \sigma_{\alpha_0} \sigma_{\alpha_1} & \sigma_{\alpha_1}^2
\end{pmatrix} <br />
\end{pmatrix} \label{eq:rmnMulti} \\
\lower 0.4em \epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2)
\end{align*}</p>
<p>Cuando $\beta$ también incluye intercepción y pendiente $\beta_0+\beta_1 x+\cdots$ al igual que los grupos, las estimaciones por grupos están centradas en cero, es decir $\mu_{\alpha_0} = 0$ y $\mu_{\alpha_1} = 0$.</p>
<p>Es importante notar que todos los pares de coeficientes $(\alpha_0^{(j)} , \alpha_1^{(j)})$ provendrán de una misma distribución. Es decir, el modelo "aprende" de toda la data para ajustar la distribución bivariada, y con ésta estima los coeficientes $\alpha_0^{(j)}$ y $\alpha_1^{(j)}$ considerando el patrón de las observaciones en cada grupo.</p>
<p>Es fácil intuir que si tuviéramos variación con otro predictor por grupo como $\alpha_2^{(j)} x_2$, entonces $(\alpha_0, \alpha_1, \alpha_2)$ tendrían un distribución normal <em>tri-variada</em> con una matriz de covarianza de tamaño $3 \times 3$, con 3 coeficientes de correlación: uno para cada par de predictores.</p>
<h2 id="toyvariandointercepcinypendiente:ajustedelmodelo">Toy | Variando Intercepción y Pendiente: Ajuste del modelo</h2>
<p>Esta vez ajustaremos, a la data Toy, un modelo con intercepción y pendiente variando por grupo. Simplemente modificamos la fórmula para usar <code>x</code> como predictor en las variaciones por grupo: <code>(1 + x | g)</code>.</p>
<div class="highlight-r"><pre class="hl">brm2_intc_slo <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>y <span class="hl opt">~</span> <span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">|</span> g<span class="hl opt">),</span> data<span class="hl opt">,</span> seed<span class="hl opt">=</span><span class="hl num">1234</span><span class="hl opt">)</span>
<span class="hl kwd">summary</span><span class="hl opt">(</span>brm2_intc_slo<span class="hl opt">)</span>
</pre></div>
<pre><code> Family: gaussian
Links: mu = identity; sigma = identity
Formula: y ~ x + (1 + x | g)
Data: data2 (Number of observations: 504)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~g (Number of levels: 49)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 1.31 0.24 0.86 1.82 1.00 1779 2269
sd(x) 0.90 0.11 0.71 1.14 1.00 1974 2983
cor(Intercept,x) 0.14 0.21 -0.26 0.54 1.01 594 1069
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 9.09 0.25 8.60 9.61 1.00 3394 2719
x 3.70 0.14 3.43 3.97 1.00 1391 1860
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.55 0.05 1.45 1.66 1.00 4885 3320
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
</code></pre>
<p>El modelo resultante es</p>
<p>\begin{align*}
y^{(i)} &= \alpha_0^{(j)} + \alpha_1^{(j)} + 9.09 + 3.70 x +\epsilon^{(i)} ,~para~i=1, \dots ,504 \\
\begin{pmatrix}
\alpha_0^{(j)} \\ \alpha_1^{(j)}
\end{pmatrix}
&\sim
\mathcal{N}
\begin{pmatrix}
\begin{pmatrix}
0 \\
0
\end{pmatrix}
,
\begin{pmatrix}
1.31^2 & 0.14 \times 1.31 \times 0.90 \\
0.14 \times 1.31 \times 0.90 & 0.90^2
\end{pmatrix} <br />
\end{pmatrix} \\
\epsilon^{(i)} &\sim \mathcal{N}(0, 1.55^2)
\end{align*}</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">print</span><span class="hl opt">(</span><span class="hl kwd">coef</span><span class="hl opt">(</span>brm2_intc_slo<span class="hl opt">),</span> digits <span class="hl opt">=</span> <span class="hl num">4</span><span class="hl opt">)</span>
</pre></div>
<pre><code>$g
, , Intercept
Estimate Est.Error Q2.5 Q97.5
1 8.388 0.7567 6.932 9.865
2 8.202 0.5326 7.154 9.215
3 8.440 0.8177 6.828 9.979
:
, , x
Estimate Est.Error Q2.5 Q97.5
1 2.464 0.2915 1.891 3.041
2 3.278 0.1615 2.960 3.598
3 4.377 0.2284 3.934 4.842
:
</code></pre>
<p>El error estándar <code>Est.Error</code> de la correlación $\rho$ es más alto que la estimación de la misma, es decir hay incertidumbre de que sea un valor positivo. Lo mismo señala el rango <code>l-95% CI</code> y <code>u-95% CI</code>, con los límites del <em>intervalo creíble</em> de la estimación [<a href="http://www.odelama.com/sandbox/#Burkner-1">Burkner-1</a>], pues contiene el cero.</p>
<p>Este <em>intervalo creíble</em> (bayesiano) es calculado de la distribución posterior de $\rho$ obtenida con <code>stan</code> y es una afirmación de que la probabilidad de que se $\rho$ encuentre en este intervalo es del 95%. En cambio, un <em>intervalo de confianza</em> (frecuentista) de este tipo dice que si obtenemos repetidas muestras, de la misma población de donde proviene la data, y estimamos el intervalo con el mismo método, en el 95% de los casos el intervalo va a contener el parámetro estimado, lo cual no es una afirmación sobre la certidumbre del valor del estimado como ocurre en el modo Bayesiano.</p>
<p>En general, estos intervalos se hacen más angostos con mayor volumen de data y más anchos con la varianza de ésta. En el caso de $\rho$, que por su naturaleza es pequeño, requiere un intervalo muy angosto para ser estadísticamente significativo, el valor <em>verdadero</em> de $\rho$ con el que se generó la data es de $0.4$. </p>
<h2 id="toyvariandointercepcinypendiente:bondaddelajuste">Toy | Variando Intercepción y Pendiente: Bondad del Ajuste</h2>
<p><em>pp_check</em>. Construimos la gráfica de la distribución posterior de $y$ como hicimos antes.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/brm2_intc_slope_dens_pars.png">
<img alt="Bondad de modelo variando Intercepción y Pendiente" class="pic-border" data-width="600" src="http://www.odelama.com/sandbox/brm2_intc_slope_dens_pars.png" title="Bondad de modelo variando Intercepción y Pendiente" />
</a>
<p class="caption-text"><em>Distribución Predictiva Posterior del Modelo Toy Variando Intercepción y Pendiente por Grupo</em>. El ajuste ahora está mucho mejor explicando el comportamiento observado en la data.</p>
</div>
<p>El ajuste ahora está mucho mejor explicando el comportamiento observado en la data.</p>
<p><em>loo()</em>. El indice <code>loo()</code> que se obtiene ahora es:</p>
<pre><code>Computed from 4000 by 504 log-likelihood matrix
Estimate SE
elpd_loo -977.9 16.3
p_loo 71.3 4.8
looic 1955.8 32.7
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 481 95.4% 502
(0.5, 0.7] (ok) 20 4.0% 100
(0.7, 1] (bad) 3 0.6% 73
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.
</code></pre>
<p>Vemos en <code>Pareto k diagnostic</code> que tenemos tres observaciones con mala clasificación y el estimado de loo podría no ser confiable. Para mayor seguridad, hacemos validación <code>k-fold</code> sobre nuestro último modelo y el anterior, y comparamos los resultados.</p>
<div class="highlight-r"><pre class="hl">brm2_intc_slo <span class="hl opt"><-</span><span class="hl kwd">add_model_metrics</span><span class="hl opt">(</span>brm2_intc_slo<span class="hl opt">)</span>
brm2_intc_slo_kfold <span class="hl opt"><-</span> <span class="hl kwd">kfold</span><span class="hl opt">(</span>brm2_intc_slo<span class="hl opt">,</span> K<span class="hl opt">=</span><span class="hl num">12</span><span class="hl opt">,</span> chains<span class="hl opt">=</span><span class="hl num">5</span><span class="hl opt">)</span>
brm2_intc_kfold <span class="hl opt"><-</span> <span class="hl kwd">kfold</span><span class="hl opt">(</span>brm2_intc<span class="hl opt">,</span> K<span class="hl opt">=</span><span class="hl num">12</span><span class="hl opt">,</span> chains<span class="hl opt">=</span><span class="hl num">5</span><span class="hl opt">)</span>
<span class="hl kwd">compare_ic</span><span class="hl opt">(</span>brm2_intc_slo_kfold<span class="hl opt">,</span> brm2_intc_kfold<span class="hl opt">)</span>
</pre></div>
<pre><code> KFOLDIC SE
brm2_intc_slo 1971.34 33.49
brm2_intc 2237.36 41.54
brm2_intc_slo - brm2_intc -266.02 40.82
</code></pre>
<p>Lo cual nos dice que los estimados <code>looic</code> que teníamos eran confiables. De todos modos, en la comparación vemos que este último modelo es definitivamente mejor que el anterior, pues la diferencia es en ese sentido y su magnitud es significativamente diferente que el error estándar de la misma.</p>
<h1 id="modeloverdaderodelconjuntodedatostoyspanidtoymodel">Modelo "Verdadero" del Conjunto de Datos Toy <span id="toyModel"></span></h1>
<p>En <a href="http://www.odelama.com/sandbox/#genToyData">este anexo</a> encontrará el código usado para generar la data "Toy", donde el modelo implementado es:</p>
<p>$$
\begin{aligned}
y^{(i)} &= \alpha_0^{(j)} + \alpha_1^{(j)} + 2.7 + 1.3 x +\epsilon^{(i)} ,~para~i=1, \dots ,504 \\
\begin{pmatrix}
\alpha_0^{(j)} \\ \alpha_1^{(j)}
\end{pmatrix}
&\sim
\mathcal{N}
\begin{pmatrix}
\begin{pmatrix}
6.2 \\
2.3
\end{pmatrix}
,
\begin{pmatrix}
1.4^2 & 0.4 \times 1.4 \times 0.90 \\
0.4 \times 1.4 \times 0.90 & 0.90^2
\end{pmatrix} <br />
\end{pmatrix} \\
\epsilon^{(i)} &\sim \mathcal{N}(0, 1.6^2)
\end{aligned}
$$</p>
<p>Como ya hemos comentado, al tener el modelo intercepción y pendiente tanto a nivel población como de grupo, la parte del modelo a nivel población ($\beta_0$ y $\beta_1$) obtiene la media de ($\alpha_0$, $\alpha_1$), y a nivel grupo las variaciones centradas en cero.</p>
<p>En nuestra data la media <em>"verdadera"</em> es $\mu_{\alpha_0}=2.7+6.2=8.9$ y $\mu_{\alpha_1}=1.3+2.3=3.6$ mientras que el modelo ha estimado $\beta_0=\mu_{\alpha_0}=9.09$ y $\beta_1=\mu_{\alpha_1}=3.70$, así que considerando las varianzas empleadas estamos razonablemente cerca de los valores <em>"verdaderos"</em>.</p>
<p>Note que el modelo de donde proviene la data <a href="https://en.wikipedia.org/wiki/Identifiability"><em>no es identificable</em></a>, porque la manera arbitraria como se ha partido la media, de $\alpha_0$ por ejemplo, en $2.7$ para población y $6.2$ por grupos, es una de las infinitas maneras arbitrarias de hacer tal partición y producir la misma data. La respuesta a esto de las RMN, usando la media en la población y variaciones en los grupos, es sencilla y elegante.</p>
<p>Si ajustamos un modelo sin intercepción ni pendiente a nivel de población <code>y ~ -1 + (1 + x | g)</code>, deberíamos obtener $\alpha_0^{(j)}$ y $\alpha_1^{(j)}$ centradas en sus medias:</p>
<div class="highlight-r"><pre class="hl">brm2_intc_slo_gr <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>y <span class="hl opt">~ -</span><span class="hl num">1</span> <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">|</span> g<span class="hl opt">),</span> data<span class="hl opt">,</span> seed<span class="hl opt">=</span><span class="hl num">1234</span><span class="hl opt">)</span>
<span class="hl slc"># Media estimada empiricamente de alpha_0</span>
<span class="hl kwd">ranef</span><span class="hl opt">(</span>brm2_intc_slo_gr<span class="hl opt">)</span>$g<span class="hl opt">[,</span><span class="hl str">"Estimate"</span><span class="hl opt">,</span><span class="hl str">"Intercept"</span><span class="hl opt">] %>%</span> <span class="hl kwd">mean</span><span class="hl opt">()</span>
<span class="hl slc"># Media estimada empiricamente de alpha_1</span>
<span class="hl kwd">ranef</span><span class="hl opt">(</span>brm2_intc_slo_gr<span class="hl opt">)</span>$g<span class="hl opt">[,</span><span class="hl str">"Estimate"</span><span class="hl opt">,</span><span class="hl str">"x"</span><span class="hl opt">] %>%</span> <span class="hl kwd">mean</span><span class="hl opt">()</span>
</pre></div>
<pre><code>[1] 9.128348
[1] 3.707252
</code></pre>
<p>Efectivamente, las medias estimadas con el promedio de $\alpha_0^{(j)}$ y $\alpha_1^{(j)}$, cualitativamente coinciden con lo que obtuvimos antes: $9.09$ y $3.70$.</p>
<h1 id="tiposderegresionesydeagrupamiento">Tipos de Regresiones y de Agrupamiento</h1>
<p><em>Agrupamiento:</em> <strong>Completo</strong> </p>
<p>$$\begin{aligned}
para~i &= 1, \ldots, N: \\
\bar{y} &= \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p <br />
\end{aligned}$$</p>
<p>Una regresión que ignora todo posible agrupamiento, y que por lo tanto considera homogéneos a todos los puntos miembros de una misma población, es denominada de <em>agrupamiento completo</em>. En este caso una única regresión sobre todos los puntos produce una intercepción y pendiente con la tendencia global de las observaciones (línea punteada en el gráfico abajo). Como el modelo no es informado de la tendencia en cada grupo, la regresión representa un sub-ajuste (<em>underfitting</em>) sobre la data: vea grupos 23 y 6 en el gráfico.</p>
<p><em>Agrupamiento:</em> <strong>Ninguno</strong></p>
<p>$$\begin{aligned}
para~j &= 1, \ldots, J: \\
\bar{y}^{(j)} &= \beta_0^{(j)} + \beta_1^{(j)} x_1 + \beta_2^{(j)} x_2 + \dots + \beta_p^{(j)} x_p
\end{aligned}$$</p>
<p>Una regresión individual sobre cada grupo de puntos, ignorando completamente a los demás, se denomina <em>sin ningún agrupamiento</em>. En este caso se producen tantas intercepciones y pendientes como grupos (con más de un punto) hay en la data (línea de guiones en el gráfico abajo). Como cada regresión solo es informada con el patrón en el grupo, e ignora la tendencia global, tiende a un sobreajuste (<em>overfitting</em>) cuando hay pocas observaciones en el grupo: grupos 12, 22 y 42. Además es incapaz de resolver grupos con solo una observación: grupo 49.</p>
<p>Cuando se hace una regresión con pendiente e intercepción variando por grupo, pero sin usar el modelo Multinivel, (con <code>lm()</code> por ejemplo) se obtienen también estas regresiones <em>sin ningún agrupamiento</em>:</p>
<div class="highlight-r"><pre class="hl">fit <span class="hl opt"><-</span> <span class="hl kwd">lm</span><span class="hl opt">(</span>y <span class="hl opt">~</span> <span class="hl num">1</span> <span class="hl opt">+</span> x <span class="hl opt">+</span> g <span class="hl opt">+</span> g<span class="hl opt">:</span>x<span class="hl opt">,</span> data<span class="hl opt">)</span>
</pre></div>
<p><em>Agrupamiento:</em> <strong>Parcial</strong> </p>
<p>$$\begin{aligned}
para~i &= 1, \ldots, N, ~j = 1, \ldots, J: \\
\bar{y}^{(i)} &= \alpha_0^{(j)} + \alpha_1^{(j)} x_1 + \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p
\end{aligned}$$</p>
<p>Una regresión Multinivel puede ser entendida como una generalización de los dos tipos de agrupamiento anteriores, en la que se da un <em>agrupamiento parcial</em>. La pendiente e intercepción de cada grupo provienen de una misma distribución, con una media aproximada a la de <em>agrupamiento completo</em>, pero cada coeficiente por grupo es estimado teniendo en cuenta la tendencia de las observaciones en el grupo. Es decir, el modelo es informado tanto por la tendencia global como la local en cada grupo (línea continua en el gráfico abajo).</p>
<p>Cuando hay pocas observaciones, la regresión del grupo tiende al patrón global (grupos 22, 42 y 49), y cuando hay más observaciones la regresión se ajusta más sobre los puntos en el grupo (grupo 23).</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/regrPorTipoAgrupa.png">
<img alt="Tipo de Regresión vs Tipo de Agrupamiento" class="pic-border" data-width="920" src="http://www.odelama.com/sandbox/regrPorTipoAgrupa.png" title="Tipo de Regresión vs Tipo de Agrupamiento" />
</a>
<p class="caption-text"><em>Tipo de Regresión vs Tipo de Agrupamiento.</em> <br />(A) <em>Agrupamiento: Completo. </em> La línea punteada gris es de una regresión sobre todos los puntos sin grupos. Por ignorar la tendencia de cada grupo, tiende a sub-ajuste (<em>underfitting</em>): grupos 23 y 6. <br />(B) <em>Agrupamiento: Ninguno. </em> La línea de guiones es de una regresión sobre cada grupo de puntos sin ninguna conexión (agrupamiento) con los demás. Tiende a un sobreajuste (<em>overfitting</em>) cuando hay pocas observaciones en el grupo: grupos 12 y 22. Además es incapaz de resolver grupos con una sola observación: grupo 49. <br />(C) <em>Agrupamiento: Parcial. </em> La línea continua es de una regresión <em>Multinivel</em> que causa un agrupamiento <em>parcial</em>, donde la pendiente y la intercepción varían entre los dos extremos anteriores. El ajuste sobre los puntos en el grupo es mayor cuando el grupo tiene más observaciones: grupo 23, mientras que es más hacia la tendencia global cuando hay pocos puntos en el grupo: grupos 22, 42 y 49.</p>
</div>
<h2 id="encogimientoenagrupamientoparcial">Encogimiento en Agrupamiento Parcial</h2>
<p>Para ver con más detalle el efecto del <em>agrupamiento parcial</em> en las RMN, compararemos en cada grupo, del modelo Toy, la pendiente obtenida con regresiones usando los tres tipos de agrupamiento que vimos en la sección anterior.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/slopePorTipoAgrupa.png">
<img alt="Pendientes por tipo de Agrupamiento" class="pic-border" data-width="920" src="http://www.odelama.com/sandbox/slopePorTipoAgrupa.png" title="Pendientes por tipo de Agrupamiento" />
</a>
<p class="caption-text"><em>Pendientes por tipo de Agrupamiento.</em> La pendiente de cada grupo está dispuesta a lo largo del eje horizontal de acuerdo con número de observaciones en cada grupo. En la mayoría de los casos las estimaciones Multinivel (punto azul) tienen un valor entre la estimación sin agrupamiento (punto rojo) y la de agrupamiento completo (trazo gris). Esto reduce el ruido en la estimación (intervalos más angostos) y mejora la estimación acercándola más a la pendiente <em>real</em>, que se muestra con un rombo verde.</p>
</div>
<p>El trazo horizontal gris grueso es la pendiente de la regresión con agrupamiento total, calculada con <code>stan</code> (inferencia Bayesiana) sobre todas las observaciones. Los puntos representan la pendiente estimada en cada grupo y las líneas verticales su <em>intervalo creíble</em> al 95%. Los grupos están dispuestos en el eje horizontal según su tamaño (número de observaciones), así el grupo 49 -en el extremo izquierdo- tiene una sola observación y el 9 -a la derecha- tiene 21. Como la data es fabricada por nosotros, conocemos el valor <em>"verdadero"</em> de la pendiente, que es mostrado con un rombo verde.</p>
<p>Los trazos en azul son de la regresión <em>Multinivel</em> (calculada con <code>brms</code>) y los de rojo de la regresión <em>sin ningún agrupamiento</em> (una regresión por grupo, calculada con <code>stan</code>). Para el caso de grupos con 3 o menos puntos, <code>stan</code> no pudo calcular la regresión por errores de diagnóstico causados por la escasa data. Es por eso por lo que los dos primeros grupos a la izquierda solo muestran estimaciones en azul de la regresión <em>Multinivel</em>, que sí es capaz de resolver estos casos. </p>
<p>En la mayoría de los grupos, las estimaciones <em>RMN</em> (punto azul) tienen un valor entre la estimación sin agrupamiento (punto rojo) y la de agrupamiento completo (trazo gris), lo cual se denomina <em>encogimiento (shrinkage)</em>. Este <em>encogimiento</em> (hacia la media) tiende a ser mayor cuando hay menos observaciones en el grupo, pues la mayor incertidumbre en grupos con pocos puntos es compensada teniendo más en cuenta el patrón global de todos los puntos. Esto se refleja en estimaciones con mayor certidumbre (intervalos más angostos), y en la mayoría de los casos, mejor aproximación al valor <em>verdadero</em>.</p>
<p><a href="http://www.odelama.com/sandbox/#Gelman-2">Gelman 2013</a> señala que este encogimiento permite comparar efectos en los grupos sin tener que preocuparse en seguir los procedimientos clásicos de comparaciones múltiples. Incluso propone construir RMNs cuando se requieran estas comparaciones:</p>
<blockquote>
<p>Proponemos la construcción de modelos multinivel en los entornos donde surgen múltiples comparaciones. Los modelos multinivel realizan agrupamientos parciales, mientras que los procedimientos clásicos suelen mantener estacionarios los centros de los intervalos, ajustando las comparaciones múltiples haciendo los intervalos más amplios [...]. Por lo tanto, los modelos multinivel abordan el problema de las comparaciones múltiples y también producen estimaciones más eficientes, especialmente en entornos con baja variación a nivel de grupo, que es donde las comparaciones múltiples son una preocupación particular.</p>
</blockquote>
<h1 id="anlisisdedatosdejuniorschoolproject">Análisis de Datos de "Junior School Project"</h1>
<p>Vamos a analizar los datos que Peter Mortimer usó en <a href="http://www.odelama.com/sandbox/#Mortimer">Un Estudio de Escuelas Primarias Eficaces</a>. Dicho conjunto de datos es conocido como "Junior School Project" (JSP). Los datos estan disponibles en el paquete R <code>faraway</code> con el nombre <code>jsp</code>.</p>
<p>El JSP fue un estudio longitudinal de alrededor de 2000 alumnos de 50 escuelas primarias elegidas al azar de las 636 escuelas de la <em>Inner London Education Authority (ILEA)</em> en 1980. Aquí contamos con un subconjunto de esos datos:</p>
<ul>
<li><em>school</em>: Id de la escuela {1,..., 50} (categórica).</li>
<li><em>class</em>: Id del salón {1,..., 4} (categórica).</li>
<li><em>student</em>: Id del alumno {1,..., 1402} (categórica).</li>
<li><em>gender</em>: Sexo del alumno {m, f} (categórica).</li>
<li><em>year</em>: Año de estudio en la escuela primaria {3, 4, 5} (categórica).</li>
<li><em>social</em>: Clasificación socioeconómica del padre del alumno (categórica). Los niveles son: <code>I</code>, <code>II</code>, <code>III nonmanual</code>, <code>III manual</code>, <code>IV</code>, <code>V</code>, <code>lterm unemp</code> (desempleado a largo plazo), <code>curr unemp</code> (actualmente desempleado), <code>father absent</code> (padre ausente).</li>
<li><em>raven</em>: Puntaje (entre 0 y 40) del alumno en <a href="https://en.wikipedia.org/wiki/Raven%27s_Progressive_Matrices">la prueba Raven</a>, que mide el razonamiento abstracto no-verbal. En varios estudios se ha encontrado esta prueba, <a href="https://openpsychometrics.org/info/wais-raven-correlation/">mediana o altamente correlacionada con el IQ del evaluado</a>. Esta prueba fue evaluada al inicio de los estudios del alumno.</li>
<li><em>math</em>: Puntaje (entre 0 y 40) del alumno, por año de estudios, en la prueba de matemáticas. </li>
<li><em>english</em>: Puntaje (entre 0 y 100) del alumno, por año de estudios, en la prueba de lenguaje (Inglés).</li>
</ul>
<p><em>Objetivo Principal</em></p>
<p>El objetivo principal del estudio es establecer si algunas escuelas fueron más eficaces que otras para promover el aprendizaje y el desarrollo de los alumnos, después de descontar las diferencias sociodemográficas de los alumnos como <code>gender</code> o <code>raven</code>.</p>
<h2 id="correccindesesgo">Corrección de Sesgo</h2>
<p>Los puntajes de matemáticas y Raven tienen distribuciones con fuerte sesgo (<em>skew</em>) negativo (cola larga a la izquierda) debido a un <em>"efecto de techo"</em> en el que la diferencia de puntaje entre los alumnos más competentes varía muy poco respecto a la habilidad medida. Asimismo, el rango de puntaje en <code>english</code> es diferente al de las otras pruebas. Para corregir estos inconvenientes hemos transformado las variable <code>math</code>, <code>english</code>, y <code>raven</code> usando <a href="https://cran.r-project.org/web/packages/bestNormalize/vignettes/bestNormalize.html">el paquete R <code>bestNormalize</code></a>, que analiza varias posibles transformaciones monotónicas (p.ej. Yeo-Johnson, Lambert) y aplica la más eficaz en cada caso. En nuestro caso, el método seleccionado fue <em>Ordered Quantile normalization</em>.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_raw_trans.png">
<img alt="JSP: Sesgo y Transformación de Corrección" class="pic-noborder" data-width="800" src="http://www.odelama.com/sandbox/jsp_raw_trans.png" title="JSP: Sesgo y Transformación de Corrección" />
</a>
<p class="caption-text"><em>JSP Corrección de Sesgo en los Puntajes.</em> (a) Histograma con los puntajes originales en <code>math</code>, y (b) luego de su normalización.</p>
</div>
<p>Esta normalización no parece tan arbitraria si se considera, en primer lugar, que ésta es una práctica que no es rara en instituciones educativas, y en segundo lugar que ellos y nosotros lo hacemos precisamente para corregir el <em>efecto techo</em> encontrado. El <em>efecto techo</em> es señal de que el examen con el que se ha evaluado al alumno no discrimina lo suficiente a los más destacados, y en consecuencia, a partir de cierto nivel de competencia en adelante están recibiendo puntajes similares. Lo ideal es agregar preguntas más difíciles a la prueba para discriminar mejor a los alumnos más competentes, pero en los casos <em>post facto</em> como éste, la normalización redistribuye las notas para corregir esto.</p>
<p>Las variables transformadas han sido además estandarizadas (centradas y escaladas) para tener el cero como media y una desviación estándar como unidad. Estas variables transformadas tienen nombres terminando en <code>_z</code> y los valores originales se conservan en variables que terminan en <code>_raw</code>. El gráfico arriba muestra <code>math_raw</code> (datos originales con el sesgo) y <code>math_z</code> (datos transformados y estandarizados sin sesgo).</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">str</span><span class="hl opt">(</span>jsp<span class="hl opt">)</span>
</pre></div>
<pre><code>'data.frame': 3236 obs. of 12 variables:
$ school : Factor w/ 49 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ student : Factor w/ 1192 levels "1","2","3","4",..: 1 1 1 2 2 3 3 3 4 4 ...
$ year : Factor w/ 3 levels "3","4","5": 1 2 3 1 2 1 2 3 1 2 ...
$ class : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$ gender : Factor w/ 2 levels "m","f": 2 2 2 1 1 1 1 1 1 1 ...
$ social : Factor w/ 9 levels "I","II","III nonmanual",..: 9 9 9 2 2 2 2 2 9 9 ...
$ raven_z : num -0.403 -0.403 -0.403 -1.588 -1.588 ...
$ math_z : num -0.528 -0.415 -0.528 -1.462 -1.825 ...
$ english_z : num 0.55 0.928 -0.424 -2.128 -1.336 ...
</code></pre>
<h2 id="influenciadelaescuelaenelprogresoacadmicodelalumno">Influencia de la Escuela en el Progreso Académico del Alumno</h2>
<p>En el gráfico a continuación vemos una comparación entre las notas de matemáticas, al quinto año respecto al tercero, para algunas escuelas seleccionadas. Las líneas corresponden a regresiones lineales clásicas ajustadas sobre las notas de alumnos en una misma escuela, mostrando la tendencia.</p>
<div class="imcap-noframe">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_math_scores.png">
<img alt="Ejemplos de Notas de Matemáticas" class="pic-noborder" data-width="500" src="http://www.odelama.com/sandbox/jsp_math_scores.png" title="Ejemplos de Notas de Matemáticas" />
</a>
</div>
<p>En algunas escuelas la regresión es una línea casi paralela, completamente por encima de la de otra escuela (p.ej. 36:celeste sobre 18:verde oliva), que corresponde a que la mayoría de los alumnos en la escuela arriba, progresaron, del 3er al 5to año, más que los alumnos en la escuela abajo. Correspondientemente, la escuela arriba fue más eficaz en la mejora del alumno que la de abajo.</p>
<p>Cuando las líneas de tendencia entre dos escuelas no son paralelas, la que está más arriba en el lado izquierdo corresponde a la que más progreso ha tenido con los alumnos más rezagados, mientras la que está más arriba a la derecha es en la que más han mejorado los alumnos más avanzados. En este sentido, las líneas de poca pendiente corresponden a escuelas donde la brecha entre alumnos rezagados y avanzados ha disminuido.</p>
<p>Como tenemos los puntajes de <code>english</code> y <code>math</code> en el tercer y quinto año podríamos analizar el efecto de la escuela en el rendimiento por materia, pero terminaríamos con dos modelos y dos clasificaciones de las escuelas (una por cada materia). Aun cuando podríamos ver la manera de consolidar esos resultados en una sola clasificación, ésta parece una ruta complicada.</p>
<p><em>Fusionar Materias</em></p>
<p>Como hemos normalizado las notas de ambas materias, aquí hemos elegido unirlas y considerarlas simplemente puntajes al tercer <code>score_3_z</code> y quinto año <code>score_5_z</code>. Esto nos ayuda a mostrar el uso de RMN sin complicarnos en exceso con el análisis. Asimismo, siguiendo el mismo criterio, hemos conservado solo los alumnos que han estudiado del 3ro al 5to año en una misma escuela.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">str</span><span class="hl opt">(</span>jsp_sc<span class="hl opt">)</span>
</pre></div>
<pre><code>tibble [1,774 × 13] (S3: tbl_df/tbl/data.frame)
$ school : Factor w/ 49 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ student : Factor w/ 1192 levels "1","2","3","4",..: 1 3 4 6 7 8 11 14 15 16 ...
$ class : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
$ gender : Factor w/ 2 levels "m","f": 2 1 1 2 1 2 2 2 2 1 ...
$ social : Factor w/ 9 levels "I","II","III nonmanual",..: 9 2 9 4 9 9 4 2 7 5 ...
$ raven_z : num [1:1774] -0.403 -0.576 -1.817 -1.401 -1.256 ...
$ score_3_z : num [1:1774] -0.351 1.655 -0.215 -0.873 -0.485 ...
$ score_5_z : num [1:1774] -1.07056 1.83063 -0.00526 -2.20746 -0.78523 ...
</code></pre>
<p>Desde esta perspectiva, la variable dependiente "$y$" es <code>score_5_z</code>, y las variables predictoras son <code>score_3_z</code>, <code>raven_z</code>, <code>school</code>, <code>gender</code>, y <code>social</code>. Como es usual, la intercepción en la regresión tendrá las unidades de $y$, que debido al escalado que hicimos, es una desviación estándar de <code>score_5_z</code>, a la que para abreviar llamaremos $\sigma_y$. La intercepción corresponde al caso cuando <code>score_3_z=0</code> y <code>raven_z=0</code>, que debido al centrado de estas variables, es el puntaje esperado en 5to año cuando el alumno ha obtenido el puntaje promedio en <code>score_3_z</code> y <code>raven_z</code>, es decir la intercepción es el puntaje del <em>alumno promedio</em>.</p>
<h3 id="diagramacausalspanidcausalnetspan">Diagrama Causal <span id="causalNet"></span></h3>
<p>La hipótesis principal del estudio es que la escuela influye en cuanto aprende el alumno, lo cual se refleja en la mejora de su puntaje de 3ro a 5to año, es decir $\text{(score_5_z - score_3_z)} \sim \text{school} + \text{otros}$ donde $\sim$ significa <em>"depende de"</em>. Por lo tanto, $\text{score_5_z} \sim \text{score_3_z} + \text{school} + \text{otros}$.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_causal_net.png">
<img alt="JSP: Diagrama Causal" class="pic-noborder" data-width="750" src="http://www.odelama.com/sandbox/jsp_causal_net.png" title="JSP: Diagrama Causal" />
</a>
<p class="caption-text"><em>JSP Diagrama Causal</em>. Factores que influyen en el desempeño del alumno.</p>
</div>
<p>EL diagrama causal arriba considera que la escuela y la habilidad cognitiva del alumno influyen en sus notas, donde dicha habilidad cognitiva es una variable no observada pero estimada con el puntaje en la prueba Raven. </p>
<p>En varias escuelas, la cantidad de alumnos con padres en las diferentes clasificaciones sociales está bastante desbalanceada hacia las clases más altas o las más bajas, por lo que parece posible que la clase social del padre influya en la selección de la escuela en la que estudia el alumno. También es posible que esta clase del padre influya en el entorno y las facilidades que tiene para estudiar y en consecuencia influya en sus notas, lo cual se muestra con líneas tenues.</p>
<p>Finalmente, el género del alumno podría tener algún rol en el desempeño del alumno lo cual se muestra también con líneas tenues.</p>
<h2 id="modelodereferencia">Modelo de Referencia</h2>
<p>Para tener una noción del progreso en los siguientes modelos, vamos a construir el modelo más simple posible, en el que las notas de 5to año son constantes.</p>
<p>$$\begin{aligned}
y^{(i)} &\equiv \textrm{score_5_z}^{(i)} = \beta_0 + \sigma_\epsilon^{(i)} \\
\epsilon^{(i)} &\sim \mathcal{N}(0, \sigma_{\epsilon}^2)
\end{aligned}$$</p>
<div class="highlight-r"><pre class="hl">jsp_null <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>score_5_z <span class="hl opt">~</span> <span class="hl num">1</span><span class="hl opt">,</span> data<span class="hl opt">=</span>jsp_sc<span class="hl opt">)</span>
jsp_null
<span class="hl kwd">loo</span><span class="hl opt">(</span>jsp_null<span class="hl opt">)</span>
</pre></div>
<pre><code> Family: gaussian
Links: mu = identity; sigma = identity
Formula: score_5_z ~ 1
Data: jsp_sc (Number of observations: 1774)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.00 0.02 -0.05 0.05 1.00 3520 2341
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.00 0.02 0.97 1.03 1.00 3591 2725
Computed from 4000 by 1774 log-likelihood matrix
Estimate SE
elpd_loo -2518.7 27.2
p_loo 1.9 0.1
looic 5037.4 54.4
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
</code></pre>
<p>La intercepción $\beta_0=0.00$ es la media de $y \equiv \textrm{score_5_z}$, que es cero por efecto de su centrado. Recordemos que $\sigma_{\epsilon}$, representado por <code>sigma</code> en el listado arriba, tiene las mismas unidades que $y$. Entonces, como <code>sigma=1.00</code>, hay una desviación estándar $\sigma_{\epsilon}=1.00 \times \sigma_y$ en las notas de 5to año que no son explicadas por el modelo.</p>
<h2 id="intercepcinporescuelayporcategorasocial">Intercepción por Escuela y por Categoría Social</h2>
<p>Para evaluar si es apropiado usar grupos por escuela <code>school</code> y por clasificación social del padre <code>social</code> en nuestro modelo, lo ajustaremos a la data solo solo con intercepciones en estos grupos y con los predictores a nivel población que usaremos en nuestro modelo final.</p>
<div class="highlight-r"><pre class="hl">jsp_m0 <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>score_5_z <span class="hl opt">~</span> gender <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">|</span> school<span class="hl opt">) + (</span><span class="hl num">1</span> <span class="hl opt">|</span> social<span class="hl opt">),</span>
data<span class="hl opt">=</span>jsp_sc<span class="hl opt">,</span> refresh<span class="hl opt">=</span><span class="hl num">0</span><span class="hl opt">)</span>
</pre></div>
<pre><code>There were 127 divergent transitions after warmup.
The largest R-hat is 1.06, indicating chains have not mixed.
Effective Samples Size (ESS) is too low, indicating posterior means and medians
and tail quantiles may be unreliable.
</code></pre>
<p><code>stan</code> tiene dificultad para muestrear <em>la posterior</em>, esto se refleja en divergencias después del calentamiento (<a href="https://translate.google.com/translate?sl=en&tl=es&u=https%3A%2F%2Fmc-stan.org%2Fmisc%2Fwarnings.html%23divergent-transitions-after-warmup">lea acerca de eso aquí</a>). En nuestro caso se requieren más iteraciones <code>iter=3000</code> y más calentamiento <code>warmup=2000</code> (por defecto son 2000 y 1000 respectivamente), y hace falta un mayor valor para los parámetros <code>adapt_delta=0.999</code> y <code>max_treedepth=12</code> (por defecto son $0.8$ y $10$ respectivamente), use el vínculo sugerido anteriormente para entender mejor qué controlan estos parámetros.</p>
<div class="highlight-r"><pre class="hl">jsp_m0 <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>
score_5_z <span class="hl opt">~</span> gender <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">+ (</span><span class="hl num">1</span> <span class="hl opt">|</span> school<span class="hl opt">) + (</span><span class="hl num">1</span> <span class="hl opt">|</span> social<span class="hl opt">),</span>
data <span class="hl opt">=</span> jsp_sc<span class="hl opt">,</span>
seed <span class="hl opt">=</span> <span class="hl num">123</span><span class="hl opt">,</span>
iter <span class="hl opt">=</span> <span class="hl num">3000</span><span class="hl opt">,</span> warmup <span class="hl opt">=</span> <span class="hl num">2000</span><span class="hl opt">,</span>
control <span class="hl opt">=</span> <span class="hl kwd">list</span><span class="hl opt">(</span>adapt_delta <span class="hl opt">=</span> <span class="hl num">0.99</span><span class="hl opt">,</span> max_treedepth <span class="hl opt">=</span> <span class="hl num">12</span><span class="hl opt">)</span>
<span class="hl opt">)</span>
jsp_m0 <span class="hl opt"><-</span> <span class="hl kwd">add_criterion</span><span class="hl opt">(</span>jsp_m0<span class="hl opt">,</span> <span class="hl str">"loo"</span><span class="hl opt">)</span>
jsp_m0
<span class="hl kwd">loo</span><span class="hl opt">(</span>jsp_m0<span class="hl opt">)</span>
</pre></div>
<p><span id="jspIntcOnly"></span></p>
<pre><code> Family: gaussian
Links: mu = identity; sigma = identity
Formula: score_5_z ~ gender + score_3_z + raven_z + (1 | school) + (1 | social)
Data: jsp_sc (Number of observations: 1774)
Samples: 4 chains, each with iter = 3000; warmup = 2000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~school (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.24 0.03 0.19 0.31 1.00 1268 2030
~social (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.08 0.04 0.02 0.19 1.01 796 822
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.01 0.05 -0.12 0.09 1.00 1188 1977
genderf 0.04 0.03 -0.02 0.10 1.00 6771 3002
score_3_z 0.62 0.02 0.59 0.66 1.00 4194 3061
raven_z 0.17 0.02 0.13 0.21 1.00 4040 2532
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.62 0.01 0.60 0.64 1.00 6589 2717
</code></pre>
<p>Estas intercepciones por grupo son el efecto medio de las escuelas y la clase social en la nota de 5to año, con $\sigma_y$ como unidad. La intercepción a nivel población corresponde a la de los hombres, y <code>genderf</code> es la intercepción para las mujeres, donde la alumna promedio tiene $0.05\,\sigma_y$ mejor puntaje que el alumno promedio.</p>
<p>La varianza entre escuelas es $0.24^2$, y entre categorías sociales es $0.08^2$ que corresponden a un $ICC=0.24^2/(0.62^2+0.24^2+0.08^2)=12.8\%$ y $1.4\%$ respectivamente. Adicionalmente, la cantidad promedio de observaciones por escuela es $37$ y por categoría social es $197$, por lo que para las escuelas tenemos $\textit{deff}=1+(37-1)\times 0.128=5.6>1.1$ y para <code>social</code> $\textit{deff}=1+(197-1)\times 0.014=3.74>1.1$ por lo que es necesario mantener el modelando de la data con esos dos grupos.</p>
<blockquote>
<p>El efecto de las escuelas $\sigma_{school}=0.24$ es grande. Si comparamos una escuela en el $5\%$ superior $+1.6\,\sigma$ con una en el $5\%$ inferior $-1.6\,\sigma$ estamos hablando de una diferencia de $2 \times 1.6 \times 0.24\,\sigma_y=0.77\, \sigma_y$. Es decir, la nota del alumno promedio en una <em>muy buena</em> escuela está $0.77 \sigma_y$ por encima de la del mismo alumno en una escuela <em>muy mala</em>. Para poner esto en perspectiva, esta diferencia corresponde a $\text{score_3_z} = 0.97$ más $\text{raven_z} = 0.97$, pues $0.97 \times 0.62\,\sigma_y + 0.97 \times 0.17\sigma_y = 0.77\, \sigma_y$. Es decir, en este modelo, el alumno en una <em>muy mala</em> escuela necesita tener un <code>score_3_z</code> y <code>raven_z</code> casi con una desviación estándar sobre el promedio para alcanzar al alumno promedio en una <em>muy buena</em> escuela.</p>
</blockquote>
<h2 id="haciaelmodelofinal">Hacia el Modelo Final</h2>
<div class="highlight-r"><pre class="hl">formula <span class="hl opt"><-</span> score_5_z <span class="hl opt">~</span> <span class="hl num">1</span> <span class="hl opt">+</span> gender <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">+</span>
<span class="hl opt">(</span><span class="hl num">1</span> <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">|</span> school<span class="hl opt">) + (</span><span class="hl num">1</span> <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">|</span> social<span class="hl opt">)</span>
</pre></div>
<p>Exploratoriamente, nuestro modelo final propone que el puntaje al quinto año <code>score_5_z</code>, depende del género del alumno <code>gender</code>, y de su habilidad de razonamiento <code>raven</code>, pues esto ayuda a su aprendizaje. Además, propone que estas variables predictoras tienen influencia diferente en el alumno dependiendo de la escuela <code>(*|school)</code> y de la clasificación socioeconómica del padre <code>(*|social)</code>.</p>
<p>En este caso <code>gender</code> es una variable categórica a nivel población. No es una buena idea crear un grupo con esta variable pues tiene apenas dos valores posibles y técnicamente no es posible encontrar la varianza entre solo dos valores con un nivel de certeza razonable. Por otro lado, al tener nuestro modelo tanto intercepción como variables categóricas a nivel población, vamos a obtener la intercepción global y solo la intercepción para uno de los géneros (como vemos [en el modelo anterior(#jspIntcOnly)]). Esto ocurre porque no es identificable tener una intercepción global más una por cada género. De manera que vamos a remover la intercepción global para poder obtener la de cada género.</p>
<div class="highlight-r"><pre class="hl">formula <span class="hl opt"><-</span> score_5_z <span class="hl opt">~</span> <span class="hl num">0</span> <span class="hl opt">+</span> gender <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">+</span>
<span class="hl opt">(</span><span class="hl num">1</span> <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">|</span> school<span class="hl opt">) + (</span><span class="hl num">1</span> <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">|</span> social<span class="hl opt">)</span>
</pre></div>
<h3 id="asignacindedistribucionesanteriores">Asignación de Distribuciones Anteriores</h3>
<p>Esta vez vamos a establecer nosotros algunas distribuciones <em>anteriores</em> (<em>priors</em>) del modelo, en vez de usar las asignaciones por defecto como hemos venido haciendo. Para ver las distribuciones <em>anteriores</em> del modelo podemos usar la función <code>get_prior</code>.</p>
<div class="highlight-r"><pre class="hl"><span class="hl opt">(</span>jsp_priors <span class="hl opt"><-</span> <span class="hl kwd">get_prior</span><span class="hl opt">(</span>formula<span class="hl opt">,</span> jsp_sc<span class="hl opt">))</span>
</pre></div>
<pre><code> prior class coef group resp dpar nlpar bound
1 b
2 b genderf
3 b genderm
4 b raven_z
5 b score_3_z
6 lkj(1) cor
7 cor school
8 cor social
9 student_t(3, 0, 2.5) sd
10 sd school
11 sd Intercept school
12 sd raven_z school
13 sd score_3_z school
14 sd social
15 sd Intercept social
16 sd raven_z social
17 sd score_3_z social
18 student_t(3, 0, 2.5) sigma
</code></pre>
<p>En <a href="https://translate.google.com/translate?sl=en&tl=es&u=https%3A%2F%2Fmc-stan.org%2Frstanarm%2Farticles%2Fpriors.html">este artículo</a> encontrará sugerencias para estas distribuciones <em>anteriores</em>. En el listado arriba vemos que como hemos omitido la intercepción global, vamos a obtener los coeficientes <code>genderf</code> y <code>genderm</code>, correspondientes a la intercepción de cada género, tal como queríamos.</p>
<p>Los parámetros con la columna <code>prior</code> en blanco tienen una distribución <em>anterior</em> asignada por defecto por <code>stan</code>, y las que no están en blanco son las asignadas por defecto por <code>brms</code>. Es posible establecer una distribución <em>anterior</em> para cada parámetro por separado, o puede establecerse una misma distribución para un subconjunto de estos, como lo denotan las filas 1, 6, y 9. Así —por ejemplo— se puede apreciar que los priors con <code>class=b</code> son asignados por <code>stan</code> (están en blanco), mientras que los de las clases <code>class=cor</code>, <code>class=sd</code> y <code>class=sigma</code> tienen distribuciones por defecto de <code>brm</code> (filas 6, 9, 18). Más precisamente, la fila 6 indica que el <em>prior</em> con <code>class=cor</code> (la correlación entre coeficientes de un mismo grupo), para ambos grupos es <code>lkj(1)</code> (<a href="http://www.odelama.com/sandbox/#lkj">LKJ</a>).</p>
<p>La <code>class=b</code> corresponde a la distribución <em>anterior</em> de los coeficientes a nivel de población. Como nuestras variables estan estandarizadas, a esos coeficientes vamos a darles la distribución $\mathcal{N}(0,2)$.</p>
<div class="highlight-r"><pre class="hl">jsp_priors<span class="hl opt">[</span><span class="hl num">1</span><span class="hl opt">,</span> <span class="hl num">1</span><span class="hl opt">] <-</span> <span class="hl kwd">prior</span><span class="hl opt">(</span><span class="hl kwd">normal</span><span class="hl opt">(</span><span class="hl num">0</span><span class="hl opt">,</span> <span class="hl num">2</span><span class="hl opt">))</span>
</pre></div>
<p>Recordemos que el modelo va a estimar un juego de coeficientes para los predictores (incluyendo la intercepción) para cada miembro de un grupo, pero todos los coeficientes de un mismo grupo provienen de una misma distribución normal multivariada. Como nuestro modelo incluye dos grupos, tenemos dos distribuciones multivariadas.</p>
<p>Las <em>anteriores</em> con la clase <code>sd</code> y con la columna <code>group</code> no en blanco, corresponden a la desviación estándar de los predictores en cada grupo, como $\sigma_{\alpha_0}$ y $\sigma_{\alpha_1}$ en $\eqref{eq:rmnMulti}$. Para estos vamos a usar una distribución Exponencial(1).</p>
<div class="highlight-r"><pre class="hl">jsp_priors<span class="hl opt">[</span><span class="hl num">9</span><span class="hl opt">,</span> <span class="hl num">1</span><span class="hl opt">] <-</span> <span class="hl kwd">prior</span><span class="hl opt">(</span><span class="hl kwd">exponential</span><span class="hl opt">(</span><span class="hl num">1</span><span class="hl opt">))</span>
</pre></div>
<p>La <code>class=sigma</code> corresponde a la desviación estándar a nivel de población $\sigma_{\epsilon}$. Ya hemos visto que ésta es menor que 1, así que vamos a asignarle una distribución Cauchy(0, 1).</p>
<div class="highlight-r"><pre class="hl">jsp_priors<span class="hl opt">[</span><span class="hl num">18</span><span class="hl opt">,</span> <span class="hl num">1</span><span class="hl opt">] <-</span> <span class="hl kwd">prior</span><span class="hl opt">(</span><span class="hl kwd">cauchy</span><span class="hl opt">(</span><span class="hl num">0</span><span class="hl opt">,</span> <span class="hl num">0.8</span><span class="hl opt">))</span>
</pre></div>
<p>Finalmente, hay que notar que como una desviación estándar solo pueden tener valores positivos, cuando le asignamos una distribución con posibles valores negativos, como la Normal o Cauchy, internamente <code>brm</code> hace el ajuste del caso para que solo se use la parte positiva (truncada) de dicha distribución.</p>
<h3 id="chequeopredictivodeladistribucinanterior">Chequeo Predictivo de la Distribución Anterior</h3>
<p>Podemos revisar la bondad de nuestras distribuciones <em>anteriores</em> calculando directamente sus predicciones. Es decir, muestreamos parámetros de la distribución <em>anterior</em> y los introducimos en la fórmula del modelo (junto con los predictores en la data) para obtener $\tilde{y}_{\textrm{rep}}$ de la distribución de muestreo. Esto es análogo al <a href="http://www.odelama.com/sandbox/#pp-check">Chequeo Predictivo de la Posterior</a> que hicimos antes, solo que en lugar de usar la distribución de los parámetros condicionados a la data usamos sus <em>anteriores</em> (distribuciones no condicionadas a la data).</p>
<p>Podemos hacer esto facilmente con <code>brm()</code> usando el parámetro <code>sample_prior = "only"</code>; luego graficamos la distribución de los valores predecidos $y_{rep}$ la compararemos con el de la data usando <code>pp_check()</code>.</p>
<div class="highlight-r"><pre class="hl">jsp_m1_prior <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>
formula<span class="hl opt">,</span>
data <span class="hl opt">=</span> jsp_sc<span class="hl opt">,</span>
seed <span class="hl opt">=</span> <span class="hl num">123</span><span class="hl opt">,</span>
prior <span class="hl opt">=</span> jsp_priors<span class="hl opt">,</span>
sample_prior <span class="hl opt">=</span> <span class="hl str">"only"</span><span class="hl opt">,</span>
iter <span class="hl opt">=</span> <span class="hl num">1200</span><span class="hl opt">,</span> warmup <span class="hl opt">=</span> <span class="hl num">200</span><span class="hl opt">,</span> refresh <span class="hl opt">=</span> <span class="hl num">0</span>
<span class="hl opt">)</span>
<span class="hl kwd">pp_check</span><span class="hl opt">(</span>jsp_m1_prior<span class="hl opt">,</span> nsamples <span class="hl opt">=</span> <span class="hl num">70</span><span class="hl opt">) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>title<span class="hl opt">=</span><span class="hl str">"JSP: Chequeo de priors"</span><span class="hl opt">) +</span>
<span class="hl kwd">theme</span><span class="hl opt">(</span>legend<span class="hl opt">.</span>position <span class="hl opt">=</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">0.8</span><span class="hl opt">,</span> <span class="hl num">0.5</span><span class="hl opt">)) +</span>
<span class="hl kwd">coord_cartesian</span><span class="hl opt">(</span>xlim<span class="hl opt">=</span><span class="hl kwd">c</span><span class="hl opt">(-</span><span class="hl num">20</span><span class="hl opt">,</span><span class="hl num">20</span><span class="hl opt">))</span>
</pre></div>
<p>Lo que deseamos encontrar es que la $y_{rep}$ de los <em>anteriores</em> tenga una distribución aproximadamente proporcional a la de $y$ observada, pero dando posibilidad a un rango razonablemente más amplio de valores.</p>
<div class="imcap-noframe">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_check_prior.png">
<img alt="JSP: Chequeo de Priors" class="pic-noborder" data-width="500" src="http://www.odelama.com/sandbox/jsp_check_prior.png" title="JSP: Chequeo de Priors" />
</a>
</div>
<p>Las predicciones de nuestras distribuciones <em>anteriores</em> se ven bien ✅.</p>
<h2 id="ajustedelmodeloaladata">Ajuste del Modelo a la Data</h2>
<p>Habiendo confirmado que nuestras distribuciones <em>anteriores</em> son razonables, ajustamos el modelo a la data.</p>
<div class="highlight-r"><pre class="hl">jsp_m1 <span class="hl opt"><-</span> <span class="hl kwd">brm</span><span class="hl opt">(</span>
formula<span class="hl opt">,</span>
data <span class="hl opt">=</span> jsp_sc<span class="hl opt">,</span>
seed <span class="hl opt">=</span> <span class="hl num">123</span><span class="hl opt">,</span>
prior <span class="hl opt">=</span> jsp_priors<span class="hl opt">,</span>
iter <span class="hl opt">=</span> <span class="hl num">2500</span><span class="hl opt">,</span> warmup <span class="hl opt">=</span> <span class="hl num">1500</span><span class="hl opt">,</span>
control <span class="hl opt">=</span> <span class="hl kwd">list</span><span class="hl opt">(</span>adapt_delta <span class="hl opt">=</span> <span class="hl num">0.99</span><span class="hl opt">)</span>
<span class="hl opt">)</span>
jsp_m1 <span class="hl opt"><-</span> <span class="hl kwd">add_criterion</span><span class="hl opt">(</span>jsp_m1<span class="hl opt">,</span> <span class="hl str">"loo"</span><span class="hl opt">)</span>
</pre></div>
<pre><code> Family: gaussian
Links: mu = identity; sigma = identity
Formula: score_5_z ~ 0 + gender + score_3_z + raven_z +
(1 + score_3_z + raven_z | school) + (1 + score_3_z + raven_z | social)
Data: jsp_sc (Number of observations: 1774)
Samples: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~school (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.25 0.03 0.19 0.31 1.00 1297 2224
sd(score_3_z) 0.11 0.02 0.07 0.16 1.00 1924 2935
sd(raven_z) 0.04 0.03 0.00 0.09 1.00 1262 1960
cor(Intercept,score_3_z) -0.75 0.15 -0.96 -0.41 1.00 1566 2115
cor(Intercept,raven_z) 0.16 0.42 -0.71 0.89 1.00 4144 2621
cor(score_3_z,raven_z) -0.24 0.45 -0.90 0.73 1.00 3282 3136
~social (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.08 0.04 0.01 0.18 1.00 1013 975
sd(score_3_z) 0.03 0.02 0.00 0.08 1.00 2548 2027
sd(raven_z) 0.06 0.04 0.00 0.16 1.00 1167 1548
cor(Intercept,score_3_z) 0.01 0.49 -0.87 0.87 1.00 4384 2707
cor(Intercept,raven_z) -0.07 0.46 -0.88 0.79 1.00 2604 2395
cor(score_3_z,raven_z) -0.12 0.49 -0.92 0.81 1.00 1797 2993
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
genderm -0.01 0.05 -0.12 0.09 1.00 1145 1799
genderf 0.03 0.05 -0.07 0.14 1.00 1058 1673
score_3_z 0.62 0.03 0.56 0.68 1.00 2056 2236
raven_z 0.17 0.03 0.10 0.23 1.00 2341 1971
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.61 0.01 0.59 0.63 1.00 6454 2912
</code></pre>
<p>La varianza de <code>raven_z</code> entre escuelas y de <code>score_3_z</code> entre clasificaciones sociales parece demasiado pequeña. Probemos sin esos predictores en dichos grupos.</p>
<div class="highlight-r"><pre class="hl">jsp_m2 <span class="hl opt"><-</span> <span class="hl kwd">update</span><span class="hl opt">(</span>
jsp_m1<span class="hl opt">,</span>
formula<span class="hl opt">=</span> score_5_z <span class="hl opt">~</span> <span class="hl num">0</span> <span class="hl opt">+</span> gender <span class="hl opt">+</span> score_3_z <span class="hl opt">+</span> raven_z <span class="hl opt">+</span>
<span class="hl opt">(</span><span class="hl num">1</span> <span class="hl opt">+</span> score_3_z <span class="hl opt">|</span> school<span class="hl opt">) + (</span><span class="hl num">1</span> <span class="hl opt">+</span> raven_z <span class="hl opt">|</span> social<span class="hl opt">)</span>
<span class="hl opt">)</span>
jsp_m2 <span class="hl opt"><-</span> <span class="hl kwd">add_criterion</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span> <span class="hl str">"loo"</span><span class="hl opt">)</span>
<span class="hl kwd">loo</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span> jsp_m1<span class="hl opt">)</span>
</pre></div>
<pre><code>:
Model comparisons:
elpd_diff se_diff
jsp_m2 0.0 0.0
jsp_m1 -1.2 0.9
</code></pre>
<p>La diferencia entre el anterior modelo <code>jsp_m1</code> y este último <code>jsp_m2</code> no es significativa, pues la magnitud de la diferencia de $\text{elpd_loo}$ entre ellos "$\text{elpd_diff}=-1.2$" además de ser muy pequeña, apenas es $1.2/0.9=1.3$ veces mayor que su diferencia estándar, por lo que nos quedamos con el más simple: el último.</p>
<h2 id="modelofinal">Modelo Final</h2>
<p>Nuestro modelo final, incluyendo las distribuciones <em>anteriores</em>, es:</p>
<p>\begin{align}
\text{--- A nivel Alumno:}
\nonumber \\
\text{score_5_z} &= \beta_0^{(gender)}
+ \beta_1 \, \text{raven_z}
+ \beta_2 \, \text{score_3_z}
+ \pmb{\alpha}^{(school)}
+ \pmb{\gamma}^{(social)}
+ \epsilon \label{eq:jspFinal}
\\
\text{prior:} \nonumber \\
\epsilon &\sim \mathrm{Cauchy}(0, 1)
\nonumber \\[0.8em]
\text{--- A nivel Escuela:}
\nonumber \\
\pmb{\alpha}^{(school)} &= \alpha_0^{(school)}
+ \alpha_1^{(school)} \, \text{score_3_z} \label{eq:jspSchool}
\\
\begin{pmatrix}
\alpha_0^{(school)} \ \alpha_1^{(school)}
\end{pmatrix}
&\sim
\mathcal{N}
\begin{pmatrix}
\begin{pmatrix}
0 \
0
\end{pmatrix}
,
\Sigma_{school}
\end{pmatrix}
\nonumber \\
\text{priors:} \nonumber \\
\Sigma_{school} &= \operatorname{diag}(\boldsymbol{\mathbf{\tau}}_\alpha) \, \boldsymbol{\mathbf{\Omega}}_\alpha \operatorname{diag}(\boldsymbol{\mathbf{\tau}}_\alpha)
\nonumber \\
\boldsymbol{\mathbf{\tau}}_\alpha^{(j)} &\sim \mathrm{Exponential}(1), \; j = 1, 2 \
\nonumber \\
\boldsymbol{\mathbf{\Omega}}_\alpha &\sim \mathrm{LKJ}(1)
\nonumber \\[1em]
\text{--- A nivel Clase Social del Padre:} <br />
\nonumber \\
\pmb{\gamma}^{(social)} &= \gamma_0^{(social)}
+ \gamma_1^{(social)} \, \text{raven_z} \label{eq:jspSocial}
\\
\begin{pmatrix}
\gamma_0^{(social)} \ \gamma_1^{(social)}
\end{pmatrix}
&\sim
\mathcal{N}
\begin{pmatrix}
\begin{pmatrix}
0 \
0
\end{pmatrix}
,
\Sigma_{social}
\end{pmatrix} <br />
\nonumber \\
\text{priors:} \nonumber \\
\Sigma_{social} &= \operatorname{diag}(\boldsymbol{\mathbf{\tau}}_\gamma) \, \boldsymbol{\mathbf{\Omega}}_\gamma \operatorname{diag}(\boldsymbol{\mathbf{\tau}}_\gamma)
\nonumber \\
\boldsymbol{\mathbf{\tau}}_\gamma^{(j)} &\sim \mathrm{Exponential}(1), \; j = 1, 2 \
\nonumber \\
\boldsymbol{\mathbf{\Omega}}_\gamma &\sim \mathrm{LKJ}(1) \nonumber
\end{align}</p>
<div class="highlight-r"><pre class="hl">jsp_m2
<span class="hl kwd">loo</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">)</span>
</pre></div>
<p><span id="jsp_m2List"></span></p>
<pre><code> Family: gaussian
Links: mu = identity; sigma = identity
Formula: score_5_z ~ gender + score_3_z + raven_z + (1 + score_3_z | school) + (1 + raven_z | social) - 1
Data: jsp_sc (Number of observations: 1774)
Samples: 4 chains, each with iter = 2500; warmup = 1500; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~school (Number of levels: 48)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.24 0.03 0.19 0.31 1.00 1103 2138
sd(score_3_z) 0.10 0.02 0.06 0.15 1.00 1444 2434
cor(Intercept,score_3_z) -0.79 0.14 -0.99 -0.44 1.00 1512 2347
~social (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.08 0.04 0.02 0.18 1.00 1071 854
sd(raven_z) 0.06 0.04 0.00 0.15 1.01 967 1664
cor(Intercept,raven_z) -0.05 0.51 -0.92 0.85 1.00 2056 2278
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
genderm -0.01 0.06 -0.11 0.11 1.00 903 1899
genderf 0.03 0.05 -0.07 0.14 1.00 902 1794
score_3_z 0.62 0.03 0.57 0.67 1.00 1252 2679
raven_z 0.17 0.03 0.11 0.23 1.00 2421 2523
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.61 0.01 0.59 0.63 1.00 6072 2520
Computed from 4000 by 1774 log-likelihood matrix
Estimate SE
elpd_loo -1671.5 33.6
p_loo 66.3 3.1
looic 3342.9 67.2
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are ok (k < 0.7).
</code></pre>
<p>Nuestro modelo final explica el $1-0.61^2/1.00^2=63\%$ de la varianza de las notas de 5to año <code>score_5_z</code> ✅. Asimismo, $loo$ ha mejorado notoriamente.</p>
<p><em>Interpretación Preliminar</em></p>
<p>Una nota que debemos poner junto a cada interpretación que hagamos de los parámetros del modelo es (*) <em>"... en un modelo de regresión de las notas en primaria descontando el efecto de la escuela, el género, la categoría social del padre, y la capacidad de razonamiento abstracto del alumno."</em></p>
<p>La intercepción a nivel población, que es parte de la nota esperada del <em>alumno promedio</em>, está separada por género (<code>genderm</code>:masculino, <code>genderf</code>:femenino) y es muy pequeña, aunque con una ligera ventaja de las chicas sobre los chicos. </p>
<p>A nivel de escuelas, hay una correlación negativa $-0.79$ entre la intercepción y la pendiente de <code>score_3_z</code>, es decir a mayor intercepción menor pendiente y viceversa.</p>
<p>A nivel población, el coeficiente para <code>raven_z</code> es la correlación entre el puntaje en dicha prueba y la nota de 5to año, descontando <em>las demás variables en el modelo</em>, con un valor estimado en $17\%$. El coeficiente para <code>score_3_z</code> tiene una interpretación análoga, con un estimado en $62\%$. Si calculamos la correlación <em>bruta</em> entre <code>score_3_z</code> y <code>score_5_z</code> obtenemos $73\%$ Si consideramos nuestro [Diagrama Causal(#causalNet)], esto sugiere que $73\%-62\%=11\%$ proviene del efecto de las otras variables predictoras en dichas notas.</p>
<h3 id="chequeodelaposterior">Chequeo de la Posterior</h3>
<p>Demos una mirada a un par de chequeos gráficos para tener una intuición de que tan bien está el modelo.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_post_checks.png">
<img alt="JSP: Chequeo del Modelo Final" class="pic-noborder" data-width="800" src="http://www.odelama.com/sandbox/jsp_post_checks.png" title="JSP: Chequeo del Modelo Final" />
</a>
<p class="caption-text"><em>JSP: Chequeo del Modelo Final.</em> (a) Chequeo Predictivo Posterior del modelo, que luce bien. (b) Chequeo de Residuos vs Valores Ajustados ($\tilde{y}$ predicha), donde no hay patrones en la distribución de los puntos ni en su variabilidad. Los patrones diagonales son artefactos debido a que la data original es discreta.</p>
</div>
<p>El modelo luce bien en ambos graficos. En el primer gráfico vemos que la distribución de las predicciones del modelo contiene a la data observada. Aun cuando hay segmentos donde la data observada ésta en el margen de las predicciones, no se aleja demasiado. En el segundo no hay patrón en los residuos ni cambio en su variabilidad. Los patrones diagonales son artefactos causados por la data que es discreta.</p>
<h3 id="unaregresinporalumno">Una Regresión por Alumno</h3>
<p>Como los puntajes están conformados por las notas en <code>math</code> e <code>english</code>, podemos trazar para cada alumno una línea uniendo las notas en esas dos materias.</p>
<div class="highlight-r"><pre class="hl">step <span class="hl opt"><-</span> <span class="hl num">0.05</span>
sel_school <span class="hl opt"><-</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">29</span><span class="hl opt">,</span> <span class="hl num">22</span><span class="hl opt">,</span> <span class="hl num">1</span><span class="hl opt">,</span> <span class="hl num">14</span><span class="hl opt">,</span> <span class="hl num">41</span><span class="hl opt">,</span> <span class="hl num">38</span><span class="hl opt">,</span> <span class="hl num">36</span><span class="hl opt">,</span> <span class="hl num">31</span><span class="hl opt">)</span>
subset <span class="hl opt"><-</span> jsp_sc <span class="hl opt">%>%</span> <span class="hl kwd">filter</span><span class="hl opt">(</span>school <span class="hl opt">%</span><span class="hl kwa">in</span><span class="hl opt">%</span> sel_school<span class="hl opt">)</span>
subset <span class="hl opt"><-</span> subset <span class="hl opt">%>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>fake <span class="hl opt">=</span> <span class="hl kwb">F</span><span class="hl opt">) %>%</span>
<span class="hl kwd">bind_rows</span><span class="hl opt">(</span>
subset <span class="hl opt">%>%</span> <span class="hl kwd">mutate</span><span class="hl opt">(</span>score_3_z <span class="hl opt">=</span> score_3_z <span class="hl opt">-</span> step<span class="hl opt">,</span> fake <span class="hl opt">=</span> <span class="hl kwb">T</span><span class="hl opt">),</span>
subset <span class="hl opt">%>%</span> <span class="hl kwd">mutate</span><span class="hl opt">(</span>score_3_z <span class="hl opt">=</span> score_3_z <span class="hl opt">+</span> step<span class="hl opt">,</span> fake <span class="hl opt">=</span> <span class="hl kwb">T</span><span class="hl opt">)</span>
<span class="hl opt">) %>%</span> <span class="hl kwd">arrange</span><span class="hl opt">(</span>student<span class="hl opt">,</span> score_3_z<span class="hl opt">)</span>
preds <span class="hl opt"><-</span> <span class="hl kwd">posterior_predict</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span> subset<span class="hl opt">) %>%</span> <span class="hl kwd">colMeans</span><span class="hl opt">()</span>
subset <span class="hl opt"><-</span> subset <span class="hl opt">%>%</span> <span class="hl kwd">mutate</span><span class="hl opt">(</span>pred <span class="hl opt">=</span> preds<span class="hl opt">)</span>
subset <span class="hl opt">%>%</span>
<span class="hl kwd">ggplot</span><span class="hl opt">(</span><span class="hl kwd">aes</span><span class="hl opt">(</span>x <span class="hl opt">=</span> score_3_z<span class="hl opt">,</span> color <span class="hl opt">=</span> school<span class="hl opt">)) +</span>
<span class="hl kwd">geom_point</span><span class="hl opt">(</span>
data <span class="hl opt">= . %>%</span> <span class="hl kwd">filter</span><span class="hl opt">(</span>fake <span class="hl opt">==</span> <span class="hl kwb">F</span><span class="hl opt">),</span> alpha<span class="hl opt">=</span><span class="hl num">0.6</span><span class="hl opt">,</span>
<span class="hl kwd">aes</span><span class="hl opt">(</span>y <span class="hl opt">=</span> score_5_z<span class="hl opt">)</span>
<span class="hl opt">) +</span>
<span class="hl kwd">geom_line</span><span class="hl opt">(</span>
<span class="hl kwd">aes</span><span class="hl opt">(</span>y <span class="hl opt">=</span> pred<span class="hl opt">,</span> group <span class="hl opt">=</span> student<span class="hl opt">,</span> color <span class="hl opt">=</span> school<span class="hl opt">)</span>
<span class="hl opt">) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>
title <span class="hl opt">=</span> <span class="hl str">"Puntajes de Escuelas Seleccionadas"</span><span class="hl opt">,</span>
x <span class="hl opt">=</span> <span class="hl str">"Puntaje al 3er año"</span><span class="hl opt">,</span> y <span class="hl opt">=</span> <span class="hl str">"Puntaje al 5to año"</span>
<span class="hl opt">) +</span>
<span class="hl kwd">theme</span><span class="hl opt">(</span>legend<span class="hl opt">.</span>position<span class="hl opt">=</span><span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl num">0.9</span><span class="hl opt">,</span><span class="hl num">0.35</span><span class="hl opt">))</span>
</pre></div>
<p></p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_plot_final.png">
<img alt="JSP: Predicciones en escuelas seleccionadas" class="pic-noborder" data-width="800" src="http://www.odelama.com/sandbox/jsp_plot_final.png" title="JSP: Predicciones en escuelas seleccionadas" />
</a>
<p class="caption-text"><em>JSP: Predicciones en escuelas seleccionadas.</em> Se incluye solo las ocho escuelas con las intercepciones más altas y las más bajas. Segmentos de línea unen las predicciones de las notas para inglés y matemáticas de cada alumno. Estas líneas convergen en la esquina superior derecha del gráfico correspondiendo a que los alumnos más destacados obtienen puntajes elevados dependiendo menos de la escuela en la que estudian, en oposición a los alumnos más rezagados (en el extremo inferior izquierdo) que obtienen notas más diferentes dependiendo de la escuela.</p>
</div>
<p>En el gráfico <em>"JSP: Predicciones en escuelas seleccionadas"</em> arriba cada línea corresponde a un alumno, la línea une sus notas en 3er año con las de 5to, y se incluyen solo las ocho escuelas con las intercepciones más altas o bajas. Estas líneas convergen aproximadamente en la esquina superior derecha del gráfico. Esto corresponde a <a href="http://www.odelama.com/sandbox/#jsp_m2List">una correlación negativa, estimada</a> en $\text{cor(Intercept,\, score_3_z)}=-0.79$ entre la intercepción y pendiente por escuela, cuya interpretación es:</p>
<blockquote>
<p>Los alumnos más destacados obtienen puntajes elevados variando menos con la escuela en la que estudian, en oposición a los alumnos más rezagados cuyas notas varían más con la escuela. Cuando la escuela tiene un efecto negativo en el alumno (respecto al efecto medio de las escuelas), el impacto es mayor en los alumnos más rezagados.</p>
</blockquote>
<h2 id="hiptesisproposicinconparmetrosestimados">Hipótesis= Proposición con Parámetros Estimados</h2>
<p>A partir de la distribución conjunta de los parámetros del modelo podemos probar diferentes hipótesis de una manera que es imposible con estimaciones puntuales clásicas. Con hipótesis nos referimos a una proposición que depende de uno o más parámetros estimados por el modelo, como por ejemplo "$\beta_0>0$" o "$\alpha_0^{(2)}> \alpha_0^{(3)}$". A continuación, veremos algunos ejemplos.</p>
<p><span id="comp_schools"></span>
<em>Ejemplo Comparando la Intercepción de dos Escuelas</em></p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_demo_intc_schools.png">
<img alt="JSP: Comparación de Escuelas 29 y 6" class="pic-noborder" data-width="700" src="http://www.odelama.com/sandbox/jsp_demo_intc_schools.png" title="JSP: Comparación de Escuelas 29 y 6" />
</a>
<p class="caption-text"><em>JSP: Comparación de Intercepción de Escuelas s29 y s6.</em> Es contraintuitivo a primera vista del gráfico, pero la Intercepción de la escuela <code>s6</code> arriba es mayor que el de la escuela <code>s29</code> abajo en el 99% de las muestras. Se puede tener una intuición de cómo es esto posible si se imagina que cada muestra de la <code>s29</code> está apareada con una de la <code>s6</code> como lo sugieren las flechas verdes, donde la mayoría de las muestras de la escuela <code>s6</code> son mayores que las de la <code>s29</code>.</p>
</div>
<p>El gráfico arriba muestra la distribución de la intercepción de dos escuelas, donde en cada base se señala en azul su <em>intervalo creíble</em> al $95\%$ y su media. Como las distribuciones se superponen tanto, a primera vista parece contraintuitivo afirmar que tenemos una certeza del $99\%$ de que la intercepción de la escuela <code>s6</code> (arriba) es mayor la de la <code>s29</code> (abajo). Puede tener una intuición de cómo es esto posible si imagina que los valores de las intercepciones están apareadas de la manera como lo sugieren las flechas verdes.</p>
<p>En una inferencia Bayesiana, para comparar dos parámetros, lo hacemos usando sus valores en cada muestra y estimamos la certeza con la proporción de los hallazgos. Por ejemplo, de las $4000$ muestras obtenidas con nuestro modelo, en $3958$ de ellas la intercepción de la escuela <code>s6</code> es mayor que la de la <code>s29</code>, entonces la probabilidad estimada de esto es $3958/4000=98.95$.</p>
<p>Con <code>brms</code> se pueden probar hipótesis con la función <code>hypothesis()</code>. En nuestro ejemplo anterior la <em>hipótesis</em> es:</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">hypothesis</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span>
<span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl str">"s6>s29"</span> <span class="hl opt">=</span> <span class="hl str">"r_school[6,Intercept] > r_school[29,Intercept]"</span><span class="hl opt">),</span>
class <span class="hl opt">=</span> <span class="hl kwb">NULL</span><span class="hl opt">)</span>
</pre></div>
<pre><code>Hypothesis Tests for class :
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 s6>s29 0.35 0.16 0.1 0.62 94.24 0.99 *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
</code></pre>
<p>Donde la columna <code>Estimate</code> muestra el valor estimado de la diferencia entre las intercepciones de las escuelas ($s6-s29$). El valor en <code>Est.Error</code> es la desviación estándar de <code>Estimate</code>. La columna <code>Evid.Ratio</code> es la proporción de la evidencia a favor respecto a la que hay en contra de la hipótesis $Evid.Ratio= 0.9895/(1−0.9895)=94.24$. Finalmente, <code>Post.Prob</code> muestra la proporción de las muestras que están a favor de la hipótesis.</p>
<p>Para entender mejor los resultados de la función <code>hypothesis()</code>, graficamos la distribución de la diferencia de la intercepción de las escuelas (<code>s6</code> menos <code>s29</code>).</p>
<div class="highlight-r"><pre class="hl">school_diff <span class="hl opt"><-</span> as<span class="hl opt">.</span>data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">) %>%</span>
<span class="hl kwd">select</span><span class="hl opt">(</span>`r_school<span class="hl opt">[</span><span class="hl num">29</span><span class="hl opt">,</span>Intercept<span class="hl opt">]</span>`<span class="hl opt">,</span> `r_school<span class="hl opt">[</span><span class="hl num">6</span><span class="hl opt">,</span>Intercept<span class="hl opt">]</span>`<span class="hl opt">) %>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>intc_29 <span class="hl opt">=</span> `r_school<span class="hl opt">[</span><span class="hl num">29</span><span class="hl opt">,</span>Intercept<span class="hl opt">]</span>`<span class="hl opt">) %>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>intc_6 <span class="hl opt">=</span> `r_school<span class="hl opt">[</span><span class="hl num">6</span><span class="hl opt">,</span>Intercept<span class="hl opt">]</span>`<span class="hl opt">) %>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>diff_6_29 <span class="hl opt">=</span> intc_6 <span class="hl opt">-</span> intc_29<span class="hl opt">) %>%</span>
<span class="hl kwd">select</span><span class="hl opt">(</span>intc_29<span class="hl opt">,</span> intc_6<span class="hl opt">,</span> diff_6_29<span class="hl opt">)</span>
s6_greater_prob <span class="hl opt"><-</span> <span class="hl kwd">mean</span><span class="hl opt">(</span>school_diff$diff_6_29 <span class="hl opt">></span> <span class="hl num">0</span><span class="hl opt">)</span>
diff_sd <span class="hl opt"><-</span> <span class="hl kwd">sd</span><span class="hl opt">(</span>school_diff$diff_6_29<span class="hl opt">)</span>
diff_mean <span class="hl opt"><-</span> <span class="hl kwd">mean</span><span class="hl opt">(</span>school_diff$diff_6_29<span class="hl opt">)</span>
school_diff <span class="hl opt">%>%</span>
<span class="hl kwd">ggplot</span><span class="hl opt">() +</span>
<span class="hl kwd">geom_histogram</span><span class="hl opt">(</span><span class="hl kwd">aes</span><span class="hl opt">(</span>x <span class="hl opt">=</span> diff_6_29<span class="hl opt">),</span> fill <span class="hl opt">=</span> <span class="hl str">"lightblue"</span><span class="hl opt">,</span> color <span class="hl opt">=</span> <span class="hl str">"gray60"</span><span class="hl opt">,</span> alpha <span class="hl opt">=</span> <span class="hl num">0.9</span><span class="hl opt">,</span> bins <span class="hl opt">=</span> <span class="hl num">40</span><span class="hl opt">) +</span>
<span class="hl kwd">geom_vline</span><span class="hl opt">(</span>xintercept <span class="hl opt">=</span> diff_mean<span class="hl opt">,</span> linetype <span class="hl opt">=</span> <span class="hl num">2</span><span class="hl opt">) +</span>
<span class="hl kwd">annotate</span><span class="hl opt">(</span><span class="hl str">"text"</span><span class="hl opt">,</span> x<span class="hl opt">=</span>diff_mean<span class="hl opt">,</span> y<span class="hl opt">=</span><span class="hl num">10</span><span class="hl opt">,</span>
label<span class="hl opt">=</span><span class="hl kwd">fround</span><span class="hl opt">(</span>diff_mean<span class="hl opt">)</span>
<span class="hl opt">) +</span>
<span class="hl kwd">annotate</span><span class="hl opt">(</span><span class="hl str">"text"</span><span class="hl opt">,</span>
x <span class="hl opt">= -</span><span class="hl num">0.1</span><span class="hl opt">,</span> y <span class="hl opt">=</span> <span class="hl num">50</span><span class="hl opt">,</span>
label <span class="hl opt">= (</span><span class="hl kwd">fround</span><span class="hl opt">(</span><span class="hl num">1</span> <span class="hl opt">-</span> s6_greater_prob<span class="hl opt">) *</span> <span class="hl num">100</span><span class="hl opt">) %</span>cat<span class="hl opt">%</span> <span class="hl str">"%"</span><span class="hl opt">,</span>
size <span class="hl opt">=</span> <span class="hl num">8</span><span class="hl opt">,</span> color <span class="hl opt">=</span> <span class="hl str">"red"</span><span class="hl opt">,</span> alpha <span class="hl opt">=</span> <span class="hl num">0.6</span>
<span class="hl opt">) +</span>
<span class="hl kwd">annotate</span><span class="hl opt">(</span><span class="hl str">"text"</span><span class="hl opt">,</span>
x <span class="hl opt">=</span> <span class="hl num">0.36</span><span class="hl opt">,</span> y <span class="hl opt">=</span> <span class="hl num">110</span><span class="hl opt">,</span>
label <span class="hl opt">= (</span><span class="hl kwd">fround</span><span class="hl opt">(</span>s6_greater_prob<span class="hl opt">,</span> <span class="hl num">4</span><span class="hl opt">) *</span> <span class="hl num">100</span><span class="hl opt">) %</span>cat<span class="hl opt">%</span> <span class="hl str">"%"</span><span class="hl opt">,</span>
size <span class="hl opt">=</span> <span class="hl num">12</span><span class="hl opt">,</span> color <span class="hl opt">=</span> <span class="hl str">"blue"</span><span class="hl opt">,</span> alpha <span class="hl opt">=</span> <span class="hl num">0.8</span>
<span class="hl opt">) +</span>
<span class="hl kwd">annotate</span><span class="hl opt">(</span><span class="hl str">"text"</span><span class="hl opt">,</span>
label<span class="hl opt">=</span> <span class="hl str">"sd="</span> <span class="hl opt">%</span>cat<span class="hl opt">%</span> <span class="hl kwd">fround</span><span class="hl opt">(</span>diff_sd<span class="hl opt">),</span>
x<span class="hl opt">=</span><span class="hl num">0.65</span><span class="hl opt">,</span> y<span class="hl opt">=</span><span class="hl num">230</span><span class="hl opt">,</span> size<span class="hl opt">=</span><span class="hl num">10</span><span class="hl opt">,</span> color<span class="hl opt">=</span><span class="hl str">"gray30"</span>
<span class="hl opt">) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>
title <span class="hl opt">=</span> <span class="hl str">"Diferencia de Intercepción entre Escuelas 6 - 29"</span><span class="hl opt">,</span>
x <span class="hl opt">=</span> <span class="hl str">"Diferencia de Intercepción 6 - 29"</span><span class="hl opt">,</span>
y <span class="hl opt">=</span> <span class="hl str">"Cantidad de Muestras"</span>
<span class="hl opt">)</span>
<span class="hl opt">)</span>
</pre></div>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_diff_intc_6_29.png">
<img alt="JSP: Diferencia de Intercepciones de Escuelas 6 menos 29" class="pic-noborder" data-width="600" src="http://www.odelama.com/sandbox/jsp_diff_intc_6_29.png" title="JSP: Diferencia de Intercepciones de Escuelas 6 menos 29" />
</a>
<p class="caption-text"><em>JSP: Diferencia de Intercepciones de Escuelas 6 y 29.</em> En el 98% de las muestras, la Intercepción de la escuela 6 es mayor que el de la 29.</p>
</div>
<p>Como esperábamos, el $98.95\%$ de las diferencias es positivo, igual al resultado de <code>hypothesis()</code> en la columna <code>Post.Prob</code>, y su desviación estándar (sd en el gráfico) es igual a <code>Est.Error</code>. La diferencia media es $0.35$, tal como es reportado en la columna <code>Estimate</code>.</p>
<p><em>Ejemplo comparando la Intercepción por Género</em></p>
<p>En <a href="http://www.odelama.com/sandbox/#jsp_m2List">el listado de nuestro modelo <code>jsp_m2</code></a>, vemos en la intercepción por género, que las chicas <code>genderf</code> (gender:fem) tienen una muy ligera ventaja (de $0.03-(-0.01)=0.04 \, \sigma_y$) respecto a los chicos <code>genderm</code> (gender:masc). Pero no hay una gran certeza de que ninguna de estas dos intercepciones sea diferente de cero, pues sus intervalos creíbles lo incluyen. ¿Qué podemos decir acerca de la diferencia entre estas dos intercepciones? ¿Tal vez la varianza de la diferencia es tan grande que nos impide afirmar nada acerca de esta diferencia con certidumbre razonable? Resolvamos esto con la ayuda de la función <code>hypothesis()</code>.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">hypothesis</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span> <span class="hl str">"genderf - genderm > 0"</span><span class="hl opt">,</span> class<span class="hl opt">=</span><span class="hl str">"b"</span><span class="hl opt">)</span>
</pre></div>
<pre><code>Hypothesis Tests for class b:
Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (genderf-genderm) > 0 0.04 0.03 -0.01 0.09 10.33 0.91
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
</code></pre>
<p>Encontramos que en el $91\%$ de las muestras la diferencia de intercepción es a favor de las chicas.</p>
<blockquote>
<p>La alumna promedio de 3er año tiene ligeramente mejor puntaje al 5to que el alumno promedio. Los datos muestran esto con una certidumbre estimada en $91\%$, y la ventaja estimada de las chicas sobre los chicos es de $0.04 \, \sigma_y$.</p>
</blockquote>
<h3 id="efectodeclasesocialdelpadredelalumno">Efecto de Clase Social del Padre del Alumno</h3>
<p>Podemos probar como hipótesis el signo de la intercepción en cada clase <code>social</code>, es decir probar si la clase social del padre tiene un efecto positivo o negativo en el aprendizaje del alumno.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">hypothesis</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">,</span> <span class="hl str">"Intercept > 0"</span><span class="hl opt">,</span> scope<span class="hl opt">=</span><span class="hl str">"ranef"</span><span class="hl opt">,</span> group<span class="hl opt">=</span><span class="hl str">"social"</span><span class="hl opt">)</span>
</pre></div>
<pre><code>Hypothesis Tests for class :
Group Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 I (Intercept) > 0 0.05 0.07 -0.05 0.17 3.45 0.78
2 II (Intercept) > 0 0.06 0.05 -0.02 0.15 7.89 0.89
3 III nonmanual (Intercept) > 0 0.03 0.05 -0.05 0.12 2.54 0.72
4 III manual (Intercept) > 0 -0.03 0.04 -0.10 0.03 0.23 0.19
5 IV (Intercept) > 0 0.00 0.05 -0.09 0.08 0.87 0.47
6 V (Intercept) > 0 -0.07 0.06 -0.18 0.02 0.12 0.11
7 lterm unemp (Intercept) > 0 -0.06 0.06 -0.17 0.03 0.19 0.16
8 curr unemp (Intercept) > 0 -0.02 0.07 -0.14 0.08 0.62 0.38
9 father absent (Intercept) > 0 0.05 0.05 -0.02 0.14 7.53 0.88
</code></pre>
<p>Recordemos que en nuestro modelo estas intercepciones son desviaciones alrededor de la media a nivel población. Es decir, están centradas en cero y los valores positivos refieren a intercepciones arriba de la media y viceversa. Es por eso que tiene sentido evaluar la certidumbre que tenemos de estas intercepciones con respecto a cero.</p>
<p>Cuando <code>Evid.Ratio</code> es menor que 1 señala que tenemos más evidencia de que el efecto de <code>social</code> es negativo que positivo. Así, en el caso de la fila $6$, <code>Evid.Ratio=0.12</code> indica que tenemos $1/0.12=8.3$ veces más evidencia de que el efecto es negativo que positivo, y que correspondientemente la probabilidad de que sea negativa es $1-0.12=88\%$.</p>
<blockquote>
<p>La influencia de la clase <code>social</code> del padre en el rendimiento del alumno es pequeña, con una magnitud promedio de $0.04 \, \sigma_y$ y un máximo estimado de $0.07 \, \sigma_y$. Tenemos más certeza de que esta influencia es positiva para alumnos con padres en las dos primeras clases sociales y para la categoría <code>father absent</code>; y de que es negativa con padres en la última categoría <code>V</code>, o que estan desempleados a largo plazo <code>lterm unemp</code>. Sin embargo, la certidumbre que tenemos de los efectos referidos está en el orden de $[78\% - 89\%]$, hace falta un mayor volumen de datos para mejorar esta certidumbre.</p>
<p>Sorprendentemente, la categoría padre ausente <code>father absent</code> tiene un efecto positivo en el aprendizaje del alumno, con una ventaja de $0.05-(-0.07)=0.12 \, \sigma_y$ cuando es comparado al efecto de un padre en la última categoría <code>V</code>.</p>
</blockquote>
<h2 id="efectodelaescuelaenelalumno">Efecto de la Escuela en el Alumno</h2>
<p>Recordemos que con la función <code>ranef()</code> podemos obtener los coeficientes de los predictores por grupo, como la intercepción y el coeficiente de <code>score_3_z</code> por escuela.</p>
<div class="highlight-r"><pre class="hl"><span class="hl kwd">ranef</span><span class="hl opt">(</span>jsp_m2<span class="hl opt">)</span>$school <span class="hl opt">%>%</span> <span class="hl kwd">print</span><span class="hl opt">(</span>digits<span class="hl opt">=</span><span class="hl num">3</span><span class="hl opt">)</span>
</pre></div>
<p></p>
<pre><code>, , Intercept
Estimate Est.Error Q2.5 Q97.5
1 -0.326838 0.0838 -0.4910 -0.1671
2 -0.162320 0.1175 -0.4076 0.0626
3 0.102519 0.1365 -0.1658 0.3700
4 0.040418 0.0846 -0.1305 0.2011
5 -0.067128 0.0974 -0.2599 0.1283
6 -0.112183 0.1147 -0.3362 0.1088
:
, , score_3_z
Estimate Est.Error Q2.5 Q97.5
1 0.12259 0.0540 0.01752 0.23031
2 0.05322 0.0664 -0.07661 0.18670
3 -0.03970 0.0776 -0.20149 0.10940
4 -0.01508 0.0529 -0.12086 0.09166
5 0.02551 0.0611 -0.09491 0.14975
6 0.05866 0.0746 -0.07529 0.22207
:
</code></pre>
<p>Si consideramos que el 52% de los alumnos en la data son chicas, podemos estimar la intercepción esperada a nivel población con el promedio ponderado de las intercepciones por género y obtener un estimado de $0.0153$. De manera similar, si usamos la distribución de la clase social del padre en la data, con las intercepciones por este grupo, obtenemos $-0.00191$ como intercepción esperada en este concepto, lo que suma $0.0133$. Entonces la intercepción total por escuela es la suma de éste último valor con la intercepción de cada escuela, que como sabemos es el puntaje esperado del alumno promedio en 5to año. Sin embargo, al ser una constante agregada a cada intercepción de escuela, no afecta a la clasificación de estas.</p>
<div class="highlight-r"><pre class="hl">interc_school_m3 <span class="hl opt"><-</span> <span class="hl kwd">ranef</span><span class="hl opt">(</span>jsp_m3<span class="hl opt">)</span>$school <span class="hl opt">%>%</span> as<span class="hl opt">.</span>data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">()</span>
interc_school_m3$school <span class="hl opt"><-</span> row<span class="hl opt">.</span><span class="hl kwd">names</span><span class="hl opt">(</span><span class="hl kwd">ranef</span><span class="hl opt">(</span>jsp_m3<span class="hl opt">)</span>$school<span class="hl opt">)</span>
interc_school_m3 <span class="hl opt">%>%</span>
<span class="hl kwd">ggplot</span><span class="hl opt">(</span><span class="hl kwd">aes</span><span class="hl opt">(</span>x <span class="hl opt">=</span> school <span class="hl opt">%>%</span> <span class="hl kwd">fct_reorder</span><span class="hl opt">(</span>Estimate<span class="hl opt">.</span>Intercept<span class="hl opt">),</span> y <span class="hl opt">=</span> Estimate<span class="hl opt">.</span>Intercept<span class="hl opt">)) +</span>
<span class="hl kwd">geom_point</span><span class="hl opt">(</span>color <span class="hl opt">=</span> <span class="hl str">"dodgerblue"</span><span class="hl opt">) +</span>
<span class="hl kwd">geom_hline</span><span class="hl opt">(</span>yintercept <span class="hl opt">=</span> <span class="hl num">0</span><span class="hl opt">,</span> color <span class="hl opt">=</span> <span class="hl str">"gray70"</span><span class="hl opt">,</span> linetype <span class="hl opt">=</span> <span class="hl num">2</span><span class="hl opt">) +</span>
<span class="hl kwd">geom_errorbar</span><span class="hl opt">(</span><span class="hl kwd">aes</span><span class="hl opt">(</span>ymin <span class="hl opt">=</span> Q2<span class="hl opt">.</span>5<span class="hl opt">.</span>Intercept<span class="hl opt">,</span> ymax <span class="hl opt">=</span> Q97<span class="hl opt">.</span>5<span class="hl opt">.</span>Intercept<span class="hl opt">),</span> color <span class="hl opt">=</span> <span class="hl str">"darkblue"</span><span class="hl opt">) +</span>
<span class="hl kwd">coord_cartesian</span><span class="hl opt">(</span>ylim <span class="hl opt">=</span> <span class="hl kwd">c</span><span class="hl opt">(-</span><span class="hl num">0.8</span><span class="hl opt">,</span> <span class="hl num">0.8</span><span class="hl opt">)) +</span>
<span class="hl kwd">labs</span><span class="hl opt">(</span>
title <span class="hl opt">=</span> <span class="hl str">"JSP: Intercepciones por Escuela"</span><span class="hl opt">,</span> subtitle <span class="hl opt">=</span> <span class="hl str">"Intervalos al 95%"</span><span class="hl opt">,</span>
x <span class="hl opt">=</span> <span class="hl str">"Id de Escuela"</span><span class="hl opt">,</span> y <span class="hl opt">=</span> <span class="hl str">"Intercepción Estimada (95%)"</span>
<span class="hl opt">) +</span>
x_text_45
</pre></div>
<p></p>
<p><span id="plot-intc-school"></span></p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_intc_school.png">
<img alt="JSP: Intercepciones por Escuela" class="pic-noborder" data-width="850" src="http://www.odelama.com/sandbox/jsp_intc_school.png" title="JSP: Intercepciones por Escuela" />
</a>
<p class="caption-text"><em>JSP: Intercepciones por Escuela.</em> El punto representa la Intercepción estimada y las barras verticales su intervalo de variación al 95%. Las escuelas están ordenadas a lo largo del eje horizontal aproximadamente por la intercepción estimada.</p>
</div>
<p>Si con las proporciones en la data, estimamos los coeficientes para la escuela con menor intercepción, tenemos:</p>
<p>$$\begin{aligned}
\text{Escuela con menor intercepción:} \\
\text{score_5_z} &= (0.01529 - 0.00191)
+ (0.1681 + 0.00437) \, \text{raven_z} \\
&+ 0.6167 \, \text{score_3_z} +
(-0.4611 + 0.1556 \, \text{score_3_z} ) \\
&= 0.01338 + 0.1725 \, \text{raven_z} + 0.6167 \, \text{score_3_z} \\
&+ (-0.4611 + 0.1556 \, \text{score_3_z} ) \\
&= −0.4477 + 0.1725 \, \text{raven_z} + 0.7723 \, \text{score_3_z}
\end{aligned}$$</p>
<p>Donde $0.01529$ es la intercepción ponderada por <code>gender</code>, $-0.00191$ es la intercepción ponderada por <code>social</code>, $+0.00437$ es la pendiente de <code>raven_z</code> ponderada por <code>social</code>, y el último término entre paréntesis son los coeficientes de la escuela con menor intercepción; los demás son coeficientes a nivel población. Procediendo de la misma manera con la escuela con mayor intercepción:</p>
<p>$$\begin{aligned}
\text{Escuela con mayor intercepción:} \\
\text{score_5_z} &= 0.01338 + 0.1725 \, \text{raven_z} + 0.6167 \, \text{score_3_z} \\
& + (0.5606 -0.1887 \, \text{score_3_z} ) \\
&= 0.5740 + 0.1725 \, \text{raven_z} + 0.4280 \, \text{score_3_z}
\end{aligned}$$</p>
<p>Donde un <em>alumno promedio</em> en la mejor escuela tiene un puntaje $0.5740 - (-0,4477)=1.02 \, \sigma_y$ mayor que uno en la peor. Para que el alumno de la peor escuela pueda alcanzar al de la mejor, se requiere que:</p>
<p>$$\begin{aligned}
0.5740 &= −0,4477 + 0.1725 \, \text{raven_z} + 0.7723 \, \text{score_3_z} \\
1.02 &= 0.1725 \, \text{raven_z} + 0.7723 \, \text{score_3_z}
\end{aligned}$$</p>
<p>Donde una solución posible es $\text{raven_z} = 1.08$ y $\text{score_3_z}=1.08$, pues $0.1725 \times 1.08 + 0.7723 \times 1.08 = 1.02 \, \sigma_y$</p>
<blockquote>
<p>La escuela tiene un gran impacto en el aprendizaje del alumno. El efecto de la escuela con mayor influencia positiva comparada con la de mayor influencia negativa en el aprendizaje del alumno tiene una ventaja estimada en $1.02$ desviaciones estándar en el puntaje obtenido por el alumno promedio al 5to año.</p>
<p>Para que un alumno en la escuela con mayor influencia negativa pueda alcanzar al <em>alumno promedio</em> en la escuela con mayor influencia positiva, requiere que tenga puntajes en 3er año y de razonamiento abstracto con una desviación estándar por encima del promedio de la población.</p>
</blockquote>
<h2 id="clasificacindelasescuelas">Clasificación de las Escuelas</h2>
<p>Vamos a hacer una comparación cruzada de la intercepción de cada escuela contra todas las demás ($49 \times (49-1)/2=1176$ comparaciones) usando la prueba de hipótesis que vimos anteriormente. En cada comparación concluiremos que una escuela tiene mayor efecto positivo —sobre el alumno— que otra, cuando la certeza sea mayor o igual al $95\%$. El clasificar a las escuelas por intercepción, equivale a hacerlo según su efecto sobre el <em>alumno promedio</em>.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/jsp_schoolComp.png">
<img alt="JSP: Comparación de cada escuela contra todas las demás" class="pic-noborder" data-width="700" src="http://www.odelama.com/sandbox/jsp_schoolComp.png" title="JSP: Comparación de cada escuela contra todas las demás" />
</a>
<p class="caption-text"><em>JSP: Comparación de cada escuela contra todas las demás.</em> Cada punto verde señala que la escuela en la fila tiene mayor efecto (95% o más) que la escuela en la columna. Cada punto rojo señala lo opuesto, también al 95%.</p>
</div>
<p>En el gráfico arriba, cada punto de color verde señala que la escuela en el eje horizontal tiene un efecto positivo mayor sobre el alumno promedio que la escuela en el eje vertical, con una certidumbre mayor o igual que $95\%$. Por ejemplo, la escuela <code>27</code> tiene un efecto mayor que las escuelas <code>29</code> y <code>1</code>. Los puntos de color rojo señalan lo opuesto, así la <code>7</code> tiene un efecto menor que las escuelas <code>36</code>, y <code>31</code>.</p>
<p>La zona intermedia, sin puntos, contiene las combinaciones de escuelas para las que no hay certidumbre (al 95%) de cuál es mejor que la otra. Para poder resolver mejor estos casos se necesita un volumen mayor de datos. Sin embargo, hemos podido discriminar en $1042$ de $2352$ comparaciones, el $44\%$ de todas ellas, que es algo que no podríamos hacer con solo los intervalos mostrados en el gráfico <a href="http://www.odelama.com/sandbox/#plot-intc-school"><em>"JSP: Intercepciones por Escuela"</em></a> como ocurre en regresiones con estimaciones puntuales.</p>
<h1 id="referencias">Referencias</h1>
<p><span id="Burkner-1"></span>
[Burkner] P. Burkner PC (2017). “brms: An R Package for Bayesian Multilevel Models using Stan.” Journal of Statistical Software, 80(1), 1–28. doi:10.18637/jss.v080.i01.</p>
<p><span id="Gelman-1"></span>
[Gelman] Andrew Gelman, 2006. Prior distributions for variance parameters in hierarchical models. <a href="http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf">http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf</a>. </p>
<p><span id="Gelman-2"></span>
[Gelman] Andrew Gelman, Jennifer Hill, Masanao Yajima. "Why We (Usually) Don’t Have to Worry About Multiple Comparisons", 2013. Journal of Research on Educational Effectiveness. <a href="http://www.stat.columbia.edu/~gelman/research/published/multiple2f.pdf">http://www.stat.columbia.edu/~gelman/research/published/multiple2f.pdf</a></p>
<p><span id="Gelman-3"></span>
[Gelman] Andrew Gelman, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.</p>
<p><span id="Gustafson"></span>
[Gustafson] Paul Gustafson, Shahadut Hossain, and Ying C. MacNab. 2006.Conservative Prior Distributions for Variance Parameters in Hierarchical Models. The Canadian Journal of Statistics. <a href="https://doi.org/10.1002/cjs.5550340302">https://doi.org/10.1002/cjs.5550340302</a>.</p>
<p><span id="Goldstein"></span>
[Goldstein] Goldstein, H. (1995). "Multilevel Statistical Models". London, Edward Arnold.</p>
<p><span id="Joox"></span>
[Joox] Joop J. Hox. 2010. Multilevel Analysis Techniques and Applications, Second Edition. Routledge</p>
<p><span id="Lai-1"></span>
[Lai et al. 2015] Lai, Mark H. C., and Oi-man Kwok. 2015. “Examining the Rule of Thumb of Not Using Multilevel Modeling: The ‘Design Effect Smaller Than Two’ Rule.” The Journal of Experimental Education 83: 423–38. <a href="https://doi.org/10.1080/00220973.2014.907229">https://doi.org/10.1080/00220973.2014.907229</a>.</p>
<p><span id="lkj"></span>
[LKJ] Lewandowski, Daniel; Dorota Kurowicka; and Harry Joe. 2009. “Generating Random Correlation Matrices Based on Vines and Extended Onion Method.” Journal of Multivariate Analysis 100: 1989–2001.</p>
<p><span id="McElreath"></span>
[McElreath] Richard McElreath. 2018. "Statistical Rethinking: A Bayesian Course with Examples in R and Stan". CRC Press.</p>
<p><span id="Mortimer"></span>
[Mortimer] Peter Mortimer et al. "A study of effective junior schools", 1989. International Journal of Educational Research. <a href="https://www.academia.edu/24104539/A_study_of_effective_junior_schools">https://www.academia.edu/24104539/A_study_of_effective_junior_schools</a>.</p>
<p><span id="stan"></span>
[stan] Carpenter B.; Gelman A.; Hoffman M. D.; Lee D.; Goodrich B.; Betancourt M.; Brubaker M.; Guo J.; Li P.; and Riddell A. 2017. "Stan: A probabilistic programming language". Journal of Statistical Software. 76(1). 10.18637/jss.v076.i01</p>
<h1 id="anexos">Anexos</h1>
<h2 id="generacindedatossimuladostoyspanidgentoydataspan">Generación de Datos Simulados Toy <span id="genToyData"></span></h2>
<p>El código con el que se ha generado la data es:</p>
<div class="highlight-r"><pre class="hl">make_data <span class="hl opt"><-</span> <span class="hl kwa">function</span><span class="hl opt">() {</span>
beta <span class="hl opt"><-</span> <span class="hl kwd">c</span><span class="hl opt">(</span>b0 <span class="hl opt">=</span> <span class="hl num">2.7</span><span class="hl opt">,</span> b1 <span class="hl opt">=</span> <span class="hl num">1.3</span><span class="hl opt">)</span>
sigma_y <span class="hl opt"><-</span> <span class="hl num">1.6</span>
mu_alpha <span class="hl opt"><-</span> <span class="hl kwd">c</span><span class="hl opt">(</span>a0 <span class="hl opt">=</span> <span class="hl num">6.2</span><span class="hl opt">,</span> a1 <span class="hl opt">=</span> <span class="hl num">2.3</span><span class="hl opt">)</span>
sigma_alpha <span class="hl opt"><-</span> <span class="hl kwd">c</span><span class="hl opt">(</span>a0 <span class="hl opt">=</span> <span class="hl num">1.4</span><span class="hl opt">,</span> a1 <span class="hl opt">=</span> <span class="hl num">0.9</span><span class="hl opt">)</span>
rho <span class="hl opt"><-</span> <span class="hl num">0.4</span>
cov_alpha <span class="hl opt"><-</span> rho <span class="hl opt">*</span> sigma_alpha<span class="hl opt">[</span><span class="hl str">"a0"</span><span class="hl opt">] *</span> sigma_alpha<span class="hl opt">[</span><span class="hl str">"a1"</span><span class="hl opt">]</span>
J <span class="hl opt"><-</span> <span class="hl num">48</span>
mu_n_j <span class="hl opt"><-</span> <span class="hl num">20</span>
sigma_n_j <span class="hl opt"><-</span> <span class="hl num">3</span>
min_x <span class="hl opt"><-</span> <span class="hl num">0</span>
max_x <span class="hl opt"><-</span> <span class="hl num">6</span>
Sigma_alpha <span class="hl opt"><-</span> <span class="hl kwd">matrix</span><span class="hl opt">(</span><span class="hl kwd">c</span><span class="hl opt">(</span>sigma_alpha<span class="hl opt">[</span><span class="hl str">"a0"</span><span class="hl opt">]</span>^<span class="hl num">2</span><span class="hl opt">,</span> cov_alpha<span class="hl opt">,</span> cov_alpha<span class="hl opt">,</span> sigma_alpha<span class="hl opt">[</span><span class="hl str">"a1"</span><span class="hl opt">]</span>^<span class="hl num">2</span><span class="hl opt">),</span> <span class="hl num">2</span><span class="hl opt">,</span> <span class="hl num">2</span><span class="hl opt">)</span>
data <span class="hl opt"><-</span> data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>x <span class="hl opt">=</span> <span class="hl kwd">numeric</span><span class="hl opt">(),</span> y <span class="hl opt">=</span> <span class="hl kwd">numeric</span><span class="hl opt">(),</span> g <span class="hl opt">=</span> <span class="hl kwd">character</span><span class="hl opt">())</span>
coefs <span class="hl opt"><-</span> data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>intercept<span class="hl opt">=</span><span class="hl kwd">numeric</span><span class="hl opt">(),</span> slope<span class="hl opt">=</span><span class="hl kwd">numeric</span><span class="hl opt">())</span>
set<span class="hl opt">.</span><span class="hl kwd">seed</span><span class="hl opt">(</span><span class="hl num">123</span><span class="hl opt">)</span>
<span class="hl kwa">for</span> <span class="hl opt">(</span>j <span class="hl kwa">in</span> <span class="hl num">1</span><span class="hl opt">:</span>J<span class="hl opt">) {</span>
n_j <span class="hl opt"><-</span> <span class="hl kwd">sample</span><span class="hl opt">(</span><span class="hl num">2</span><span class="hl opt">:</span><span class="hl num">21</span><span class="hl opt">,</span> size<span class="hl opt">=</span><span class="hl num">1</span><span class="hl opt">,</span> replace <span class="hl opt">=</span> <span class="hl kwb">T</span><span class="hl opt">,</span> prob<span class="hl opt">=</span><span class="hl kwd">seq</span><span class="hl opt">(</span><span class="hl num">15</span><span class="hl opt">,</span><span class="hl num">10</span><span class="hl opt">,</span>length<span class="hl opt">.</span>out <span class="hl opt">=</span> <span class="hl num">20</span><span class="hl opt">))</span>
alpha_j <span class="hl opt"><-</span> MASS<span class="hl opt">::</span><span class="hl kwd">mvrnorm</span><span class="hl opt">(</span>n <span class="hl opt">=</span> <span class="hl num">1</span><span class="hl opt">,</span> mu <span class="hl opt">=</span> mu_alpha<span class="hl opt">,</span> Sigma <span class="hl opt">=</span> Sigma_alpha<span class="hl opt">)</span>
<span class="hl kwd">names</span><span class="hl opt">(</span>alpha_j<span class="hl opt">) <-</span> <span class="hl kwd">c</span><span class="hl opt">(</span><span class="hl str">"a0"</span><span class="hl opt">,</span> <span class="hl str">"a1"</span><span class="hl opt">)</span>
intercept <span class="hl opt"><-</span> beta<span class="hl opt">[</span><span class="hl str">"b0"</span><span class="hl opt">] +</span> alpha_j<span class="hl opt">[</span><span class="hl str">"a0"</span><span class="hl opt">]</span>
slope <span class="hl opt"><-</span> beta<span class="hl opt">[</span><span class="hl str">"b1"</span><span class="hl opt">] +</span> alpha_j<span class="hl opt">[</span><span class="hl str">"a1"</span><span class="hl opt">]</span>
coefs <span class="hl opt"><-</span> coefs <span class="hl opt">%>%</span> <span class="hl kwd">bind_rows</span><span class="hl opt">(</span>data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>intercept<span class="hl opt">=</span>intercept<span class="hl opt">,</span> slope<span class="hl opt">=</span>slope<span class="hl opt">,</span> g<span class="hl opt">=</span>as<span class="hl opt">.</span><span class="hl kwd">character</span><span class="hl opt">(</span>j<span class="hl opt">)))</span>
<span class="hl kwa">for</span> <span class="hl opt">(</span>i <span class="hl kwa">in</span> <span class="hl num">1</span><span class="hl opt">:</span>n_j<span class="hl opt">) {</span>
x <span class="hl opt"><-</span> <span class="hl kwd">runif</span><span class="hl opt">(</span>n<span class="hl opt">=</span><span class="hl num">1</span><span class="hl opt">,</span> min <span class="hl opt">=</span> min_x<span class="hl opt">,</span> max <span class="hl opt">=</span> max_x<span class="hl opt">)</span>
y <span class="hl opt"><-</span> <span class="hl kwd">rnorm</span><span class="hl opt">(</span>n <span class="hl opt">=</span> <span class="hl num">1</span><span class="hl opt">,</span> mean <span class="hl opt">=</span> intercept <span class="hl opt">+</span> slope <span class="hl opt">*</span> x<span class="hl opt">,</span> sd <span class="hl opt">=</span> sigma_y<span class="hl opt">)</span>
data <span class="hl opt"><-</span> data <span class="hl opt">%>%</span> <span class="hl kwd">bind_rows</span><span class="hl opt">(</span>data<span class="hl opt">.</span><span class="hl kwd">frame</span><span class="hl opt">(</span>x<span class="hl opt">=</span>x<span class="hl opt">,</span> y<span class="hl opt">=</span>y<span class="hl opt">,</span> g<span class="hl opt">=</span>as<span class="hl opt">.</span><span class="hl kwd">character</span><span class="hl opt">(</span>j<span class="hl opt">)))</span>
<span class="hl opt">}</span>
<span class="hl opt">}</span>
row<span class="hl opt">.</span><span class="hl kwd">names</span><span class="hl opt">(</span>coefs<span class="hl opt">) <-</span> <span class="hl kwb">NULL</span>
<span class="hl kwd">return</span><span class="hl opt">(</span><span class="hl kwd">list</span><span class="hl opt">(</span>data<span class="hl opt">,</span> coefs<span class="hl opt">))</span>
<span class="hl opt">}</span>
data_list <span class="hl opt"><-</span> <span class="hl kwd">make_data</span><span class="hl opt">()</span>
data <span class="hl opt"><-</span> data_list<span class="hl opt"><span class="createlink"><a href="http://www.odelama.com/ikiwiki.cgi?do=create&from=sandbox%2Ftest&page=__60__%2Fspan__62____60__span_class__61____34__hl_num__34____62__1__60__%2Fspan__62____60__span_class__61____34__hl_opt__34____62__" rel="nofollow">?</a>span><span class="hl opt"></span></span>
true_coefs <span class="hl opt"><-</span> data_list<span class="hl opt"><span class="createlink"><a href="http://www.odelama.com/ikiwiki.cgi?do=create&from=sandbox%2Ftest&page=__60__%2Fspan__62____60__span_class__61____34__hl_num__34____62__2__60__%2Fspan__62____60__span_class__61____34__hl_opt__34____62__" rel="nofollow">?</a>span><span class="hl opt"></span></span>
data <span class="hl opt"><-</span> data <span class="hl opt">%>%</span> <span class="hl kwd">add_row</span><span class="hl opt">(</span>x<span class="hl opt">=</span>data$x<span class="hl opt">[</span><span class="hl num">73</span><span class="hl opt">],</span> y<span class="hl opt">=</span>data$y<span class="hl opt">[</span><span class="hl num">73</span><span class="hl opt">],</span> g<span class="hl opt">=</span><span class="hl str">"49"</span><span class="hl opt">)</span>
data$g <span class="hl opt"><-</span> <span class="hl kwd">factor</span><span class="hl opt">(</span>data$g<span class="hl opt">,</span> levels<span class="hl opt">=</span>as<span class="hl opt">.</span><span class="hl kwd">character</span><span class="hl opt">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl num">49</span><span class="hl opt">))</span>
true_coefs <span class="hl opt"><-</span> true_coefs <span class="hl opt">%>%</span>
<span class="hl kwd">add_row</span><span class="hl opt">(</span>intercept<span class="hl opt">=</span>true_coefs$intercept<span class="hl opt">[</span><span class="hl num">6</span><span class="hl opt">],</span> slope<span class="hl opt">=</span>true_coefs$slope<span class="hl opt">[</span><span class="hl num">6</span><span class="hl opt">],</span> g<span class="hl opt">=</span><span class="hl str">"49"</span><span class="hl opt">) %>%</span>
<span class="hl kwd">mutate</span><span class="hl opt">(</span>g <span class="hl opt">=</span> <span class="hl kwd">factor</span><span class="hl opt">(</span>true_coefs$g<span class="hl opt">,</span> levels<span class="hl opt">=</span>as<span class="hl opt">.</span><span class="hl kwd">character</span><span class="hl opt">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl num">49</span><span class="hl opt">)))</span>
</pre></div>
TK Luminosity Masks Analysishttp://www.odelama.com/sandbox/tk-lumask-analysis/Oscar de Lama2020-06-08T03:52:15Z2015-04-14T20:24:29Z
<div class="options" data-disqus="yes" data-eocline="yes" data-uuid="fded5f6d-6134-4fd6-835e-78e65180a26c" style="display:none"></div>
<p>Hi Tony.</p>
<p>This page <em>is hidden</em> in the sense that do not exists any link to this page from any part of this or any other web site. The objective of this page is only to share with you some findings in my TK Luminosity mask analysis.</p>
<p>I learned about luminosity masks from <a href="http://www.amazon.com/The-Digital-Zone-System-Control/dp/1937538133">Tony Fisher's <em>"The Digital Zone System"</em></a>, from your site, and later from the material I bought from you since 2012. I recently updated the panel to the CS6 PS version.</p>
<p>I have done this analysis —as usual in my case— mainly for pure curiosity and also looking for actionable knowledge that can help me to develop better photos. You can find this kind of analysis about different photography subjects in my site. I am a software engineer actually studying data science (which is very useful for digital photograph analysis). Some of <a href="http://photos.odelama.com/">my photos are here</a> and <a href="https://www.flickr.com/photos/roofwalker/">here</a>.</p>
<p>I always avoid as possible the use of technical jargon and terms, to make this content readable for any person with just some basic knowledge about digital images. If the text starts to seem too technical, please bear with me and keep reading, because probably the matter will be explained in better way few lines later.</p>
<p>This analysis considers the three most common color spaces: <em>Adobe RGB</em>, <em>ProPhoto</em> and <em>sRGB</em>. They are also the ones suggested in the TK package documents.</p>
<p>Finally, my English is not yet as good as I would like. Please feel free to ask me better explanations about the topics you want.</p>
<p><strong>License</strong></p>
<p>By sharing this with you, I will not expect, and much less demand, any retribution from you; I just want a better tool. If you enhance the TK tool somehow based on this material, it would be nice. Otherwise, I will just implement and try ideas from this analysis by myself.</p>
<p>I am sharing this page with you under a <a href="https://creativecommons.org/licenses/by/3.0/">CC Attribution license</a>: You can freely include, in your web page or product, content from this page: If you use a graph, code, ideas or anything —directly or derived — from this page, I will appreciate the attribution or mention of my name (and/or my site).</p>
<div class="ContentTable">
<header class="tocHeader">Contents</header>
</div>
<h1 id="tonalvalues">Tonal Values</h1>
<p>We have made the analysis building luminosity masks for a gradient of neutral colors. A neutral color is encoded having the same pixel value on each RGB channel (e.g. as in <code>(127, 127, 127)</code>). In turn, those pixel values are encoded in each <em>color space</em> (<em>color profile</em> in <em>Photoshop</em> terms) using different tonal scales. This way, two different color spaces may have similar, or even exactly the same tonal scale, and however depict the colors with different RGB pixel values. As this may seem confusing; we will make a short review of this subject.</p>
<p>The tonal scales encode gray colors, from black (<code>0%</code>) to white (<code>100%</code>). The root of all the tonal scales is the <em>linear</em> tonal scale, which is the scale used in the camera raw image data, and has the dark tones very compressed in the lower values (to the left side of the scale). This is so notorious that the perceptual middle gray; this is the gray that seems halfway between black and white; is just the <code>18.4% linear</code>, very low for a "middle gray" tone, which should have a value closer to <code>50%</code> as the <em>middle</em> term suggests.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/tones.png">
<img alt="Tonal Scales" class="pic-noborder" data-width="766" src="http://www.odelama.com/sandbox/tones.png" title="Tonal Scales" />
</a>
<p class="caption-text"><strong>(Tonal Scales)</strong> The white ticks mark the same perceptual middle gray.</p>
</div>
<p>These clumped-to-the-left linear dark tones, cause technical issues that are addressed by the different <em>non-linear</em> tonal scales: They <em>map</em> the perceptual tones referred by the linear tone values to a different tonal scale, using a function which is basically a <em>gamma curve</em>, that look like this:</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/tonal-conversion.png">
<img alt="Tonal Scale from linear" class="pic-noborder" data-width="365" src="http://www.odelama.com/sandbox/tonal-conversion.png" title="Tonal Scale from linear" />
</a>
<p class="caption-text">Close to the origin the curve may have a linear segment, represented by the blue line.</p>
</div>
<p>For example, the <em>Adobe RGB</em> and <em>sRGB</em> color spaces, use a curve based on a <code>2.2</code> gamma; the <em>ProPhoto</em> space uses <code>1.8</code> gamma, and the <em>L</em> one (The L —Lightness— component in the <em>Lab</em> color space) uses a <code>3</code> gamma. The greater the gamma value, the higher is pulled the curve to above.</p>
<p>In general, this curve produces tonal scales with more room for the dark tones —as desired— <em>taking the space</em> from the lighter ones: the curve stretches the linear dark area while compressing the lighter one. The higher the gamma value, the more are stretched the dark tones to the right. We can notice this in the <em>tonal scales</em> graph above, where the middle gray is moved from the far left to the middle of the scale in the <em>L</em> tonal scale.</p>
<p>As we know, each RGB channel can be seen as a gray scale, and correspondingly, a color space uses these non-linear tonal scales to encode the tonal values in their RGB channels, representing from <code>0%</code> to <code>100%</code> the channel color intensity. However, they do a step ahead, mapping each channel <code>100%</code> value with a specific hue of the RGB colors, which makes them encompass different color spaces. Nonetheless —internally— each RGB channel intensity value is mapped from linear to non-linear values using the tonal scales we have seen above.</p>
<h1 id="introductiontothekindsofgraphswewilluse">Introduction to the Kinds of Graphs we Will Use</h1>
<p>We will use graphs plotting the source tonal values (image pixel values) in the <code>x</code> axis and the luminosity-mask transparency in the <code>y</code> axis.</p>
<p>In the <code>x</code> axis, <code>0%</code> refers to the minimum and <code>100%</code> refers to the maximum pixel value intensity. Of course, this is relative to the channel we are talking about; in the red channel a <code>100%</code> tonal value means the maximum intensity of red.</p>
<p>In this analysis I am considering only neutral colors. For this reason, a <code>0%</code> value in <code>x</code> corresponds to black, and <code>100%</code> to white. These tonal values are expressed in the tonal scale referred by the color space in use (e.g. <em>Adobe RGB</em> tonal scale).</p>
<p>I am using the percentage format for tonal values, because we are dealing with 16-bit images, whose pixel values are far beyond of <code>[0, 255]</code>. This way —for example— <code>50% ProPhoto</code> is easier to grasp than <code>8191 ProPhoto</code>.</p>
<p>Higher values in the <code>y</code> axis mean the luminosity mask is exposing more the source pixel: a <code>100%</code> value corresponds to <em>total revealing</em> and <code>0%</code> means <em>total concealing</em>.</p>
<p>For example, the following graph shows the Lights masks in <em>Adobe RGB (aRGB)</em>.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/adobeRGB-Light-FFF.png" id="aRGB-Light-FFF">
<img alt="Adobe RGB Light Masks" class="pic-noborder" data-width="598" src="http://www.odelama.com/sandbox/adobeRGB-Light-FFF.png" title="Adobe RGB Light Masks" />
</a>
<p class="caption-text">Adobe RGB Lights Masks</p>
</div>
<p>You can make a reality check of this values in Photoshop <em>(PS)</em> using the <code>Information</code> panel showing the RGB values in 32 bits (floating-point) while checking the channel containing the luminosity masks.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/check-values.png" id="aRGB-Light-FFF">
<img alt="Checking Values" class="pic-noborder" data-width="544" src="http://www.odelama.com/sandbox/check-values.png" title="Checking Values" />
</a>
<p class="caption-text">Checking the highest mask value of (6) Mid-tone Lights: 0.463 aRGB</p>
</div>
<p>In the example above, the highest <code>(6) Mid-tone Lights</code> mask value is <code>46.3%</code> which pretty much corresponds with the value shown in the plot.</p>
<h2 id="graphshowingthetonevaluesinllab">Graph Showing the Tone Values in L (Lab)</h2>
<p>A given tonal value, in different tonal scales, may correspond to different perceptual tones. For example <code>18% Adobe RGB</code> (which we will denote using <code>18% aRGB</code>) does not correspond to the same tone we see (hence "perceptual") with <code>18% proPhoto RGB</code>.</p>
<p>Considering this, is useful to show the plot with the <code>x</code> tonal values marks expressed in <em>L</em> (the <em>L</em> <em>"Lightness"</em> component in the <em>Lab</em> space).</p>
<p>Notice in the following plot the <code>x</code> axis is still labeled as <code>Adobe RGB Tonal Scale</code> but it shows mark values converted to <em>L</em>.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/adobeRGB-Light-TFF.png" id="aRGB-Light-TFF">
<img alt="Adobe RGB Light Masks" class="pic-noborder" data-width="598" src="http://www.odelama.com/sandbox/adobeRGB-Light-TFF.png" title="Adobe RGB Light Masks" />
</a>
<p class="caption-text">Adobe RGB Light Masks. The marks are showing equivalent L(ab) values. However, as the x label says they are spaced as Adobe RGB</p>
</div>
<p>Notice this is not —at all— the same as converting all the source tonal values from <em>Adobe RGB</em> to <em>L</em>. This is exactly the same plot shown in the previous section. However, the <code>x</code> tick labels show the <em>L</em> equivalent value. For example, we can see, the second mark for <code>10% aRGB</code> corresponds to <code>9.0 L</code></p>
<p>This is similar to the case when —for example— the <code>x</code> marks are logarithmically spaced, as in photo stops for example. In this example, the <code>x</code> marks on the plot are really located and spaced at positions proportional to <code>0, 1, 2, 3 ...</code>. However, the plot may show under those marks <code>1, 2, 4, 8 ...</code> as labels: The <code>x</code> axis is logarithmically spaced, but shows the <em>natural</em> <code>x</code> values (instead their logarithms) under each mark.</p>
<p>In the same sense, in the plot above, the <code>x</code> axis is drawn and spaced as <em>Adobe RGB</em> tonal values. Nonetheless, is showing the equivalent <em>L</em> tonal values.</p>
<h2 id="graphplottedwiththetonevaluesinllab">Graph Plotted with the Tone Values in L (Lab)</h2>
<p>Sometimes we will show a plot with the <code>x</code> tonal values converted to the <em>L</em> (Lab) tonal scale. In this case, the <code>x</code> axis is labeled as <code>L (Lab) Tonal Scale</code>. Because of this transformation, some parts of the non-converted plot are shown horizontally stretched or compressed in the converted plot, horizontally "deforming" the plot without such transformation.</p>
<p>The following plot is the same shown in previous sections, but with the <code>x</code> axis tonal values converted to <em>L</em>. As the <em>Adobe RGB</em> tonal scale is similar to the <em>L</em> one, except for the darker values, the shape of this image is similar, but not exactly equal to the previous ones; and is more deformed close to the axes origin.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/adobeRGB-Light-TTF.png" id="aRGB-Light-TFF">
<img alt="Adobe RGB Light Masks" class="pic-noborder" data-width="598" src="http://www.odelama.com/sandbox/adobeRGB-Light-TTF.png" title="Adobe RGB Light Masks" />
</a>
<p class="caption-text">Adobe RGB Light Masks. The <code>x</code> axis is show with the tonal values in L.</p>
</div>
<p><span id="tonal-scales-sum"></span></p>
<h3 id="whywewantplotswiththellabtonalequivalents">Why we Want Plots with the L (Lab) Tonal Equivalents</h3>
<p>The <em>Lab</em> color space was built with the intention to get color values correlated with the perceptual color (the color we see). This means, the more different two colors look, the greater will be the difference between their <em>Lab</em> color values. This desired property is something that does not occur with most color spaces, at least not along all their possible values.</p>
<p>For the case of neutral colors, this means the <em>L</em> tonal scale describes better how different or similar will look two tonal values. Additionally, the ubiquitous <code>18% linear</code> photographic middle-gray corresponds exactly to <code>50L</code>, the middle of the <em>L</em> scale. As consequence, this kind of <em>L</em> plot, describes better which perceptual tones (lights, mid-tones, darks) will be really affected by a given luminosity mask, regardless (or despite of) the tonal scale related to the color space we are using.</p>
<p>As we have already seen, the <em>Adobe RGB</em> tone values correspond pretty much to the <em>L</em> ones, and consequently the masks will affect the targeted perceptual tones.</p>
<p>This match between <em>Adobe RGB</em> and <em>L</em> occurs because the <em>L</em> tonal scale is based on a gamma of <code>3</code>, while <em>Adobe RGB</em> and <em>sRGB</em> are both based on a gamma of <code>2.2</code>. Even when <code>3</code> doesn't seems very close to <code>2.2</code>, they are close enough to make both, the <em>sRGB</em> and the <em>Adobe RGB</em> tonal scale values close to the <em>L</em> ones, as we will see in the following plot.</p>
<p>On one hand, the <em>Adobe RGB</em> and <em>sRGB</em> tonal scales are very similar to each other (because of the fact that both are based on a <code>2.2</code> gamma) except for the darkest values. On the other hand, the <em>ProPhoto</em> color space uses tonal values based on a <code>1.8</code> gamma, which is not close at all to the <em>L</em> <code>3</code> gamma. This will cause some issues in the luminosity masks that we will address later.</p>
<p>The following plot shows how the <em>Adobe RGB</em> and <em>sRGB</em> tonal scales are pretty close to the <em>L</em> one (represented by the dark gray dashed diagonal line) and how they are close to each other except for the darkest tones. However, the <em>ProPhoto</em> tonal curve is notoriously different to all the other ones.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/sRGB-vsAdobeRGB.png">
<img alt="Adobe RGB vs sRGB" class="pic-noborder" data-width="441" src="http://www.odelama.com/sandbox/sRGB-vsAdobeRGB.png" title="Adobe RGB vs sRGB" />
</a>
<p class="caption-text">Tonal scales versus the L one. The L tonal scale is represented by the dark gray dashed diagonal line.</p>
</div>
<p>Notice how for the darkest tones, the <em>sRGB</em> tonal values are closer to the <em>L</em> ones than those in the other curves. This nice property, and the fact that the <em>ProPhoto</em> tonal scale is too far from <em>L</em>, is the reason why <em>Adobe Lightroom</em> uses internally the <em>ProPhoto</em> color space but with the <em>sRGB</em> tonal scale (instead the <code>1.8</code> gamma tonal scale prescribed by the <em>ProPhoto</em> definition), which sometimes is informally called the <em>Melissa</em> color space.</p>
<h2 id="normalizingthemaskscalescurveheight">Normalizing the Mask Scales Curve Height</h2>
<p>The beauty of each luminosity masks is in how they target the pixels within some luminosity range and very softly diminishes its transparency to protect the pixels with luminosity different to the targeted one. This relative mask values —masking profile— is the more relevant characteristic of the TK masks.</p>
<p>On the contrary, the absolute mask value for each pixel luminosity is not so relevant. We can say that because if the mask exposes too much the targeted image pixels, the image developer can use softer adjustment intensity or diminish the adjustment layer opacity; and if the mask exposes the targeted pixels too little, greater adjustment intensity or the cloning of that adjustment solves the issue.</p>
<p>All of this means that if two given masks have the same profile, but just different scales, one of those masks is redundant; because the same adjustment effect can be produced with the other mask by just varying its intensity as we described above.</p>
<p>To fairly compare the profile of two masks, it is required to normalize their scale so they both have the same maximum transparency value. For example, the following two masks may seem different:</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/same-profile.png">
<img alt="Adobe RGB Mid-tone masks" class="pic-noborder" data-width="451" src="http://www.odelama.com/sandbox/same-profile.png" title="Adobe RGB Mid-tone masks" />
</a>
<p class="caption-text">Both mask curves have their natural scale.</p>
</div>
<p>But if we scale them up, normalizing their heights, we will find something different:</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/same-profile-normal-height.png">
<img alt="Adobe RGB Mid-tone masks" class="pic-noborder" data-width="451" src="http://www.odelama.com/sandbox/same-profile-normal-height.png" title="Adobe RGB Mid-tone masks" />
</a>
<p class="caption-text">Both mask curves have their height scaled to match at 100%.</p>
</div>
<p>We can see they have almost the very same profile. So, they are equivalent and consequently one is redundant!</p>
<p>The subtitle <code>With normalize y Scale</code> will be an indication that the mask curves have a normalized height.</p>
<h2 id="howdidwegetthedatafortheseplots">How did we get the Data for these Plots?</h2>
<p>Except when indicated otherwise, these plots were built from images prepared with <em>PS</em> using the settings in the recommendations given in the "AA--Setting Up the Color Working Space.pdf" document, which is included in the <em>TKActions</em> package.</p>
<p>All the tests were done using <em>all the time</em> 16-bit color depth. All the tests for each color space was done in a single <code>.psd</code> document, with the <em>PS</em> color settings adjusted before the test and never changed again for the resulting <code>.psd</code> working document.</p>
<p>The basic test procedure is:</p>
<ol>
<li>Create a new document with 512 bits of width in the <code>Lab</code> color mode (16-bits).</li>
<li>Paint a horizontal gradient from black to white and from end to end, in the luminosity channel.</li>
<li>Convert the document (menu option <em>"/Image/Mode/RGB color"</em>) to the current RGB working space.</li>
<li>This <em>source</em> document is exported as a (16-bit) <code>.tif</code> file. Notice the three channels will contain the same information.</li>
<li>Build a luminosity masks as a channel using the TK panel.</li>
<li>Create new image layer and copy the luminosity channels into the image layer RGB channels (e.g. the red channel).</li>
<li>The image created in the previous step is exported (the <em>masks</em>) as a (16-bit) <code>.tif</code> file.</li>
<li>The pixel <em>source</em> values are plotted in <code>x</code> axis and the <em>mask</em> value corresponding to each <em>source</em> value is plotted in the <code>y</code> axis.</li>
</ol>
<p>I am including a section with all the R-Language code used for this analysis. Probably you don't know how to use it, but you can share it with a friend who can handle it, either for validation or for further analysis. If you want I can help you with further analysis.</p>
<h1 id="adobergbluminositymask">Adobe RGB Luminosity mask</h1>
<p>In this analysis we will make remarks and suggestions valid for the masks in all the other color spaces.</p>
<h2 id="adobergblightsmasks">Adobe RGB Lights Masks</h2>
<p>As prescribed, the <em>PS</em> color settings are <code>RGB: Adobe RGB</code> and <code>Gray: Gray Gamma 2.2</code>.</p>
<p>In the following gadget, please click on the titles at the right side to see the graph corresponding to each of them.</p>
<div class="slidorion" id="aRGB-Light_slider">
<div class="slider" data-height="477" data-width="500">
<div class="slide"><img data-src="aRGB-Light-FFF.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Light-FFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Light-TFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Light-TTT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The Mid-tone and the 1/4 Lights are normalized to the same height.</div>
<div class="header">Normalized and showing L values</div><div class="content">Showing original tones labeled with equivalent L values.</div>
<div class="header">Normalized with Tones in L</div><div class="content">With the tonal values converted to L.</div>
</div>
<p></p></div>
<p>It is very nice to see how each mask slowly conceals more the pixels with their luminosity out of the mask target.</p>
<p>Surprisingly, the <em>"Lights"</em> mask is just a diagonal line, which means the mask strength is linear with respect to each pixel luminosity.</p>
<p>The <em>"Lights Lights"</em> mask corresponds to a curve with <code>1/2</code> gamma and the following <em>"Bright Lights"</em> mask corresponds to one with <code>1/4</code> gamma. These approximations were found accurately (not eye balled for example). However, I won't bore you with the technical stuff.</p>
<p>The gamma pattern in these curves, means that they have lower and lower values from light to dark tones, but only really reach to zero in the black tone <code>0%</code>. The lower the gamma, the steep is the dive to the <code>0%</code> transparency level, but the curve hits that level only at the <code>0%</code> black tone: the axes origin. This is the mathematical explanation of how they gradually protect the darkest tones in these Lights masks.</p>
<p>In the <em>PS</em> reality, this is limited by the image bit depth. So, in 16-bit images the channels and hence the mask values are saved as 16-bit integers, which means the <em>real</em> channel values are quantized to <code>1/(2^16-1) = 1/32767 = 0.003%</code> steps, this means values below this limit are considered zero, which of course is still a great precision.</p>
<div class="leftIndent">
<h3 id="aboutthemid-tonelightsand14lightsmasksheight">About the "Mid-tone Lights" and "1/4 Lights" Masks Height</h3>
<p>It would be great to have <em>"Mid-tone Lights"</em> and <em>"1/4 Lights"</em> masks with transparency values above <code>50%</code>. As we have commented before, a change of mask scale won't change its profile and consequently won't damage its overall adjustment effect; and this will bring the benefit to see <em>the PS marching ants</em> when using those selections. I have normalized them to <code>61.8% aRGB</code> —the <em>golden ratio</em>— just for not to pick a totally random value somewhat above <code>50%</code>.</p>
<p>A greater transparency scale will be also beneficial, because, even if it is a little greater than you need, it is easy to soften by just lowering the adjustment layer opacity. Also, avoiding <em>too weak</em> masks is beneficial because this way we need less often to clone an adjustment layer to raise its strength, which is rather annoying and hits the document size.</p>
<p>However, I don't think the aforementioned masks are <em>too weak</em>. The above comments are included for the sake of completeness about the benefits of having not <em>too weak</em> masks.</p>
<p>The scale of the masks can be changed in PS converting the mask to quick mask mode and applying an image curve adjustment, where the curve is a straight line. For example, to raise the scale in a factor of <code>1.335</code> the adjustment line would start at <code>(0,0)</code> and finish at <code>(255/1.335=191, 255)</code>. After the scale adjustment, the mask is reverted to normal mask mode, and the next steps can be resumed as normal.</p>
<p>As software engineering myself, I understand that changing too much the product from one version to another, might upset some old users. However, in this case, both aforementioned masks have their peak almost to the same height and not too far to the <code>50%</code> limit. So, I would consider in first place to set them to the same height and in second place raise them a little above the <code>50%</code> limit.</p>
</div>
<p>The <em>"Mid-tone Lights"</em> mask is intended to be focused to pixels with <code>160/255 = 63%</code> luminosity. The observed mask peak is around <code>68% aRGB</code> which is a negligible difference.</p>
<p>The <em>"1/4 lights"</em> mask is intended to be focused in pixels with <code>192/255 = 75%</code> of luminosity. However in <em>Adobe RGB</em> and in <em>L</em> tonal values its peak is in <code>83% aRGB</code> and <code>85 L</code> respectively. Not very far, but nevertheless too far from the <code>75%</code> target. To solve this issue I have as suggestion the mask shown in the following graphs with the name "(8+) 1/4 Mid-Lights", which has a peak exactly at <code>75% aRGB</code>.</p>
<div class="slidorion" id="aRGB-Light-Sug_slider">
<div class="slider" data-height="477" data-width="500">
<div class="slide"><img data-src="AdobeRGB-Lights-mid-suggest-nat.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="AdobeRGB-Lights-mid-suggestion.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The mask curves are normalized to have the same height.</div>
</div>
<p></p></div>
<p>The <em>"(6) Mid-tone Lights (ML)"</em>, <em>"(7) 1/4 Lights (1/4)"</em> and <em>"(8+) 1/4 Mid-Lights"</em> mask curves are normalized to the same height using as corresponding factor <code>1.335573</code>, <code>1.459237</code> and <code>1.3586</code>.</p>
<p>The proposed <em>"(8+) 1/4 Mid-Lights"</em> mask was built with the following formula, where the <em>"Super Lights"</em> mask is subtracted twice from the "<em>Lights"</em> mask.</p>
<pre><code>1/4 Mid-Lights = Lights - Super Lights - Ultra Lights - Darks - Super Lights
</code></pre>
<p>The proposed curve profile is very similar to the "(6) Mid-tone Lights (ML)" curve but with its peak shifted over the <code>75% aRGB</code> mark; those masks protect the <em>"Super Lights"</em> far more than the <em>"1/4 lights"</em> mask does.</p>
<p>The proposed mask is a good addition, if not a replacement, to the original <em>"1/4 lights"</em> mask. However, I haven't test yet this mask in my actual photo development.</p>
<h2 id="adobergbmid-tonemasks">Adobe RGB Mid-Tone Masks</h2>
<p>There are three mid-tones masks that when their height is normalized, they look pretty much the same.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/AdobeRGB-Mids-redundant.png">
<img alt="Redundant Mid-Tone masks." class="pic-noborder" data-width="500" src="http://www.odelama.com/sandbox/AdobeRGB-Mids-redundant.png" title="Redundant Mid-Tone masks." />
</a>
<p class="caption-text">Redundant Mid-Tone masks.</p>
</div>
<p>These masks are not exactly equal to each other. However, they are too much alike to have a mask devoted to each curve.</p>
<p>The following graphs show the actual set of mid-tone masks curves.</p>
<div class="slidorion" id="aRGB-Mid_slider">
<div class="slider" data-height="477" data-width="500">
<div class="slide"><img data-src="aRGB-Mid-FFF.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Mid-FFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Mid-TFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Mid-TTT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The masks curves are normalized to have the same height.</div>
<div class="header">Normalized and showing L values</div><div class="content">Showing original tones labeled with equivalent L values.</div>
<div class="header">Normalized with Tones in L</div><div class="content">With the tonal values converted to L. The curves are slightly tilted to the light side.</div>
</div>
<p></p></div>
<p>When plotted with the tonal values converted to <em>L</em>, the curves are slighted tilted to the light side, but this tilt is negligible.</p>
<p>To handle the issue of having three very similar mid-tones masks, I suggest one that I thought I will find among those mid-tone curves: The sum of the <em>"Mid-tone Lights"</em> and <em>"Mid-tone Darks"</em> which is denoted as <em>"(7+) Mid Mid-tones"</em> in the following graphs.</p>
<pre><code>Mid Mid-tones = Mid-tone Lights + Mid-tone Darks
</code></pre>
<p>To target the mid-tones in a tighter way, I have another suggestion which is denoted as "(6+) Inner Mid-tones". This curve is obtained like the previous one , but also subtracting the <em>"Lights Lights"</em> and <em>"Darks Darks"</em> masks:</p>
<pre><code>Inner Mid-tones = Mid Mid-tones - Lights Lights - Darks Darks
</code></pre>
<p>From the three similar mid-tones masks, in the following height-normalized graph, I am leaving the central one <em>"Basic Mid-tones"</em> and including the above suggested mid-tones masks.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/AdobeRGB-Mids-suggestions.png">
<img alt="Mid-Tone masks." class="pic-noborder" data-width="520" src="http://www.odelama.com/sandbox/AdobeRGB-Mids-suggestions.png" title="Mid-Tone masks." />
</a>
<p class="caption-text">Mid-Tone masks with two added suggestions: (6+) and (7+).</p>
</div>
<p>Notice how the additional curves are gradually more bell shaped, protecting more the lightest and darkest tones, as expected from a mid-tone curve.</p>
<div class="leftIndent">
<h3 id="aboutthemid-tonescurvesheight">About the Mid-Tones Curves Height</h3>
<p>I really think the <em>"Basic Mid-tones"</em> mask, as-is (natural), is <em>too weak</em>. Also, I think that realizing that three of these mid-tone curves have almost the same profile is an unavoidable call to enhance them. I would take advantage of this situation to raise the height of the <em>weak</em> <em>"Basic Mid-tones"</em> mask curve above the <code>50%</code> limit and I would also introduce the new masks matching that height, leaving the <em>"Super Mid-tones"</em> and <em>"Ultra mid-tones"</em> masks at their current similar level.</p>
<p>This way, the suggested set of mid-tone masks would have the following curves:</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/aRGB-Mid-Enhanced.png">
<img alt="Mid-Tone masks." class="pic-noborder" data-width="531" src="http://www.odelama.com/sandbox/aRGB-Mid-Enhanced.png" title="Mid-Tone masks." />
</a>
<p class="caption-text">Suggested Mid-Tone curves.</p>
</div>
<p>In the above graph the <em>"(1) Basic Mid-tones"</em>, <em>"(6+) Inner Mid-tones"</em> and <em>"(7+) Mid Mid-tones"</em> have been normalized to the same height using as corresponding factor <code>2.462463</code>, <code>1.876162</code> and <code>0.8208377</code>.</p>
</div>
<h2 id="adobergbdarksmasks">Adobe RGB Darks Masks</h2>
<p>The Darks masks curves are equivalent to the Lights masks but flipped horizontally. All the suggestions and remarks to the Lights curves are valid for these Darks curves.</p>
<h2 id="adobergbzonemasks">Adobe RGB Zone Masks</h2>
<p>The following graphs show the actual set of Zone masks curves.</p>
<div class="slidorion" id="aRGB-Mid_slider">
<div class="slider" data-height="492" data-width="575">
<div class="slide"><img data-src="aRGB-Zone-FFF.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Zone-FFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Zone-TFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="aRGB-Zone-TTT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The masks curves are normalized to have the same height.</div>
<div class="header">Normalized and showing L values</div><div class="content">Showing original tones labeled with equivalent L values.</div>
<div class="header">Normalized with Tones in L</div><div class="content">With the tonal values converted to L. The curves are slightly tilted to the light side.</div>
</div>
<p></p></div>
<p>I am still a rookie using the Zone system. However I understand very well the schema and found it very useful in my photo development. Nonetheless, I don't look for selecting image pixels in a zone; I work with the regular luminosity masks and every now and then I check if some image areas are in the intended zone.</p>
<p>In general, the zones 2 to 6 are the most important because they will contain most of the image detail, and in particular the middle zones, 3 to 5, are still more important, because there will be the most clearly visible detail. Of course this may change with the user aesthetic intent, but is valid in the general case targeted by the zone system. Therefore, these zones must be crafted with the best quality.</p>
<p>Having said that, I see the zone 3 to 5 masks calling for enhancement for the following reasons:</p>
<ol>
<li>They don't protect enough the darkest and lightest tones.</li>
<li>The zones 3 and 5 are not centered in the expected luminosity zones.</li>
</ol>
<p>There is a too big change of protection of the darkest tones between the zone 5 and 6, mainly because the zone 5 doesn't protect enough the dark tones. Symmetrically, the same occurs between the zones 3 and 4 with respect to the lightest tones. This is flagged in the following image with (1).</p>
<p>In similar sense, the zone 5 should protect more those darkest and lightest tones with a shape more bell alike. Perhaps for this zone is more useful the proposed mid-tones mask <em>"(7+) Mid Mid-tones"</em>; this is flagged with (2).</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/aRGB-Zone-suggestions.png">
<img alt="Zone masks." class="pic-noborder" data-width="575" src="http://www.odelama.com/sandbox/aRGB-Zone-suggestions.png" title="Zone masks." />
</a>
<p class="caption-text">Areas with enhancement potential.</p>
</div>
<p>Another aspect that would enhance the zone masks is making them more targeted to the <em>L</em> tonal values. Even when there are proponents of different sets of specific RGB values defining the zones, the most commonly used, is the one with the zones centered in the multiples of <code>10L</code> (e.g. this is proposed by Lee Varis in his book <a href="https://books.google.com.pe/books?id=GVALAAAAQBAJ&printsec=frontcover&dq=lee+varis+mastering+exposure+and+the+zone+system&hl=en&sa=X&ei=gtsvVYXDC4aiNoSFgdAG&redir_esc=y#v=onepage&q=lee%20varis%20mastering%20exposure%20and%20the%20zone%20system&f=false">"mastering exposure and the zone system"</a>. This way, zone 1 is centered in <code>10L</code>, the zone 2 in <code>20L</code>, etc. Each zone covers ±<code>5L</code>. So, at the end the zone 0 is in <code>[0, 5]L</code>, zone 1 in <code>[5, 15]L</code> and so on, at least approximately.</p>
<p>This makes a lot of sense, because being based in <em>L</em>, the zones are evenly spaced in the perceived luminosity, which is the gist of the zones concept. This is also great to get consistent results regardless which specific color space you are using in the photo development.</p>
<p>From this perspective, there is a huge gap between the zones 3 and 4: The zone 3 is centered around <code>25L</code> (which is acceptable close to the desired <code>30L</code>), but the zone 4 is centered in <code>45L</code>: this is a two zones difference! Of course the same happens with the zones 6 and 7. Also from this perspective, the zones 1 and 2 are too dark and the zones 8 and 9 too light.</p>
<h2 id="adobergbmasks:allcurvesanddata">Adobe RGB Masks: All Curves and Data</h2>
<p>With the following gadget you can explore, zoom-in (double click), zoom-out and even download the data behind the plots. When you hover the mouse over the curves, a tool-tip window will show you curve attributes related to the mouse position.</p>
<p>You can select:</p>
<ol>
<li>The <em>type of mask</em> you want to see <em>(Darks, Lights, Mid-tones, Zone)</em>.</li>
<li>If the tonal scale in the <code>x</code> axis is in the <em>Adobe RGB</em> tonal scale (<em>false</em>) or in <em>L</em> (<em>true</em>).</li>
<li>If the curves must have their height normalized.</li>
<li>Which individual masks you want in the plot.</li>
</ol>
<p>The different masks are encoded by color. At the bottom right corner there is the mask color legend.</p>
<div class="Tableau-content">
<script src="https://public.tableau.com/javascripts/api/viz_v1.js" type="text/javascript"></script>
<div class="tableauPlaceholder" style="width: 804px; height: 669px;"><noscript><a href="http://www.odelama.com/sandbox/#"><img alt="Dashboard 1 " src="https://public.tableau.com/static/images/lu/lumaskAdobeRGB-AllMasks/Dashboard1/1_rss.png" style="border: none" /></a></noscript><object class="tableauViz" height="669" style="display:none;" width="804"><param name="host_url" value="https%3A%2F%2Fpublic.tableau.com%2F" /> <param name="site_root" value="" /><param name="name" value="lumaskAdobeRGB-AllMasks/Dashboard1" /><param name="tabs" value="no" /><param name="toolbar" value="yes" /><param name="static_image" value="https://public.tableau.com/static/images/lu/lumaskAdobeRGB-AllMasks/Dashboard1/1.png" /> <param name="animate_transition" value="yes" /><param name="display_static_image" value="yes" /><param name="display_spinner" value="yes" /><param name="display_overlay" value="yes" /><param name="display_count" value="yes" /><param name="showVizHome" value="no" /><param name="showTabs" value="y" /><param name="bootstrapWhenNotified" value="true" /></object></div>
<p></p></div>
<h1 id="prophotorgbluminositymasks">ProPhoto RGB Luminosity Masks</h1>
<p>In the section <a href="http://www.odelama.com/tonal-scales-sum"><em>Why we Want Plots with the L (Lab) Tonal Equivalents</em></a> we saw the <em>PropPhoto RGB</em> color space (<em>pRGB</em>) uses a tonal scale based in <code>1.8</code> gamma and that makes its tones a little too far from the <em>L</em> scale; and that brings some issues when working with luminosity masks.</p>
<p>When the luminosity masks are plotted against the <em>L</em> tonal scale, it is evident their target is shifted to lighter tones:</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/prophoto-right-shift.png">
<img alt="ProPhoto Right Shift" class="pic-noborder" data-width="800" src="http://www.odelama.com/sandbox/prophoto-right-shift.png" title="ProPhoto Right Shift" />
</a>
<p class="caption-text">The luminosity target of each mask is shifted to lighter tones.</p>
</div>
<p>Notice how the peak of the Mid-tones masks is at <code>60L</code>. The same happens with the zones masks, the zones 4, 5, and 6 have peaks at <code>52L</code>, <code>60L</code> and <code>65</code>, which is not acceptable for the zone system.</p>
<p>A solution to deal with this is use of the <em>Melissa RGB</em> color space; which is the ProPhoto gamut of colors but encoded with the sRGB tonal scale. I have found <a href="https://sites.google.com/site/chromasoft/icmprofiles">the ICC/ICM profile of the <em>Melissa</em> space, among others, in this page</a>, and they are free. I haven't tested the profile colors, but being easy to build with the correct data —which in this case is in the public domain— very probably they are fine. The profile tonal scale, as we will see below, is correct.</p>
<p>Using this profile, the previous graphs now look similar to what we got with the other color spaces.</p>
<div class="imcap-frame">
<a class="pic-link cboxed" href="http://www.odelama.com/sandbox/Melissa-correction.png">
<img alt="Testing Melissa correction" class="pic-noborder" data-width="800" src="http://www.odelama.com/sandbox/Melissa-correction.png" title="Testing Melissa correction" />
</a>
<p class="caption-text">ProPhoto issues solved with Melissa RGB space.</p>
</div>
<p>For the sake of completeness I will include the graphs from ProPhoto RGB luminosity masks. Even when in their own RGB space they look exactly as waht we get with the other color spaces. As prescribed, the PS color settings are <code>RGB: ProPhoto RGB</code> and <code>Gray: Gray Gamma 1.8</code>.</p>
<h2 id="prophotorgblightsmasks">ProPhoto RGB Lights Masks</h2>
<div class="slidorion" id="pfRGB-Light_slider">
<div class="slider" data-height="477" data-width="500">
<div class="slide"><img data-src="pfRGB-Light-FFF.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Light-FFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Light-TFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Light-TTT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The masks curves are normalized to have the same height.</div>
<div class="header">Normalized and showing L values</div><div class="content">Showing original tones labeled with equivalent L values.</div>
<div class="header">Normalized with Tones in L</div><div class="content">With the tonal values converted to L. The curves are slightly tilted to the light side.</div>
</div>
<p></p></div>
<h2 id="prophotorgbmid-tonemasks">ProPhoto RGB Mid-tone Masks</h2>
<div class="slidorion" id="pfRGB-Mid_slider">
<div class="slider" data-height="477" data-width="500">
<div class="slide"><img data-src="pfRGB-Mid-FFF.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Mid-FFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Mid-TFT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<div class="slide"><img data-src="pfRGB-Mid-TTT.png" src="http://www.odelama.com/_img/1px.dummy.gif" /></div>
<p></p></div>
<div class="accordion" data-width="280">
<div class="header">Natural Masks</div><div class="content">With their actual (natural) mask values.</div>
<div class="header">Normalized Height</div><div class="content">The masks curves are normalized to have the same height.</div>
<div class="header">Normalized and showing L values</div><div class="content">Showing original tones labeled with equivalent L values.</div>
<div class="header">Normalized with Tones in L</div><div class="content">With the tonal values converted to L. The curves are tilted to the light side.</div>
</div>
<p></p></div>
<h1 id="srgbluminositymasks">sRGB Luminosity Masks</h1>
<p>I haven't made the tests in the <em>sRGB</em> space. However, we have seen the tonal scale of this space is very similar to the <em>Adobe RGB</em> one except for the darkest tones. Therefore we can expect results similar to the ones we have seen in <em>Adobe RGB</em> just with some variation in the darkest tones visible only in the plots with the tones converted to <em>L</em>. Anyway, it would be interesting to have this graphs which will also portray the luminosity masks in the <em>Melissa RGB</em> space.</p>
<p><a href="http://www.odelama.com/sandbox/TonyKuyper-Lumask.zip">Here is the main code used in this analysis</a>.</p>
wiki1http://www.odelama.com/sandbox/wiki1/2013-02-23T23:49:20Z2013-02-23T23:42:53Z
<p>usa page.tmpl</p>
<h1 id="pginadeiniciodewiki1">Página de inicio de Wiki1</h1>
<p>In this post we didn’t try to showcase some eye-candy (although sometimes eye-candy is indeed offered); the designs listed below were selected for their attention to small details. Pretty and colourful header-graphics doesn’t make a good blog.</p>
<p>The blog needs a solid visual structure, a profound hierarchy of site elements; it also has to be able to build some kind of a bridge between the content and its presentation. To do this, you need to think about precision, minimalism and sound use of illustration. These criteria were the ones we’ve used to select the designs listed below. All these aspects make the designs we’ve selected look… well, not always beautiful, but outstanding, almost excellent in their own kind. Mostly it’s the idea the designers used to make the weblog as usable as possible – not the implementation of this idea – which we’ve been after.</p>
Lastest Sandbox blogshttp://www.odelama.com/sandbox/blog/2014-08-31T01:31:48Z2013-02-23T22:33:22Z
<div class="options" data-discussion="no" data-title-icon="calendar-month" style="display:none"></div>
<div class="feedlink">
<a class="feedbutton" href="http://www.odelama.com/sandbox/blog/index.rss" rel="alternate" title="odelama (RSS feed)" type="application/rss+xml">RSS</a>
</div>
<p><article class="blogpost inlinepage">
<section class="inlineheader"></section></article></p>
<div class="inline-title">
<a href="http://www.odelama.com/sandbox/blog/2013/2013-02/2013-02-17/">My Second Test post</a>
</div>
<p><span class="pagedate">
Posted: <time datetime="2013-02-17T13:05:57Z" pubdate="pubdate">Sun Feb 17 13:05:57 2013</time><br />
Last edited: <time datetime="2014-05-30T16:58:12Z">Fri May 30 16:58:12 2014</time><br />
</span></p>
<p></p>
<p><section class="inlinecontent">
Notice that the screenshots we’ve provided may give you a wrong impression about the whole design of the sites; in doubt you should take a closer look at headers, footers, comment-areas, site structure and further site elements. Please also notice that you can click on screenshots to get to the sites from which the screenshots have been taken
</section></p>
<p><footer class="inlinefooter"></footer></p>
<p><nav class="tags"></nav></p>
<p><i class="icon-tag"> </i><a href="http://www.odelama.com/tag/software/" rel="tag">software</a> </p>
<p></p>
<p>
<article class="blogpost inlinepage">
<section class="inlineheader"></section></article></p>
<div class="inline-title">
<a href="http://www.odelama.com/sandbox/blog/2013/2013-02/2013-02-15/">My first test post</a>
</div>
<p><span class="pagedate">
Posted: <time datetime="2013-02-11T21:05:57Z" pubdate="pubdate">Mon Feb 11 21:05:57 2013</time><br />
Last edited: <time datetime="2014-08-24T01:04:39Z">Sun Aug 24 01:04:39 2014</time><br />
</span></p>
<p></p>
<p><section class="inlinecontent">
Grumpy wizards make toxic brew for the evil Queen and Jack. One morning, when Gregor Samsa woke from troubled dreams, he found himself transformed in his bed into a horrible vermin. He lay on his armour-like back, and if he lifted his head a little he could see his brown belly, slightly domed and divided by arches into stiff sections.
</section></p>
<p><footer class="inlinefooter"></footer></p>
<p><nav class="tags"></nav></p>
<p><i class="icon-tag"> </i><a href="http://www.odelama.com/tag/event/" rel="tag">event</a> </p>
<p><i class="icon-tag"> </i><a href="http://www.odelama.com/tag/photography/" rel="tag">photography</a> </p>
<p></p>
<p>
</p>