5.1 Análisis empírico
Este enfoque, criticado por algunos bayesianos radicales, se centra en la escogencia de una estimación \(\hat{\boldsymbol \eta}\) para \(\boldsymbol \eta\), obtenida como el valor que maximiza la verosimilitud marginal previa dada por
\[\begin{align} \tag{5.1} p(\mathbf{Y} \mid \boldsymbol \eta)=\int p(\mathbf{Y} \mid \boldsymbol \theta)p(\boldsymbol \theta\mid \boldsymbol \eta) \ d\boldsymbol \theta \end{align}\]
Por lo tanto todo el andamiaje inferencial está supeditado a la distribución posterior estimada, \(p(\boldsymbol \theta\mid Y,\hat{\boldsymbol \eta})\). Una vez que esta está bien definida, el proceso de estimación puntual, estimación por intervalo y pruebas de hipótesis sigue su curso bayesiano idénticamente como en los capítulos anteriores.
En términos prácticos suponga que se tiene un modelo en dos etapas para cada una de las observaciones. Se asume que existen \(n\) observaciones que, si bien no conforman una muestra aleatoria, conservan la característica de intercambiabilidad y están definidas en los siguientes términos \[\begin{equation*} Y_i \sim p(Y_i \mid \theta_i) \ \ \ \ \ \ \ i=1,\ldots,n \end{equation*}\]
La segunda etapa comienza con la asignación de una distribución11 previa para los parámetros de interés \(\theta_i\). \[\begin{equation*} \theta_i \sim p(\theta_i \mid \boldsymbol \eta) \ \ \ \ \ \ \ i=1,\ldots,n \end{equation*}\]
Nótese que detrás de la asignación de la estructura probabilística para cada uno de los parámetros \(\theta_i\), se supone que éstos últimos determinan una muestra aleatoria de la distribución \(p(\boldsymbol \theta\mid \boldsymbol \eta)\). El objetivo de este enfoque es encontrar estimadores que maximicen la verosimilitud marginal previa la cual, para este caso particular y considerando independencia marginal entre las observaciones y el vector de hiperparámetros, es \[\begin{align} p(Y_i \mid \boldsymbol \eta)&=\int p(Y_i,\theta_i \mid \boldsymbol \eta) \ d\theta_i \notag \\ &=\int p(Y_i \mid \theta_i,\boldsymbol \eta)p(\theta_i \mid \boldsymbol \eta) \ d\theta_i \notag \\ &=\int p(Y_i \mid \theta_i)p(\theta_i \mid \boldsymbol \eta) \ d\theta_i \end{align}\]
De lo anterior, la verosimilitud marginal previa del vector de observaciones dada por la expresión (5.1) queda convertida en \[\begin{align} p(Y \mid \boldsymbol \eta)&=\prod_{i=1}^np(Y_i \mid \boldsymbol \eta) \notag \\ &=\prod_{i=1}^n\int p(Y_i \mid \theta_i)p(\theta_i \mid \boldsymbol \eta) \ d\theta_i \end{align}\]
A continuación, examinamos algunas casos que son de interés para los investigadores.
5.1.1 Modelo Binomial-Beta
Suponga el siguiente modelo binomial (intercambiable) en una primera etapa \[\begin{equation*} Y_i \mid \theta_i \sim Binomial(n_i,\theta_i) \ \ \ \ \ \ \ i=1,\ldots,P. \end{equation*}\]
Para la segunda etapa, se supone una muestra aleatoria proveniente de una misma distribución, tal que \[\begin{equation*} \theta_i \sim Beta(\alpha, \beta) \ \ \ \ \ \ \ i=1,\ldots,p \end{equation*}\]
Puesto que cada \(\theta_i\) se encuentra en el intervalo \((0,1)\), sería apropiado asignarle una distribución Beta.
Análisis preliminar
Es bien sabido que la distribución posterior para cada uno de los parámetros de interés involucrados en el anterior contexto está dada por \[\begin{equation*} \theta_i \mid Y_i \sim Beta(\alpha+Y_i, \beta+n_i-y_i) \end{equation*}\]
para todo \(i=1,\cdots,p\). Sin embargo, como se desconoce totalmente el valor de los hiperparámetros \(\alpha\) y \(\beta\), entonces se debe encontrar una estimación de estos para poder proseguir normalmente con la inferencia bayesiana, pero esta vez enfocados en la estimación de la distribución posterior dada por \[\begin{equation*} \theta_i \mid Y_i \sim Beta(\hat{\alpha}+Y_i, \hat{\beta}+n_i-y_i) \end{equation*}\]
Para tal fin, nótese que la esperanza y la varianza previa de \(\theta_i\) están dadas por las siguientes expresiones \[\begin{align} \tag{5.2} E(\theta_i)&=\frac{\alpha}{\alpha+\beta}\label{ecua2}\\ Var(\theta_i)&=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} \end{align}\]
De donde se tiene que \[\begin{align} \tag{5.3} \alpha=E(\theta_i)(\alpha+\beta) \end{align}\]
y también que \[\begin{align} \tag{5.4} 1-E(\theta_i)=\frac{\beta}{\alpha+\beta} \end{align}\]
Por lo tanto \[\begin{align} \tag{5.5} \beta=(1-E(\theta_i))(\alpha+\beta) \end{align}\]
Reemplazando (5.3) y (5.4) en (5.2), se concluye que: \[\begin{align*} Var(\theta_i)&=\frac{E(\theta_i)(1-E(\theta_i))}{(\alpha+\beta+1)} \end{align*}\]
Por consiguiente, \[\begin{align} \alpha+\beta=\frac{E(\theta_i)(1-E(\theta_i))}{Var(\theta_i)}-1 \end{align}\]
Con el anterior razonamiento, es posible encontrar los estimadores de los hiperparámetros utilizando el método frecuentista de los momentos, los cuales corresponden a \[\begin{align} \widehat{\alpha+\beta}&=\frac{\bar{Y}(1-\bar{Y})}{S^2}-1 \end{align}\]
Donde \(\bar{Y}\) y \(S^2\) es el promedio y la varianza de las cantidades \(Y_1/n_1, Y_2/n_2,\ldots, Y_P/n_P\), respectivamente. Ahora, teniendo en cuenta (5.3) y (5.5), se tiene que: \[\begin{align} \hat{\alpha}&=\left(\frac{\bar{Y}(1-\bar{Y})}{S^2}-1\right)\bar{Y}\\ \hat{\beta}&=\left(\frac{\bar{Y}(1-\bar{Y})}{S^2}-1\right)(1-\bar{Y}) \end{align}\]
Con las anteriores estimaciones es posible ahora conectarlas a la distribución posterior de \(\theta_i\).
Análisis legítimo
Según A. Gelman et al. (2003, pág. 119), el anterior análisis no implica simplemente un punto de partida que da pie a la exploración de la idea de la estimación de los parámetros de la distribución posterior y, de ninguna manera, constituye un cálculo bayesiano, puesto que no está basado en ningún modelo de probabilidad. Sin embargo, el análisis empírico de esta situación, hace uso del esperanza y varianza condicional de la distribución beta de los parámetros \(\theta_i\) (\(i=1,\ldots, p\)).
Para realizar este tipo de análisis, vamos a suponer que contamos con una variable \(Y\), distribuida de forma binomial en \(n\) ensayos y con probabilidad de éxito \(\theta\). De esta manera, se tiene que el primer momento está dado por
\[\begin{align} E_{binom}\left(\frac{Y}{n}\right)&=E_{beta}\left(E_{binom}\left(\frac{Y}{n} \mid \theta\right)\right) \notag \\ &=E_{beta}\left(\theta\right) \notag \\ \tag{5.6} &=\frac{\alpha}{\alpha+\beta} \end{align}\]
Por otro lado, se tiene que la varianza, que es función del primer y segundo momento, está dada por
\[\begin{align*} Var_{binom}\left(\frac{Y}{n}\right) &=E_{beta}\left(Var_{binom}\left(\frac{Y}{n} \mid \theta\right)\right) + Var_{beta}\left(E_{binom}\left(\frac{Y}{n} \mid \theta\right)\right) \\ &=E_{beta}\left( \frac{1}{n}\theta(1-\theta)\right) + Var_{beta}\left(\theta\right) \\ &=\frac{1}{n}E_{beta}(\theta) - \frac{1}{n}E_{beta}(\theta^2)+ Var_{beta}\left(\theta\right) \\ &=\frac{1}{n}E_{beta}(\theta) - \frac{1}{n}Var_{beta}(\theta)- \frac{1}{n}(E_{beta}\theta)^2+ Var_{beta}\left(\theta\right) \\ &=\frac{n-1}{n}Var_{beta}(\theta) + \frac{1}{n}E_{beta}(\theta)(1-E_{beta}(\theta)) \\ &=\frac{n-1}{n}\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)} + \frac{1}{n}\frac{\alpha\beta}{(\alpha+\beta)^2} \\ &=\frac{1}{n}\frac{\alpha}{\alpha+\beta}\frac{\beta}{\alpha+\beta} \left(\frac{n-1}{\alpha+\beta+1}+1\right) \\ &=\frac{1}{n}E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right) \left(\frac{n-1}{\alpha+\beta+1}+1\right) \end{align*}\]
De esta última expresión, y despejando \(\alpha +\beta\), se tiene que
\[\begin{align} \alpha+\beta&=\frac{(n-1)E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)} {nVar_{binom}\left(\frac{Y}{n}\right)- E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)}-1\notag \\ &=\frac{E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)-Var_{binom}\left(\frac{Y}{n}\right)} {Var_{binom}\left(\frac{Y}{n}\right)-\frac{1}{n}E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)} \end{align}\]
Ahora, despejando \(\alpha\) de la expresión (5.6) se tiene que \[\begin{align} \alpha=E_{binom}\left(\frac{Y}{n}\right)\frac{E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right) -Var_{binom}\left(\frac{Y}{n}\right)}{Var_{binom}\left(\frac{Y}{n}\right) -\frac{1}{n}E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)} \end{align}\]
Además, también despejando \(\beta\) de (5.6) se tiene que \[\begin{align} \beta&=\frac{\alpha\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)}{E_{binom}\left(\frac{Y}{n}\right)} \notag \\ &=\frac{E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)-Var_{binom}\left(\frac{Y}{n}\right)} {Var_{binom}\left(\frac{Y}{n}\right)-\frac{1}{n}E_{binom}\left(\frac{Y}{n}\right)\left(1-E_{binom}\left(\frac{Y}{n}\right)\right)}\left(1-E_{binom}\left(\frac{Y}{n}\right)\right) \end{align}\]
El anterior enfoque nos ha llevado a poder expresar los parámetros de interés en términos de \(E_{binom}\left(\frac{Y}{n}\right)\), \(Var_{binom}\left(\frac{Y}{n}\right)\) y \(n\). Una vez que podamos estimar las anteriores cantidades, es posible realizar la inferencia bayesiana empírica de la manera correcta. Para lo anterior, es necesario observar la naturaleza de las variables que, aunque no representan una muestra aleatoria, sí son una sucesión de variables aleatorias intercambiables. Por lo tanto, teniendo en cuenta que la inferencia se realiza con las cantidades \(Y_1/n_1, Y_2/n_2, \ldots, Y_P/n_P\), es posible proponer los siguientes estimadores \[\begin{align} \widehat{E}_{binom}\left(\frac{Y}{n}\right)&=\bar{Y} \\ \widehat{Var}_{binom}\left(\frac{Y}{n}\right)&=S^2 \\ \widehat{n}&=\frac{1}{p}\sum_{i=1}^pn_i \end{align}\]
Con base en lo anterior, las estimaciones empíricas de los parámetros \(\alpha\) y \(\beta\) son \[\begin{align} \hat{\alpha}=\bar{Y}\left(\frac{\bar{Y}\left(1-\bar{Y}\right) -S^2}{S^2-\frac{1}{\hat{n}}\bar{Y}\left(1-\bar{Y}\right)}\right) \end{align}\]
y \[\begin{align} \hat{\beta}&=(1-\bar{Y})\frac{\bar{Y}\left(1-\bar{Y}\right)-S^2}{S^2-\frac{1}{\hat{n}}\bar{Y}\left(1-\bar{Y}\right)} \end{align}\]
respectivamente. Cuando la cantidad de ensayos \(n_i\) es diferente en cada experimento, existen otras formas de obtener estimaciones para los parámetros \(\alpha\) y \(\beta\) (Carlin y Louis 1996, pág. 81).
library(pscl)
data(EfronMorris)
attach(EfronMorris)
y <- p # Porcentaje de bateo de los 18 jugadores
y## [1] 0.346 0.298 0.276 0.222 0.273 0.270 0.263 0.210 0.269 0.230 0.264 0.256
## [13] 0.303 0.264 0.226 0.285 0.316 0.200y.bar <- mean(y)
S2 <- var(y)
n.hat <- mean(n)
alfa <- y.bar * (y.bar * (1 - y.bar) - S2)/
(S2 - y.bar * (1 - y.bar)/n.hat)
beta <- (1 - y.bar) * (y.bar * (1 - y.bar) - S2)/
(S2 - y.bar * (1 - y.bar)/n.hat)
alfa## [1] 56.99536beta## [1] 158.0364De donde podemos concluir que la distribución previa para cada \(\theta_i\) es la distribución \(Beta(57, 158)\), Nótese que la esperanza de esta distribución es 0.26, que coincide con el porcentaje de bateo promedio de los datos de los 18 jugadores. Ahora, podemos calcular los parámetros de la distribución para cada \(\theta_i\) con \(i=1,\cdots,18\) como sigue:
alfa.new <- alfa + p * n
beta.new <- beta + (1 - p) * n
head(alfa.new)## [1] 183.9774 183.9434 200.7914 118.0454 171.1094 182.8154head(beta.new)## [1] 398.0544 457.0884 535.2404 371.9864 461.9224 498.2164De esta manera, podemos realizar inferencias para cualquier \(\theta_i\). Por ejemplo, para primer jugador, Roberto Clemente, la distribución posterior para su porcentaje de bateo es \(Beta(184, 398)\); por consiguiente la estimación para el porcentaje de bateo de este jugador es \(184/(184+398)=0.3162\) y su intervalo de credibilidad será \((0.279,0.354)\).
Referencias
En esta etapa la distribución previa no está completamente especificada puesto que se desconocen los hiperparámetros que la indexan.↩︎