1 Introducción

El tipo más simple de datos multinivel tiene dos niveles, en los que un nivel consiste de grupos y el otro consiste en unidades dentro de los grupos.

Se denota con \(y_{i,j}\) la observación de la unidad \(i\) en el grupo \(j\), para \(i = 1,\ldots,n_j\) y \(j=1,\ldots,m\), donde \(m\) es el número de grupos y \(n = \sum_{j=1}^m n_j\) es el tamaño de la muestra.

El conjunto de datos es \(\boldsymbol{y} = (\boldsymbol{y}_1,\ldots,\boldsymbol{y}_m)\), donde \(\boldsymbol{y}_j=(y_{1,j},\ldots,y_{n_j,j})\) son las observaciones asociadas con el grupo \(j\), para \(j=1,\ldots,m\).

2 Modelo jerarquico general de dos niveles

Caracterización dentro de los grupos: \[ y_{i, j}\mid\boldsymbol{\phi}_j \stackrel{\text {iid}}{\sim} p\left(y_{i,j} \mid \boldsymbol{\phi}_{j}\right) \]
Caracterización entre los grupos: \[ \boldsymbol{\phi}_{j}\mid\boldsymbol{\psi} \stackrel{\text {iid}}{\sim} p\left(\boldsymbol{\phi}_j \mid \boldsymbol{\psi}\right) \]
Distribución previa: \[ \boldsymbol{\psi} \sim p(\boldsymbol{\psi}) \]
Los parámetros del modelo son \(\boldsymbol{\theta}=(\boldsymbol{\phi}_1,\ldots,\boldsymbol{\phi}_m,\boldsymbol{\psi})\).

3 Modelo jerarquico Normal de dos niveles

Un modelo popular para describir la heterogeneidad de las medias en varias poblaciones es el modelo jerárquico normal, en el cual la variabilidad dentro y entre se modela usando distribuciones normales:

Caracterización dentro de los grupos: \[ y_{i, j}\mid\theta_j,\sigma^2 \stackrel{\text {iid}}{\sim} \textsf{N}\left(\theta_j,\sigma^2\right) \]
Caracterización entre los grupos: \[ \theta_{j}\mid\mu,\tau^2 \stackrel{\text {iid}}{\sim} \textsf{N}\left(\mu,\tau^2\right) \]
Distribución previa: \[ p(\sigma^2,\mu,\tau^2) = p(\sigma^2)\,p(\mu)\,p(\tau^2) \] donde \[ \sigma^2\sim\textsf{GI}\left(\tfrac{\nu_0}{2},\tfrac{\nu_0\,\sigma^2_0}{2}\right)\qquad\mu\sim\textsf{N}(\mu_0,\gamma_0^2)\qquad\tau^2\sim\textsf{GI}\left(\tfrac{\eta_0}{2},\tfrac{\eta_0\,\tau^2_0}{2}\right) \]
Los parámetros del modelo son \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_m,\sigma^2,\mu,\tau^2)\).
Los hiper-parámetros del modelo son \((\nu_0,\sigma^2_0,\mu_0,\gamma_0^2,\eta_0,\tau^2_0)\).

3.1 Estimación

Construir un muestreador de Gibbs para obtener muestras de la distribución posterior \(p(\boldsymbol{\theta}\mid\boldsymbol{y})\).

Distribución posterior: \[ \begin{aligned} p(\boldsymbol{\theta} \mid \boldsymbol{y}) &\propto p(\boldsymbol{y} \mid \boldsymbol{\theta})\, p(\boldsymbol{\theta}) \\ &= p(\boldsymbol{y}_1,\ldots,\boldsymbol{y}_m \mid \theta_1,\ldots\theta_m,\sigma^2)\times p(\theta_1,\ldots,\theta_m\mid\mu,\tau^2)\times p(\sigma^2,\mu,\tau^2) \\ &=\prod_{j=1}^m\prod_{i=1}^{n_j} \mathrm{N}\left(y_{i, j} \mid \theta_{j}, \sigma^{2}\right) \times \prod_{j=1}^m \mathrm{N}\left(\theta_{j} \mid \mu, \tau^{2}\right) \times \operatorname{GI}\left(\sigma^{2} \mid \tfrac{\nu_{0}}{2}, \tfrac{\nu_{0}\,\sigma_{0}^{2}}{2} \right) \times \mathrm{N}\left(\mu \mid \mu_{0}, \gamma_{0}^{2}\right) \times \operatorname{GI}\left(\tau^{2} \mid \tfrac{\eta_{0}}{2}, \tfrac{\eta_{0}\,\tau_{0}^{2}}{2}\right) \end{aligned} \]
Distribuciones condicionales completas: \[ \begin{aligned} p\left(\theta_{j} \mid \text { resto }\right) &=\textsf{N}\left(\theta_{j} \,\Big|\, \frac{\mu / \tau^{2} + n_{j} \bar{y}_{j} / \sigma^{2}}{1 / \tau^{2} + n_{j} /\sigma^{2}}, \frac{1}{1 / \tau^{2} + n_{j} /\sigma^{2}}\right) \\ p\left(\sigma^{2} \mid \text { resto }\right) &=\textsf{GI}\left(\sigma^{2} \,\Big|\, \frac{\nu_{0}+\sum_{j=1}^m n_{j}}{2}, \frac{\nu_{0} \sigma_{0}^{2}+\sum_{j=1}^m \sum_{i=1}^{n_j}\left(y_{i, j}-\theta_{j}\right)^{2}}{2}\right) \\ p(\mu \mid \text { resto }) &=\textsf{N}\left(\mu \,\Big|\, \frac{\mu_{0} / \gamma_{0}^{2} + m \bar{\theta} / \tau^{2}}{1 / \gamma_{0}^{2} + m / \tau^{2}}, \frac{1}{1 / \gamma_{0}^{2} + m / \tau^{2}}\right) \\ p\left(\tau^{2} \mid \text { resto }\right) &=\textsf{ GI }\left(\tau^{2} \,\Big|\, \frac{\eta_{0}+m}{2}, \frac{\eta_{0} \tau_{0}^{2}+\sum_{j=1}^m\left(\theta_{j}-\mu\right)^{2}}{2}\right) \end{aligned} \]

4 Ejemplo: Puntajes de matemáticas

Estudio observacional en US en 2002 que incluyó a estudiantes de décimo grado de 100 colegios públicos con una matrícula en décimo grado de 400 estudiantes o más.

Los puntajes de matemáticas se basaron en los resultados de un examen nacional estandarizado, diseñado para producir una media de 50 y una desviación estándar de 10.

4.1 Datos

# directorio de trabajo
setwd("C:/Users/Juan Camilo/Dropbox/UN/estadistica_bayesiana_2021_2/")
# cargar algunas funciones
source("funciones_auxiliares.R")
# cargar data
load("math.RData")

# NOTACION
# ids  : identificador de los colegios (c)
# J    : numero de grupos (colegios)
# m    : numero de grupos (colegios)
# n    : numero de estudiantes en cada colegio (c)
# YM   : puntajes  de los estudiantes (matrix)
# Y    : puntajes  de los estudiantes (list)
# ybar : promedios de los puntajes (c)
# ymed : medianas  de los puntajes (c)
# s2   : varianzas de los puntajes (c)
# sv   : varianzas de los puntajes (c)

4.2 Análisis explotarotio

# representacion de los puntajes brutos
par(mfrow=c(1,1),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(c(1,J), range(Y), type="n", ylab="Puntaje", xlab="Ranking (promedio muestral)")
for(l in 1:J) {
  j<-order(ybar)[l]
  points( rep(l, n[j]), Y[[j]], pch=16, cex=.5 )
  segments( l, min(Y[[j]]), l, max(Y[[j]]) )
}
abline(h=mean(ybar))

Es frecuente que los grupos con promedios muestrales muy altos o muy bajos correspondan a aquellos grupos con tamaños muestrales bajos, ya que \(\textsf{Var}(\bar{y}) = \sigma^2/n\).

par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0)) 
hist(ybar, freq = F, main="", xlab="Promedio", ylab="Densidad", border="gray90", col="gray90")
plot(n, ybar, xlab="Tamaño del grupo", ylab="Promedio", pch = 16, col="gray80")

4.3 Distribución previa

El examen de matemáticas se diseñó para dar una varianza nacional de 100. Dado que esta varianza incluye la varianza tanto dentro de la escuela como entre escuelas, la varianza dentro de la escuela debe ser como máximo de 100, por lo que se toma \(\sigma_0^2 = 100\).
Se concentra débilmente la distribución previa de \(\sigma^2\) alrededor de \(\sigma_0^2 = 100\) tomando \(\nu_0 = 1\).
De manera similar, la varianza entre escuelas no debe ser superior a 100, por lo que se toma \(\tau_0^2 = 100\) y \(\eta_0 = 1\).
La media nacional de todas las escuelas es 50. Aunque la media de las grandes escuelas públicas urbanas puede ser diferente del promedio nacional, no debería diferir demasiado. Por eso se toma \(\mu_0 = 50\) y \(\gamma_0 = 5\), de forma que la probabiliad previa de que \(\mu_0\) este en el intervalo \((40, 60)\) es aprox. 95%.

# hiper-parametros para obtener previas no informativos
nu0  <- 1  
s20  <- 100
eta0 <- 1  
t20  <- 100
mu0  <- 50 
g20  <- 25

4.4 Estimación

# 0. setup
S     <- 50000                            # numero de iteraciones
ncat  <- floor(S/10)                      # progreso
THETA <- matrix(data=NA, nrow=S, ncol=m)  # almacenar thetas 
MST   <- matrix(data=NA, nrow=S, ncol=3)  # almacenar mu, sigma^2, tau^2

# 1. valores inciales
theta  <- ybar         # rep(50, m)
sigma2 <- mean(sv)     # 100
mu     <- mean(theta)  # 50
tau2   <- var(theta)   # 100

# 2. MCMC
set.seed(1)
for (s in 1:S) {
        # 2.1 actualizar theta
        for(j in 1:m) {
                vtheta   <- 1/(n[j]/sigma2+1/tau2)
                etheta   <- vtheta*(ybar[j]*n[j]/sigma2+mu/tau2)
                theta[j] <- rnorm(1,etheta,sqrt(vtheta))
        }
        # 2.2 actualizar sigma^2
        nun <- nu0+sum(n)
        ss  <- nu0*s20
        for (j in 1:m) {
                ss <- ss+sum((Y[[j]]-theta[j])^2)
        }
        sigma2 <- 1/rgamma(1,nun/2,ss/2)
        # 2.3 actualizar mu
        vmu <- 1/(m/tau2+1/g20)
        emu <- vmu*(m*mean(theta)/tau2 + mu0/g20)
        mu  <- rnorm(1,emu,sqrt(vmu)) 
        # 2.4 actualizar tau2
        etam <- eta0+m
        ss   <- eta0*t20 + sum( (theta-mu)^2 )
        tau2 <- 1/rgamma(1,etam/2,ss/2)
        # 2.5 almacenar valores
        THETA[s,] <- theta
        MST[s,]   <- c(mu,sigma2,tau2)
        # 2.6 progreso
        if (s%%ncat == 0) cat(100*round(s/S, 1), "% completed ... \n", sep = "" )
}

## 10% completed ... 
## 20% completed ... 
## 30% completed ... 
## 40% completed ... 
## 50% completed ... 
## 60% completed ... 
## 70% completed ... 
## 80% completed ... 
## 90% completed ... 
## 100% completed ...

# FINAL DEL ALGORITMO

4.5 Convergencia

# cadenas
par(mfrow=c(3,1),mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
plot(THETA[,1], xlab="Iteración", type ="l", ylab=expression(theta[1]))
plot(THETA[,2], xlab="Iteración", type ="l", ylab=expression(theta[2]))
plot(THETA[,3], xlab="Iteración", type ="l", ylab=expression(theta[3]))

# autocorrelaciones
par(mfrow=c(1,3), mar=c(3,3,2.5,1),mgp=c(1.75,.75,0))
acf(MST[,1], main=expression(mu))
acf(MST[,2], main=expression(sigma^2))
acf(MST[,3], main=expression(tau^2))

# tamaños efectivos de muestra
neff_MST <- coda::effectiveSize(MST)
neff_MST

##     var1     var2     var3 
## 36653.64 43956.19 25601.83

neff_THETA <- coda::effectiveSize(THETA)
summary(neff_THETA)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   43470   48672   50000   49232   50000   51925

# errores estandar de MC : DE(parametro)/sqrt( n_eff )
MCERR <- apply(MST,2,sd)/sqrt(neff_MST)
MCERR

##        var1        var2        var3 
## 0.002824017 0.013213703 0.027452664

TMCERR <- apply(THETA,2,sd)/sqrt(neff_THETA)
summary(TMCERR)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.006913 0.007853 0.008486 0.008891 0.009587 0.016596

4.6 Inferencia

# distribuciones posteriores para mu, sigma^2 and tau^2
par(mfrow=c(1,3),mar=c(2.75,2.75,.5,.5),mgp=c(1.7,.7,0))
# mu
plot(density(MST[,1],adj=2),xlab=expression(mu),main="",lwd=2,
ylab=expression(paste(italic("p("),mu,"|",italic(y[1]),"...",italic(y[m]),")")))
abline( v=quantile(MST[,1],c(.025,.5,.975)),col=c(2,4,2),lty=c(3,2,3) )
# sigma^2
plot(density(MST[,2],adj=2),xlab=expression(sigma^2),main="", lwd=2,
ylab=expression(paste(italic("p("),sigma^2,"|",italic(y[1]),"...",italic(y[m]),")")))
abline( v=quantile(MST[,2],c(.025,.5,.975)),col=c(2,4,2),lty=c(3,2,3) )
# tau^2
plot(density(MST[,3],adj=2),xlab=expression(tau^2),main="",lwd=2,
ylab=expression(paste(italic("p("),tau^2,"|",italic(y[1]),"...",italic(y[m]),")")))
abline( v=quantile(MST[,3],c(.025,.5,.975)),col=c(2,4,2),lty=c(3,2,3) )

# proporcion intravarianza
eta <- 100*MST[,2]/(MST[,2] + MST[,3])
plot(density(eta,adj=2), xlab=expression(eta),main="",lwd=2, ylab=expression(paste(italic("p("),eta,"|",italic(y[1]),"...",italic(y[m]),")")))
abline(v=quantile(eta,c(.025,.5,.975)),col=c(2,4,2),lty=c(3,2,3))

# medias posteriores
mean(MST[,1])

## [1] 48.12686

mean(sqrt(MST[,2]))

## [1] 9.209248

mean(sqrt(MST[,3]))

## [1] 4.97293

Aproximadamente el 99% de los puntajes dentro de una escuela están dentro de 4 × 9.21 = 37 puntos entre sí, mientras que el 99% de los puntajes promedio de las escuelas están dentro de 4 × 4.97 = 20 puntos el uno del otro.
Además se observa que la media global es aproximadamente 50, lo cual coincide con el diseño de la prueba.

# encogimiento (shrinkage)
par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
theta.hat<-apply(THETA,2,mean)
plot(ybar,theta.hat,xlab=expression(bar(italic(y))),ylab=expression(hat(theta)))
abline(0,1)
plot(n,ybar-theta.hat,ylab=expression( bar(italic(y))-hat(theta) ),xlab="Tamaño de la muestra")
abline(h=0)

La relación sigue aproximadamente una línea con una pendiente menor que uno, lo que indica que los valores altos de \(\bar{y}_j\) corresponden a valores ligeramente menos altos \(\hat{\theta}_j\).
Los grupos con tamaños de muestra bajos son los que más se “reducen” o “encogen” (shrunk) mientras que los grupos con tamaños de muestra grandes casi no se “reducen” en absoluto.
Cuanto mayor sea el tamaño de la muestra de un grupo, se tiene más información para ese grupo y menos información es necesita “tomar prestada” del resto de la población.

5 Ranking

# ranking de acuerdo con theta
theta.order<-order(theta.hat)
theta.order[1:20]

##  [1]  5 72 49  6 78 46 17 74  7 10 82 50 60 64 55 53 70 84 57 19

# ranking de acuerdo con el promedio muestral
ybar.order<-order(ybar)
ybar.order[1:20]

##  [1]  5 72 17 49 82  6 78 46 74  7 10 50 60 53 55 64 19 27 70 84

# estmaciones
theta.hat[c(46,82)]

## [1] 41.30664 42.58871

ybar[c(46,82)]

## [1] 40.17619 38.76400

# Tamaños de muetra
n[c(46,82)]

## [1] 21  5

# Pr (theta_46 < theta_82 )
mean( THETA[,46] < THETA[,82])

## [1] 0.63572

# ranking de colegios 46 y 82
par(mfrow=c(1,1),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot(density(THETA[,46],adj=2),col="black",xlim=range(c(Y[[46]],Y[[82]],THETA[,c(46,82)])),lwd=3,main="",xlab="Puntaje de matemáticas",ylim=c(-.05,.21),ylab="Densidad",yaxt="n")
axis(side=2,at=c(0,0.10,0.20))
lines(density(THETA[,82],adj=2),col="gray",lwd=3)
abline(h=0)
#
points(Y[[46]],rep(-0.01666,n[46]), col="black",pch=16)
points(ybar[46],-.01666,col="black",pch=16 ,cex=1.5)
abline(h=-.01666,col="black")
#
points(Y[[82]],rep(-0.0333,n[82]), col="gray",pch=16)
points(ybar[82],-.0333,col="gray",pch=16 ,cex=1.5)
abline(h=-.0333,col="gray")
#
segments(mean(MST[,1]), 0,mean(MST[,1]),1,lwd=2,lty=2 )
#
legend(52.5,.15, legend=c("Colegio 46","Colegio 82",expression(paste("E[", mu,"|",italic(y[1]),"...",italic(y[m]),"]"))),lwd=c(3,3),lty=c(1,1,2),col=c("black","gray"),bty="n")

Los datos brutos para las dos escuelas se muestran en gráficos de puntos debajo de las densidades posteriores; los puntos grandes representan las medias del grupo.
La densidad posterior del colegio 46 es más elevada que la del colegio 82. Esto se debe a que el tamaño de la muestra de la escuela 46 es de 21 estudiantes, mientras que la de la escuela 82 es de solo 5 estudiantes. Por lo tanto, el grado de certeza (cuantificable por medio del CV posterior, por ejemplo) acerca de \(\theta_{46}\) es mayor que el de \(\theta_{82}\).
Rankings basados en medias muestrales y medias posteriores no coinciden necesariamente.
Hay más evidencia de que \(\theta_{46}\) sea menor que \(\theta_{82}\).

6 Modelo con varianzas específicas

Caracterización dentro de los grupos: \[ y_{i, j}\mid\theta_j,\sigma_j^2 \stackrel{\text {iid}}{\sim} \textsf{N}\left(y_{i,j} \mid \theta_j,\sigma_j^2\right) \]
Caracterización entre los grupos: \[ \theta_{j}\mid\mu,\tau^2 \stackrel{\text {iid}}{\sim} \textsf{N}\left(\theta_j \mid \mu,\tau^2\right) \] \[ \sigma_j^2\mid\nu_0,\sigma^2_0 \stackrel{\text {iid}}{\sim}\textsf{GI}\left(\tfrac{\nu_0}{2},\tfrac{\nu_0\,\sigma^2_0}{2}\right) \]
Distribución previa: \[ p(\mu,\tau^2,\nu_0,\sigma^2_0) = p(\mu)\,p(\tau^2)\,p(\nu_0)\,p(\sigma^2_0) \] donde \[ \mu\sim\textsf{N}(\mu_0,\gamma_0^2)\qquad\tau^2\sim\textsf{GI}\left(\tfrac{\eta_0}{2},\tfrac{\eta_0\,\tau^2_0}{2}\right)\qquad p(\nu_0)\propto e^{-\alpha\nu_0}\qquad \sigma^2_0\sim\textsf{Gamma}(a,b) \]
Los parámetros del modelo son \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_m,\sigma^2_1,\ldots,\sigma^2_m,\mu,\tau^2,\nu_0,\sigma^2_0)\).
Los hiper-parámetros del model son \((\mu_0,\gamma_0^2,\eta_0,\tau^2_0, \alpha,a,b)\).

6.1 Estimación

Construir un muestreador de Gibbs para obtener muestras de la distribución posterior \(p(\boldsymbol{\theta}\mid\boldsymbol{y})\).

Distribuciones condicionales completas: \[ \begin{aligned} p\left(\theta_{j} \mid \text { resto }\right) &=\mathrm{N}\left(\theta_{j} \,\Big|\, \frac{n_{j} \bar{y}_{j} / \sigma_{j}^{2}+\mu / \tau^{2}}{n_{j} / \sigma_{j}^{2}+1 / \tau^{2}}, \frac{1}{n_{j} / \sigma_{j}^{2}+1 / \tau^{2}}\right) \\ p\left(\sigma_{j}^{2} \,\Big|\, \text { resto }\right) &=\operatorname{GI}\left(\sigma^{2} \mid \frac{\nu_{0}+n_{j}}{2}, \frac{\nu_{0} \sigma_{0}^{2}+\sum_{i=1}^{n_j}\left(y_{i, j}-\theta_{j}\right)^{2}}{2}\right) \\ p\left(\sigma_{0}^{2} \mid \text { resto }\right) &=\operatorname{G}\left(\sigma_{0}^{2} \,\Big|\, a+\frac{m \nu_{0}}{2}, b+\frac{\nu_{0}}{2} \sum_{j=1}^m \sigma_{j}^{-2}\right) \\ p\left(\nu_{0} \mid \text { resto }\right) & \propto\left[\frac{\left(\nu_{0} \sigma_{0}^{2} / 2\right)^{\nu_{0} / 2}}{\Gamma\left(\nu_{0} / 2\right)}\right]^{m}\left[\prod_{j=1}^m \sigma_j^{-2}\right]^{\nu_{0} / 2} {\exp}\left\{-\nu_{0}\left(\alpha+\frac{\sigma_{0}^{2}}{2} \sum_{j=1}^m \sigma_{j}^{-2}\right)\right\} \\ p(\mu \mid \text { resto }) &=\mathrm{N}\left(\mu \,\Big|\,\frac{m \bar{\theta} / \tau^{2}+\mu_{0} / \gamma_{0}^{2}}{m / \tau^{2}+1 / \gamma_{0}^{2}}, \frac{1}{m / \tau^{2}+1 / \gamma_{0}^{2}}\right) \\ p\left(\tau^{2} \mid \text { resto }\right) &=\operatorname{GI}\left(\tau^{2} \,\Big|\, \frac{\eta_{0}+m}{2}, \frac{\eta_{0} \tau_{0}^{2}+\sum_{j=1}^m\left(\theta_{j}-\mu\right)^{2}}{2}\right) \end{aligned} \]

Para muestrear de \(p(\nu_0 \mid \text{resto})\) se calcula esta distribución (en escala log) para un rango de valores de \(\nu_0\), se normaliza, y luego se muestra de la distribución discreta resultante.

Si \(\nu_0 = \infty\), entonces este modelo es igual al anterior donde todas las varianzas son iguales. De otro lado, si \(\nu_0 = 1\), entonces todas las varianzas serian muy diferentes y muy poco información acerca de las varianzas se compartiría a través de los grupos.

7 Ejemplo: Puntajes de matemáticas (cont.)

7.1 Estimación

######## hiper-parametros para obtener una previa no informativa
nu0  <- 1  ; s20 <- 100
eta0 <- 1  ; t20 <- 100
mu0  <- 50 ; g20 <- 25
a0   <- 1  ; b0  <-  1/100 ; wnu0 <- 1
########


######## valores iniciales
theta  <- ybar
sigma2 <- sv 
mu     <- mean(theta)
tau2   <- var(theta)
s20    <- 100
nu0    <- 1
########


######## setup
S      <- 50000
ncat   <- floor(S/10)
THETA  <- matrix(data=NA, nrow=S, ncol=m)
SIGMA2 <- matrix(data=NA, nrow=S, ncol=m)
MTSN   <- matrix(data=NA, nrow=S, ncol=4)  # mu, tau2, sigma02, nu0
nu0s   <- 1:100 # rango de valores para muestrear en p(nu_0 | rest)
########


######## MCMC
set.seed(1)
for(s in 1:S) {
        # actualizar thetas
        for(j in 1:m) {
                vtheta   <- 1/(n[j]/sigma2[j]+1/tau2)
                etheta   <- vtheta*(ybar[j]*n[j]/sigma2[j]+mu/tau2)
                theta[j] <- rnorm(1,etheta,sqrt(vtheta))
        }
        # actualizar sigma2s
        for(j in 1:m) { 
                nun       <- nu0+n[j]
                ss        <- nu0*s20+ sum((Y[[j]]-theta[j])^2)
                sigma2[j] <- 1/rgamma(1,nun/2,ss/2)
        }
        # actualizar s20
        s20 <- rgamma(1,a0+m*nu0/2,b0+nu0*sum(1/sigma2)/2)
        # actualizar nu0
        lpnu0 <- .5*nu0s*m*log(s20*nu0s/2)-m*lgamma(nu0s/2)+(nu0s/2-1)*sum(log(1/sigma2)) - nu0s*s20*sum(1/sigma2)/2 - wnu0*nu0s
        nu0   <- sample(x = nu0s, size = 1, prob=exp( lpnu0-max(lpnu0) ))
        # actualizar mu
        vmu <- 1/(m/tau2+1/g20)
        emu <- vmu*(m*mean(theta)/tau2 + mu0/g20)
        mu  <- rnorm(1,emu,sqrt(vmu))
        # actualizar tau2
        etam <-eta0+m
        ss   <- eta0*t20 + sum( (theta-mu)^2 )
        tau2 <-1/rgamma(1,etam/2,ss/2)
        # almacenar
        THETA[s,]  <- theta
        SIGMA2[s,] <- sigma2
        MTSN[s,]   <- c(mu,tau2,s20,nu0)
        ### progreso
        if (s%%ncat == 0) cat(100*round(s/S, 1), "% completed ... \n", sep = "" )
}

## 10% completed ... 
## 20% completed ... 
## 30% completed ... 
## 40% completed ... 
## 50% completed ... 
## 60% completed ... 
## 70% completed ... 
## 80% completed ... 
## 90% completed ... 
## 100% completed ...

######## FINAL DEL ALGORITMO

7.2 Inferencia

# comparar valores de la media posterior de sigma^2_j con las varianzas muestrales
sigma2.hat <- apply(SIGMA2,2,mean)
# grafico
par(mfrow=c(1,2),mar=c(3,3,1,1)+.3,mgp=c(1.75,.75,0))
#
plot(sv,sigma2.hat,xlab=expression(s^2),ylab=expression(hat( sigma^2)) )
abline(0,1)
#
plot(n, sv-sigma2.hat,xlab="sample size",ylab=expression(s^2-hat(sigma^2)))  
abline(h=0)

par(mfrow=c(2,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
#
plot(density(MST[,1],adj=2),lwd=2,main="",col="gray",xlab=expression(mu),
ylab=expression(paste(italic("p("),mu,"|",italic(y[1]),"...",italic(y[m]),")")))
lines(density(MTSN[,1],adj=2),lwd=2)
#
plot(density(MST[,3],adj=2),lwd=2,main="",col="gray",xlab=expression(tau^2), 
ylab=expression(paste(italic("p("),tau^2,"|",italic(y[1]),"...",italic(y[m]),")")))
lines(density(MTSN[,2],adj=2),lwd=2)
#
plot(table(MTSN[,4]),xlab=expression(nu[0]),
ylab=expression(paste(italic("p("),nu[0],"|",italic(y[1]),"...",italic(y[m]),")")))
#
plot(density(MTSN[,3],adj=2),lwd=2,main="",xlab=expression(sigma[0]^2),
ylab=expression(paste(italic("p("),sigma[0]^2,"|",italic(y[1]),"...",italic(y[m]),")")))

7.3 Ranking Bayesiano

# ranking bayesiano
ids  <- 1:m
that <- colMeans(THETA)
ic1  <- apply(X = THETA, MARGIN = 2, FUN = function(x) quantile(x, c(0.025,0.975)))
ic2  <- apply(X = THETA, MARGIN = 2, FUN = function(x) quantile(x, c(0.005,0.995)))

ranking <- order(that, decreasing = T) 
ids  <- ids [ranking]
that <- that[ranking]
ic1  <- ic1[,ranking]
ic2  <- ic2[,ranking]

k    <- 20
ids  <- ids [1:k]
that <- that[1:k]
ic1  <- ic1[,1:k]
ic2  <- ic2[,1:k]
colo <- c("blue","black")[as.numeric((ic2[1,] < 50) & (50 < ic2[2,]))+1]

# grafico
par(mfrow=c(1,1),mar=c(3,3,1.5,1),mgp=c(1.75,.75,0))
plot(NA, NA, xlab = "Colegio", ylab = "Puntaje", main = paste0("Ranking Top ", k, " Bayesiano"), 
     xlim = c(1,k), ylim = range(THETA), cex.axis = 0.75, xaxt = "n")
axis(side = 1, at = 1:k, labels = ids, cex.axis = 0.75)
abline(h = 50, col = "gray", lwd = 2)
for (j in 1:k) {
        segments(x0 = j, y0 = ic1[1,j], x1 = j, y1 = ic1[2,j], lwd = 3, col = colo[j])
        segments(x0 = j, y0 = ic2[1,j], x1 = j, y1 = ic2[2,j], lwd = 1, col = colo[j])
        lines(x = j, y = that[j], type = "p", pch = 16, cex = 0.8, col = colo[j])
}

7.4 Ranking Frecuentista

ids  <- 1:m
that <- ybar
ic1  <- lapply(X = Y, FUN = function(x) mean(x) + c(-1,1)*qnorm(0.975)*sd(x)/sqrt(length(x)))
ic2  <- lapply(X = Y, FUN = function(x) mean(x) + c(-1,1)*qnorm(0.995)*sd(x)/sqrt(length(x)))
ic1  <- matrix(unlist(ic1), 2, m) 
ic2  <- matrix(unlist(ic2), 2, m) 

ids  <- ids [ranking]
that <- that[ranking]
ic1  <- ic1[,ranking]
ic2  <- ic2[,ranking]

k    <- 20
ids  <- ids [1:k]
that <- that[1:k]
ic1  <- ic1[,1:k]
ic2  <- ic2[,1:k]
colo <- c("blue","black")[as.numeric((ic2[1,] < 50) & (50 < ic2[2,]))+1]

par(mfrow=c(1,1),mar=c(3,3,1.5,1),mgp=c(1.75,.75,0))
plot(NA, NA, xlab = "Colegio", ylab = "Puntaje", main = paste0("Ranking Top ", k, " Frecuentista"), 
     xlim = c(1,k), ylim = range(THETA), cex.axis = 0.75, xaxt = "n")
axis(side = 1, at = 1:k, labels = ids, cex.axis = 0.75)
abline(h = 50, col = "gray", lwd = 2)
for (j in 1:k) {
        segments(x0 = j, y0 = ic1[1,j], x1 = j, y1 = ic1[2,j], lwd = 3, col = colo[j])
        segments(x0 = j, y0 = ic2[1,j], x1 = j, y1 = ic2[2,j], lwd = 1, col = colo[j])
        lines(x = j, y = that[j], type = "p", pch = 16, cex = 0.8, col = colo[j])
}

Comparación de grupos

Juan Sosa PhD

jcsosam@unal.edu.co