## Warning: package 'rjags' was built under R version 4.2.3
## Loading required package: coda
## Linked to JAGS 4.3.1
## Loaded modules: basemod,bugs
Con frecuencia nos topamos con el problema de “estimar” la distribución que sigue el comportamiento de una determinada característica en un grupo diverso de individuos de los que contamos con observaciones en forma individual, pero de los que sabemos que pertenecen a una “colectividad”, es decir, sabemos que cada individuo en lo particular tiene un comportamiento específico pero también que podemos esperar cierta similitud al resto de los miembros del grupo.
En este contexto cobran relevancia los modelos jerárquicos, también llamados modelos multinivel.
En términos generales, los modelos jerárquicos no son otra cosa que una forma de modelar el comportamiento de un conjunto de variables aleatorias, condicionando este comportamiento al hecho de que todas estas variables aleatorias están gobernadas por un conjunto de parámetros que impactan en el comportamiento individual de todas por igual.
Un ámbito en el que podemos aplicar el uso de este tipo de modelos es el de la detección de observaciones atípicas aplicado a las operaciones de los clientes de un intermediario financiero (e.g., para el alertamiento de operaciones sospechosas). Este caso se ajusta a la perfección al contexto anteriormente descrito ya que queremos determinar si las operaciones de un determinado cliente pueden considerar anómalas, pero esto se da en función del comportamiento del cliente en lo individual en combinación con el comportamiento del cliente como miembro de un grupo (el grupo de clientes de la entidad o un subgrupo o categoría de clientes de la misma).
Considérese un conjunto de experimentos \(j = 1, \dots, J\), en los que cada experimento \(j\) tiene asociado un vector de observaciones \(y_j\) y un vector de parámetros \(\theta_j\), con una función de verosimilitud \(p(y_j|\theta_j).\)
Si no se cuenta con información para diferencias a los parámetros \(\theta_j\) y no se puede llevar a cabo ningún tipo de ordenamiento o agrupación de los parámetros, entonces se debe asumir simetría en su distribución previa. La simetría es representada probabilísticamente por la intercambiabilidad. La forma más simple de una distribución intercambiable trata a los parámetros \(\theta_j\) como muestras independientes de una distribución previa poblacional gobernada por un parámetro \(\phi\):
\[p(\theta|\phi)= \prod_{j = 1}^J p(\theta_j|\phi)\]
Y como generalmente se asume desconocida:
\[p(\theta) = \int \left[ \prod_{j = 1}^J p(\theta_j|\phi) \right]p(\phi)d\phi.\]
Estadísticamente, el modelo de la mezcla independiente e idénticamente distribuida caracteriza a los parámetros como obtenidos de una muestra de una “superpoblación” común determinada por los parámetros desconocidos .
La parte “jerárquica” de estos modelos se encuentra en el hecho de que es desconocida y, por lo tanto, tiene su propia distribución previa. De esta manera:
\[p(\phi, \theta) = p(\theta|\phi)p(\phi).\]
Y la probabilidad posterior conjunta:
\[p(\phi, \theta|y) = \frac{p(\theta,\phi,y)}{p(y)} \propto p(y|\theta, \phi)p(\phi,\theta)\].
Para poder derivar las distribuciones condicionales y marginales en forma analítica es necesario seguir los siguientes pasos:
Escribir la distribución posterior conjunta \(p(\phi, \theta|y)\) como el producto de la distribución hyperprevia \(p(\phi)\), la distribución de la población \(p(\theta|\phi)\) y la verosimilitud \(p(y|\phi, \theta)\): \(p(\theta|\phi,y)\).
Determinar analíticamente la densidad condicional posterior de \(\theta\), dados los hyperparámetros \(\phi\) y las observaciones y.
Estimar \(\phi\) obteniendo su distribución marginal posterior \(p(\phi|y)\).
Computacionalmente:
Sacamos muestras del vector de hyperparámetros \(\phi\) de su distribución marginal posterior \(p(\phi|y)\).
Sacamos muestras del vector de parámetros \(\theta\) de su distribución marginal posterior \(p(\theta|\phi, y)\), usando las muestras obtenidas de \(\phi\).
Obtener valores predictivos \(\tilde{y}\) the la distribución predictiva posterior dados los valores obtenidos de \(\theta\).
Tenemos 71 experimentos realizados en diferentes muestras de ratas. Los experimentos se suponen provenientes de distribuciones binomiales independientes tales que:
\[y_j|\theta_j \sim Bin(n_j,\theta_j)\]
con el número de ratas de cada experimento, \(n_j\), conocido. Los parámetros \(\theta_j\) se asumen como muestras independientes provenientes de una distribución beta:
\[\theta_j|\alpha,\beta \sim Beta(\alpha,\beta).\]
Obtención de las distribuciones posteriores: conjunta, condicional y marginal
\[p(\theta, \alpha, \beta|y) \propto p(y|\theta, \alpha, \beta)p(\theta, \alpha, \beta) \propto p(y|\theta, \alpha, \beta)p(\theta|\alpha, \beta)p(\alpha,\beta);\]
\[p(y_j|\theta,\alpha,\beta) = \binom{n_j}{y_j}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j};\]
Por lo que, considerando que las \(y_j\) son independientes:
\[p(y|\theta,\alpha,\beta) = \prod_{j=1}^{71}\binom{n_j}{y_j}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j}.\]
Sustituyendo:
\[p(\theta, \alpha, \beta|y) \propto p(\theta|\alpha,\beta)p(\alpha,\beta)\prod_{j=1}^{71}\binom{n_j}{y_j}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j} \propto p(\theta|\alpha,\beta)p(\alpha,\beta)\prod_{j=1}^{71}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j}.\]
Ahora, dado que asumimos que los parámetros era independientes, sabemos que:
\[p(\theta|\alpha,\beta) = \prod_{j=1}^{71}p(\theta_j|\alpha,\beta);\]
\[\theta_j|\alpha,\beta \sim Beta(\alpha,\beta);\]
\[p(\theta_j|\alpha,\beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}{\theta_j}^{\alpha-1}(1-{\theta_j})^{\beta-1} \Rightarrow\]
\[p(\theta|\alpha,\beta) = \prod_{j=1}^{71}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}{\theta_j}^{\alpha-1}(1-{\theta_j})^{\beta-1} \Rightarrow;\]
\[p(\theta, \alpha, \beta|y) \propto p(\alpha,\beta)\prod_{j=1}^{71}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}{\theta_j}^{\alpha-1}(1-{\theta_j})^{\beta-1}\prod_{j=1}^{71}\theta_j^{y_j}(1-\theta_j)^{n_j-y_j} =\]
\[p(\alpha,\beta)\left[\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right]^{71}\prod_{j=1}^{71}\theta_j^{y_j+\alpha-1}(1-\theta_j)^{n_j-y_j+\beta-1}\]
Ahora, para la obtención de la marginal posterior de \(\theta\), observamos que:
\[p(\theta,\alpha,\beta|y) = p(\theta|\alpha,\beta,y)p(\alpha,\beta|y).\]
En otras palabras:
\[p(\theta|\alpha,\beta,y) \propto p(\theta,\alpha,\beta|y) = p(\alpha,\beta)\left[\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\right]^{71}\prod_{j=1}^{71}\theta_j^{y_j+\alpha-1}(1-\theta_j)^{n_j-y_j+\beta-1} \propto\]
\[\prod_{j=1}^{71}\theta_j^{y_j+\alpha-1}(1-\theta_j)^{n_j-y_j+\beta-1}.\]
Si ahora observamos la parte interna del producto:
\[\theta_j^{y_j+\alpha-1}(1-\theta_j)^{n_j-y_j+\beta-1}\]
podemos claramente ver que se trata del kernel de una distribución \(Beta(\alpha+y_j,\beta+n_j-y_j)\), por lo que podemos completar la expresión analítica de la densidad posterior conjunta de \(\theta\):
\[p(\theta|\alpha,\beta,y) = \prod_{j=1}^{71}\frac{\Gamma(\alpha+\beta+n_j)}{\Gamma(\alpha+y_j)\Gamma(\beta+n_j-y_j)}\theta_j^{y_j+\alpha-1}(1-\theta_j)^{n_j-y_j+\beta-1}.\]
\[ \begin{aligned} y_{ij} &= \beta_{0j} + \beta_{1j} x_{1ij} + \beta_{2j} x_{2ij} + \epsilon_{ij} & \epsilon_{ij} \sim N(0,\sigma^2)\\ \beta_{pj} &\sim N(B_{pj}, \tau^{2}_p) & p = 1,2\\ B_{pj} &= \gamma_{p0} + \gamma_{p1}z_j &\\ \sigma &\sim U(0,1000) &\\ \tau_p &\sim U(0,1000) &\\ \end{aligned} \]
## gam0.1. gam0.2. gam0.3. gam1.1. gam1.2. gam1.3.
## -1.166769867 1.258455770 0.795596145 0.224221103 -0.001153464 -0.024365477
## s1 t1 tau.1. tau.2. tau.3.
## 0.747093920 1.793014555 0.521724658 0.078944927 0.026963460
## gam0.1. gam0.2. gam0.3. gam1.1. gam1.2. gam1.3.
## 0.268671459 0.101902383 0.038713458 0.016946325 0.006721013 0.002501302
## s1 t1 tau.1. tau.2. tau.3.
## 0.011976333 0.057393944 0.049660897 0.048668588 0.015690908
El monto de las operaciones de un cliente tiene una densidad \(p(y_j | \theta_j, x_j)\). A su vez, los parámetros \(\theta_j\) pueden considerarse como una observación proveniente de \(p(\theta_j | \phi_i)\); finalmente, los parámetros \(\phi_i\) pueden considerarse como provenientes de \(p(phi_i | \rho).\)
Pendientes: distribución marginal posterior de \(\alpha, \beta\). Página 128 Gelman et al. (2004)
Una variable aleatoria \(X\) tiene una distribución binomial con parámetros \(n\) y \(p\) si su distribución discreta es:
\[f(x|n,p) = \begin{cases} \binom{n}{x} p^x q^{n-x} & \text{Para x = 0, 1, ..., n,}\\ 0 & \text{E.O.C.} \end{cases} \]
Donde \(n \in \mathbb{Z}^+\) y \(0 \leq p \leq 1\).
Se dice que una variable aleatoria X tiene una distribución beta con parámetros \(\alpha > 0\) y \(\beta > 0\) si su función de densidad es la siguiente:
\[f(x | \alpha, \beta) = \begin{cases} \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} x^{\alpha - 1}(1 - x)^{\beta - 1} & \text{Para 0 < x < 1}\\ 0 & \text{E.O.C.} \end{cases} .\]
Donde \(\Gamma(\alpha), \alpha > 0,\) es la función Gamma:
\[\Gamma(\alpha) = (\alpha - 1)\Gamma(\alpha - 1).\]