10.10 Modelos jerárquicos

Los modelos jerárquicos involucran varios parámetros de tal manera que las creencias de unos de los parámetros dependen de manera significativa de los valores de otros parámetros. Por ejemplo, consideremos el caso en el que tenemos varias monedas acuñadas en la misma casa de monedas, es razonable pensar que una fábrica sesgada a águilas tenderá a producir monedas con sesgo hacia águilas. La estimación del sesgo de una moneda depende de la estimación del sesgo de la fábrica que a su vez está influido por los datos de todas las monedas. Veremos que la estructura de dependenica a lo largo de los parámetros generan estimaciones mejor informadas para todos los parámetros.

Si pensamos únicamente en dos monedas que provienen de la misma casa de moneda tenemos:

  1. Conocimientos iniciales de los posibles valores de los parámetros (sesgos de las monedas).

  2. Tenemos conocimiento inicial de la dependencia de los parámetros por provenir de la misma fábrica.

Cuando observamos lanzamientos de las monedas actualizamos nuestras creencias relativas a los sesgos de las monedas y también actualizamos nuestras creencias acerca de la dependencia de los sesgos.

Recordemos el caso de lanzamientos de una moneda, le asignamos una inicial beta, recordemos que la distribución beta tienen dos parámetors \(a\) y \(b\):

\[p(\theta)=\frac{1}{B(a,b)}\theta^{a-1}(1-\theta)^{b-1}\]

Con el fin de hacer los parámetros más intuitivos los podemos expresar en términos de la media \(\mu\) y el tamaño de muestra \(K\), donde \(\mu\) es la media de nuestro conocimiento inicial y la confianza está reflejada en el tamaño de muestra \(K\). Entonces los parámetros correspondientes en la distribución beta son:

\[a=\mu K, b = (1-\mu)K\]

Ahora introducimos el modelo jerárquico. En lugar de especificar un valor particular para \(\mu\), consideramos que \(\mu\) puede tomar distintos valores (entre \(0\) y \(1\)), y definimos una distribución de probabilidad sobre esos valores. Podemos pensar que esta distribución describe nuestra incertidumbre acerca de la construcción de la máquina que manufactura las monedas.

Veamos que en el caso de más de una moneda el modelo permite que cada moneda tenga un sesgo distinto pero ambas tenderán a tener un sesgo cercano a \(\mu\), algunas aleatoriamente tendrán un valor de \(\theta\) mayor a \(\mu\) y otras menor. Entre más grande \(K\) mayor será la consistencia de la acuñadora y los valores \(\theta\) serán más cercanos a \(\mu\). Si observamos varios lanzamientos de una moneda tendremos información tanto de \(\theta\) como de \(\mu\).

Para hacer un análisis bayesiano aún nos hace falta definir la distribución inicial sobre los parámetros \(\mu\), usemos una distribución Beta:

\[p(\mu)=beta(\mu|A_{\mu}, B_{\mu})\]

donde \(A_{\mu}\) y \(B_{\mu}\) se conocen como hiperparámetros y son constantes. En este caso, consideramos que \(\mu\) se ubica típicamente cerca de \(A_{\mu}/(A_{\mu} + B_{\mu})\) y \(K\) se considera constante.

Modelo jerárquico una moneda

Recordemos que en el ejemplo de una moneda teníamos que la verosimilitud era Bernoulli:

\[p(x|theta) = \theta^x(1-\theta)^{1-x}\]

Y si utilizamos las iniciales Beta para \(\mu\) y \(\theta\) como discutimos arriba, solo nos resta aplicar la regla de Bayes con nuestros dos parámetros desconocidos \(\mu\) y \(\theta\):

\[p(\theta, \mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\]

Hay dos aspectos a considerar en el problema:

  1. La verosimilitud no depende de \(\mu\) por lo que \[p(x|\theta, \mu)=p(x|\theta)\]

  2. La distribución inicial en el espacio de parámetros bivariado se puede factorizar:

\[p(\theta,\mu)=p(\theta|\mu)p(\mu)\]

Por lo tanto

\[p(\theta,\mu|x)=\frac{p(x|\theta,\mu)p(\theta,\mu)}{p(x)}\] \[=\frac{p(x|\theta)p(\theta|\mu)p(\mu)}{p(x)}\]

¿Cuál es el modelo gráfico?

Aproximación por cuadrícula

En el caso jerárquico, no se puede derivar la distribución posterior de manera analítica pero si los parámetros e hiperparámetros toman un número finito de
valores y no hay muchos de ellos, podemos aproximar la posterior usando aproximación por cuadrícula, utilizar este enfoque en el ejemplo introductorio ayuda a visualizar el proceso.

A continuación graficamos las distribuciones correspondientes al caso en que la distribución del hiperparámetro \(\mu\) tiene la forma de una distribución \(Beta(2, 2)\), es decir creemos que la media de la acuñadora \(\mu\) es \(0.5\), pero existe bastante incertidumbre acerca del valor.

La distribución de \(\theta\), esto es la distribución inicial que refleja la dependencia entre \(\theta\) y \(\mu\) se expresa por medio de otra distribución beta, en el ejemplo usamos \(K=100\):

\[\theta|\mu \sim beta(\mu 100, (1-\mu)100)\]

Esta inicial expresa un alto grado de certeza que una acuñadora con hiperparámetro \(\mu\) genera monedas con sesgo cercano a \(\mu\)

Verosimilitud.

Posterior.

Código para el análisis y gráficas de arriba.

A_mu <- 2; B_mu <- 2; K <- 100
# Distribución inicial conjunta p(theta, mu)
p_conjunta <- function(mu, theta, A_mu, B_mu, K){
    # marginal p(mu)
    p_mu <- dbeta(mu, A_mu, B_mu)
    # condicional p(theta | mu)
    p_theta_mu <- dbeta(theta, mu * K, (1 - mu) * K)
    p_mu * p_theta_mu
}

grid <- expand.grid(theta = seq(0.01, 0.99, 0.005), 
    mu = seq(0.01, 0.99, 0.005))

grid_inicial <- grid %>% 
    mutate(p_inicial = p_conjunta(mu, theta, A_mu, B_mu, K))

plot_conj <- ggplot(grid_inicial, aes(x = theta, y = mu, z = p_inicial)) + 
    stat_contour(binwidth = 2, aes(color = ..level..)) +
    scale_x_continuous(expression(theta), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 

grid_mu <- grid_inicial %>%
    group_by(mu) %>%
    summarise(p_mu = sum(p_inicial)) %>%
    ungroup() %>%
    mutate(p_mu = p_mu / sum(p_mu))

plot_mu <- ggplot(grid_mu, aes(x = mu, y = p_mu)) +
    geom_path() +
    labs(y = expression(p(mu)), x = expression(mu)) + 
    ylim(0, 0.015) + 
    coord_flip()

mus <- rbeta(5000, A_mu, B_mu)
sims_marg <- data.frame(sims = 1:5000, 
    sims_marg = rbeta(5000, mus * 100, (1 - mus) * 100))

plot_theta <- ggplot(sims_marg, aes(x = sims_marg)) + 
    geom_density(adjust = 2) +
    labs(y = expression(p(theta)), x = expression(theta))

# p(theta|mu=0.75)
sims_cond_1 <- rbeta(5000, 0.75 * 100, (1 - 0.75) * 100)
# p(theta|mu=0.25)
sims_cond_2 <- rbeta(5000, 0.25 * 100, (1 - 0.25) * 100)

# marginal = sims_marg

sims <- data.frame(sim = 1:5000, 
    cond_1 = sims_cond_1, 
    cond_2 = sims_cond_2) %>%
    gather(dist, values, -sim) %>%
    mutate(mu = ifelse(dist == "cond_1", 0.75, 0.25))


plots_marginales <- ggplot(sims, aes(x = values)) +
    labs(y = expression(p(paste(theta, l, mu))), x = expression(theta)) + 
    geom_density(adjust = 2) + 
    facet_wrap(~ mu, ncol = 1)

grid.arrange(plot_conj, plot_mu, plot_theta, plots_marginales, ncol = 2)
likeBern <- function(z, N){
    function(theta){
      theta ^ z * (1 - theta) ^ (N - z)
    }
}

# Valores observados
z <- 9; N <- 12
mi_like <- likeBern(z, N)

grid_like <- expand.grid(x = seq(0.01, 1, 0.05), y = seq(0.01, 1, 0.05))
grid_like <- grid_like %>% 
    mutate(z = mi_like(x))

ggplot(grid_like, aes(x = x, y = y, z = z)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(L(theta,mu)), guide = FALSE) +
    coord_fixed()

# cuadrícula
grid_post <- grid_inicial %>%
    mutate(
      prior = p_inicial / sum(p_inicial),
      Like = theta ^ z * (1 - theta) ^ (N - z), # verosimilitud 
      posterior = (Like * prior) / sum(Like * prior)
    )  

post_conj <- ggplot(grid_post, aes(x = theta, y = mu, z = posterior)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 
grid_mu <- grid_post %>%
    group_by(mu) %>%
    summarise(post_mu = sum(posterior)) %>%
    ungroup() %>%
    mutate(post_mu = post_mu / sum(post_mu))

post_mu <- ggplot(grid_mu, aes(x = mu, y = post_mu)) +
    geom_path() +
    labs(y = expression(p(paste(mu, l, x))), x = expression(mu)) + 
    coord_flip()

grid_theta <- grid_post %>%
    group_by(theta) %>%
    summarise(post_theta = sum(posterior)) %>%
    ungroup() %>%
    mutate(post_theta = post_theta / sum(post_theta))

post_theta <- ggplot(grid_theta, aes(x = theta, y = post_theta)) +
    geom_path() +
    labs(y = expression(p(paste(theta, l, x))), x = expression(theta))

grid_theta_m <- grid_post %>%
    filter(mu == 0.75 | mu == 0.25) %>%
    group_by(theta, mu) %>%
    summarise(post_theta = sum(posterior)) %>%
    group_by(mu) %>%
    mutate(post_theta = post_theta / sum(post_theta))

post_marginales <- ggplot(grid_theta_m, aes(x = theta, y = post_theta)) +
    geom_path() + 
    facet_wrap(~ mu, ncol = 1) +
    labs(y = expression(p(paste(theta, l, x, mu))), 
         x = expression(theta))

Multiples monedas de una misma fábrica

La sección anterior considera el escenario en que lanzamos una moneda y hacemos inferencia del parámetro de sesgo \(\theta\) y del hiperparámetro \(\mu\). Ahora consideramos recolectar información de múltiples monedas, si cada moneda tiene su propio sesgo \(\theta_j\) entonces debemos estimar un parámetro distinto para cada moneda.

Suponemos que todas las monedas provienen de la misma fábrica, esto implica que tenemos la misma información inicial \(\mu\) para todas las monedas. Suponemos también que cada moneda se acuña de manera independiente, esto es, que condicional al parámetro \(\mu\) los parámetros \(\theta_j\) son independientes en nuestros conocimientos iniciales.

Posterior vía aproximación por cuadrícula

Supongamosque tenemos dos monedas de la misma fábrica. El objetivo es estimar los sesgos \(\theta_1\), \(\theta_2\) de las dos monedas y estimar simultáneamente el parámetro \(\mu\) correspondiente a la casa de moneda que las fabricó.

Inicial.

Verosimilitud.

Posterior.

Y el código para el análisis y gráficas de arriba.

A_mu <- 2; B_mu <- 2; K <- 1000

# Distribución inicial p(theta, mu)
p_conjunta <- function(mu, theta, A_mu, B_mu, K){
    # marginal p(mu)
    p_mu <- dbeta(mu, A_mu, B_mu)
    # condicional p(theta | mu)
    p_theta_mu <- dbeta(theta, mu * K, (1 - mu) * K)
    p_mu * p_theta_mu
}

grid_inicial <- expand.grid(theta = seq(0.01, 0.99, 0.005),
    mu = seq(0.01, 0.99, 0.005))

grid_inicial <- grid_inicial %>% 
    mutate(p_inicial = p_conjunta(mu, theta, A_mu, B_mu, K))

plot_conj_1 <- ggplot(grid_inicial, aes(x = theta, y = mu, z = p_inicial)) + 
    stat_contour(binwidth = 1, aes(color = ..level..)) +
    scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 

plot_conj_2 <- ggplot(grid_inicial, aes(x = theta, y = mu, z = p_inicial)) + 
    stat_contour(binwidth = 1, aes(color = ..level..)) +
    scale_x_continuous(expression(theta[2]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta,mu)), guide = FALSE) 

grid_mu <- grid_inicial %>%
    group_by(mu) %>%
    summarise(p_mu = sum(p_inicial)) %>%
    ungroup() %>%
    mutate(p_mu = p_mu / sum(p_mu))

plot_mu <- ggplot(grid_mu, aes(x = mu, y = p_mu)) +
    geom_path() +
    labs(y = expression(p(paste(mu))), x = expression(mu)) + 
    ylim(0, 0.015) + coord_flip()

p_marg <- grid_inicial %>%
    group_by(theta) %>%
    summarise(marg = sum(p_inicial)) %>%
    ungroup() %>%
    mutate(marg = marg / sum(marg))

plot_marg_1 <- ggplot(p_marg, aes(x = theta, y = marg)) + 
    geom_path() +
    labs(y = expression(p(theta[1])), x = expression(theta[1])) +
    ylim(0, 0.02)

plot_marg_2 <- ggplot(p_marg, aes(x = theta, y = marg)) + 
    geom_path() +
    labs(y = expression(p(theta[2])), x = expression(theta[2])) +
    ylim(0, 0.02)

grid.arrange(plot_conj_1, plot_conj_2, plot_mu, plot_marg_1, plot_marg_2, 
    ncol = 3)
z_1 <- 3; N_1 <- 15
z_2 <- 4; N_2 <- 5

mi_like_1 <- likeBern(z_1, N_1)
mi_like_2 <- likeBern(z_2, N_2)

grid_like <- grid_like %>% 
    mutate(L_1 = mi_like_1(x), L_2 = mi_like_2(x))

plot_like_1 <- ggplot(grid_like, aes(x = x, y = y, z = L_1)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

plot_like_2 <- ggplot(grid_like, aes(x = x, y = y, z = L_2)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(L(theta,mu)), guide = FALSE)

grid.arrange(plot_like_1, plot_like_2, ncol = 3)
grid_post <- grid_inicial %>%
    mutate(
        prior = p_inicial / sum(p_inicial),
        Like_1 = theta ^ z_1 * (1 - theta) ^ (N_1 - z_1), # verosimilitud 
        posterior_1 = (Like_1 * prior) / sum(Like_1 * prior),
        Like_2 = theta ^ z_2 * (1 - theta) ^ (N_2 - z_2),
        posterior_2 = (Like_2 * prior) / sum(Like_2 * prior)
      )  

post_conj_1 <- ggplot(grid_post, aes(x = theta, y = mu, z = posterior_1)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta[1]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta[1],mu)), guide = FALSE) 

post_conj_2 <- ggplot(grid_post, aes(x = theta, y = mu, z = posterior_2)) + 
    stat_contour(aes(color = ..level..)) +
    scale_x_continuous(expression(theta[2]), limits = c(0, 1)) +
    scale_y_continuous(expression(mu), limits = c(0, 1)) +
    scale_color_gradient(expression(p(theta[2],mu)), guide = FALSE) 

grid_mu <- grid_post %>%
    group_by(mu) %>%
    summarise(post_mu = sum(posterior_1) * sum(posterior_2)) %>%
    ungroup() %>%
    mutate(post_mu = post_mu / sum(post_mu))

post_mu <- ggplot(grid_mu, aes(x = mu, y = post_mu)) +
    geom_path() +
    labs(y = expression(p(paste(mu, l, x))), x = expression(mu)) + 
    ylim(0, 0.015) + 
    coord_flip()

grid_theta <- grid_post %>%
    group_by(theta) %>%
    summarise(post_theta_1 = sum(posterior_1), 
        post_theta_2 = sum(posterior_2)) %>%
    ungroup() %>%
    mutate(post_theta_1 = post_theta_1 / sum(post_theta_1))

post_theta_1 <- ggplot(grid_theta, aes(x = theta, y = post_theta_1)) +
    geom_path() +
    labs(y = expression(p(paste(theta[1], l, x))), x = expression(theta[1])) + 
    ylim(0, 0.020)

post_theta_2 <- ggplot(grid_theta, aes(x = theta, y = post_theta_2)) +
    geom_path() +
    labs(y = expression(p(paste(theta[2], l, x))), x = expression(theta[2])) + 
    ylim(0, 0.020)

grid.arrange(post_conj_1, post_conj_2, post_mu, post_theta_1, post_theta_2, 
             ncol = 3)

Stan

La sección anterior utilizó un modelo simplificado con el objetivo de poder visualizar el espacio de parámetros. Ahora incluiremos más parámetros para hacer el problema más realista. En los ejemplos anteriores fijamos el grado de dependencia entre \(\mu\) y \(\theta\) de manera arbitraria, a través de \(K\), de tal manera que si \(K\) era grande, los valores \(\theta_j\) individuales se situaban cerca de \(\mu\) y cuando \(K\) era pequeña se permitía más variación.

En situaciones reales no conocemos el valor de \(K\) por adelantado, por lo que dejamos que los datos nos informen acerca de valores creíbles para \(K\). Intuitivamente, cuando la proporción de águilas observadas es similar a lo largo de las monedas, tenemos evidencia de que \(K\) es grande, mientras que si estas proporciones difieren mucho, tenemos evidencia de que \(K\) es pequeña. Debido a que \(K\) pasará de ser una constante a ser un parámetro lo llamaremos \(\kappa\).

La distribución inicial para \(\kappa\) será una Gamma.

\(\kappa \sim Gamma(1, 0.1)\), tiene media \(10\) y desviación estándar \(10\).

Los datos consisten en \(5\) monedas, cada una se lanza \(5\) veces, resultando \(4\) de ellas en \(1\) águila y \(4\) soles y otra en \(3\) águilas y \(2\) soles.

Graficamos la distribución posterior de \(\mu\) y \(\theta_i\).

Y para el resto de los parámetros.

Ejemplo: estimación de tasas de mortalidad

En esta sección veremos un problema de estimación de tasas de mortalidad, la referencia de este ejemplo es Albert (2009). Plantearemos \(3\) alternativas de modelación para resolver el problema: 1) modelo de unidades iguales, 2) modelo de unidades independientes, y 3) modelo jerárquico.

Los datos consisten en información de todas las cirugías de transplante de corazón llevadas a cabo en Estados Unidos en un periodo de \(24\) meses, entre octubre de \(1987\) y diciembre de \(1989\). Para cada uno de los \(131\) hospitales, se registró:

  • el número de cirugías de transplante de corazón,

  • y el número de muertes durante los \(30\) días posteriores a la cirugía, \(y\).

  • Además, se cuenta con una predicción de la probabilidad de muerte de cada paciente individual. Esta predicción esta basada en un modelo logístico que incluye información a nivel paciente como condición médica antes de la cirugía, género, sexo y raza. En cada hospital se suman las probabilidades de muerte de sus pacientes para calcular el número esperado de muertes, \(e\), conocido como la exposición del hospital. \(e\) refleja el riesgo de muerte debido a la mezcla de pacientes que componen un hospital particular.

La tabla de datos que analizaremos considera únicamente los \(94\) hospitales que llevaron a cabo \(10\) o más cirugías, y contiene las duplas \((y_{j}, e_{j})\) que corresponden al número observado de muertes y número de expuestos para el \(j\)-ésimo hospital, con \(j = 1,...,94\).

El objetivo es obtener buenas estimaciones de las tasas de mortalidad de cada hospital. Antes de comenzar a ajustar modelos complejos vale la pena observar los datos.

  • El cociente \(\{y_{j}/e_{j}\}\) es el número observado de muertes por unidad de exposición y se puede ver como una estimación de la tasa de mortalidad.

  • La siguiente figura grafica, en el eje vertical los cocientes \(\{y_{j}/e_{j}\}\) multiplicados por \(1000\) (con la intención de que la tasa indique número de muertes por \(1000\) expuestos), y en el eje horizontal el número de expuestos \(\{e_{j}\}\) -en escala logarítmica- para los \(94\) hospitales. Cada punto representa un hospital y esta etiquetado con el número de muertes observadas \(\{y_{j}\}\).

  • En la gráfica podemos notar que las tasas estimadas son muy variables, especialmente para programas con baja exposición.

  • Observemos también, que la mayoría de los programas que no experimentaron muertes tienen bajo número de expuestos.

La variabilidad de las tasas y los hospitales sin muertes sugieren un problema de tamaño de muestra:

Consideremos un hospital con \(700\) expuestos, la muerte por transplante de corazón no es común, el promedio nacional es de \(9\) por cada \(10,000\) expuestos, por lo que con \(700\) expuestos es probable que no se presente ninguna muerte,

  • en este caso el hospital pertenece al \(10\%\) de los hospitales con menor tasa de mortalidad.

  • Sin embargo, existe la posibilidad de que se observe una muerte, con lo cual el hospital tendrá una tasa de mortalidad que es lo suficientemente alta para que quede en el \(25\%\) de los hospitales con mayor tasa de mortalidad observada.

Una vez reconocido el problema de utilizar los datos crudos para estimar las tasas de mortalidad plantearemos \(3\) alternativas de modelación:

  1. Unidades iguales. Considera que los estudios son replicas idénticas una de otra, en este sentido vemos a las observaciones de todos los estudios como muestras independientes de una misma población y consecuentemente todas las \(\lambda_{j}\) se consideran iguales,

\[ \begin{eqnarray} \nonumber y_{j} &\sim& f(y_{j}|\lambda),\\ \nonumber \lambda &\sim& g(\lambda). \end{eqnarray} \]

El modelo gráfico correspondiente a este enfoque (omitiendo la distribución de \(\lambda\)) es:

  1. Unidades independientes. El extremo opuesto a unidades iguales considera que los estudios son tan diferentes que los resultados de cada uno no proveen información acerca de los resultados de ningún otro y por tanto realizamos estimaciones individuales para cada \(\lambda_{j}\),

\[ \begin{align} \nonumber y_{j} &\sim f(y_{j}|\lambda_{j}),\\ \nonumber \lambda_j &\stackrel{\mathrm{iid}}{\sim} g(\lambda_j). \end{align} \]

Con el siguiente modelo gráfico:

  1. Jerárquico. Es un análisis intermedio que considera los parámetros de interés \(\lambda_{j}\) como provenientes de una distribución común,

\[ \begin{eqnarray} \nonumber y_{j} &\sim& f(y_{j}|\lambda_{j}),\\ \nonumber \lambda_{j} &\sim& g(\lambda_j|\theta),\\ \theta &\sim& h(\theta). \end{eqnarray} \]

De esta manera, la probabilidad conjunta del modelo refleja una dependencia entre los parámetros al mismo tiempo que permite variaciones entre los estudios. Como resultado se estima una \(\lambda_{j}\) diferente para cada estudio usando información de todos.

Modelo de unidades iguales

Supongamos que las tasas de mortalidad son iguales a lo largo de los hospitales. Estimaremos la tasa de mortalidad con el modelo,

\[ \begin{eqnarray} \nonumber y_{j} \sim Poisson(e_{j}\lambda), \end{eqnarray} \]

donde \(y_{j}\) es el número de muertes observadas por transplante corazón en el hospital \(j\), \(e_{j}\) es el número de expuestos, y \(\lambda\) es la tasa de mortalidad, medida en número de muertes por unidad de exposción y común para todos los hospitales.

Debido a que no contamos con información inicial acerca de la tasa de mortalidad, asignamos a \(\lambda\) una distribución inicial no informativa, de la forma

\[ \begin{eqnarray} \nonumber g(\lambda)\propto\frac{1}{\lambda}. \end{eqnarray} \]

Sea \(y=(y_1,...,y_{94})\), usamos el Teorema de Bayes para calcular la densidad posterior de \(\lambda\),

\[ \begin{eqnarray} \nonumber g(\lambda|y) &\propto& g(\lambda)f(y|\lambda) \\ \nonumber &=& g(\lambda)\prod_{i=1}^{94}f(y_{i}|\lambda)\\ \nonumber &=& \frac{1}{\lambda}\prod_{i=1}^{94}\bigg(\frac{exp(-e_{i}\lambda)(e_{i}\lambda)^{y_{i}}}{y_{i}!}\bigg)\\ \nonumber &\propto& \lambda^{(\sum_{i=1}^{94}y_{j}-1)}exp\bigg(-\sum_{i=1}^{94}e_{j}\lambda\bigg), \end{eqnarray} \]

identificamos la densidad posterior como una \(Gamma(\sum_{i=1}^{94}y_{i},\sum_{i=1}^{94}e_{i})\), donde expresamos la función de densidad de una distribución \(Gamma(\alpha, \beta)\) usando el parámetro de forma \(\alpha\) y el inverso del parámetro de escala \(\beta\) de manera que la función de densidad es,

\[ \begin{align} \nonumber &f(x|\alpha, \beta) = \beta^{\alpha}\frac{1}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x}I_{(0,\infty)}(x)\\ \nonumber &\mbox{para }x \geq 0 \mbox{ y }\alpha \mbox{, }\beta > 0. \end{align} \]

Para verificar el ajuste del modelo utilizaremos la distribución predictiva posterior. Denotemos \(y_{j}^*\) el número de muertes por transplante en el hospital \(j\) con exposición \(e_{j}\) en una muestra futura. Condicional a \(\lambda\), \(y_{j}^*\) se distribuye \(Poisson(e_{j}\lambda)\), no conocemos el verdadero valor de \(\lambda\), sin embargo nuestro conocimiento actual está contenido en la densidad posterior \(g(\lambda|y)\). Por tanto, la distribución predictiva posterior de \(y_{j}^*\) esta dada por:

\[ \begin{eqnarray} \nonumber f(y_{j}^*|e_{j},y) = \int f(y_{j}^*|e_{j}\lambda)g(\lambda|y)d\lambda, \end{eqnarray} \]

donde \(f(y_j^*|e_{j}\lambda)\) es una densidad Poisson con media \(\lambda\). La densidad predictiva posterior representa la verosimilitud de observaciones futuras basadas en el modelo ajustado.

En este ejemplo, la densidad \(f(y_{j}^*|e_{j},y)\) representa el número de transplantes que se predecirían para un hospital con exposición \(e_{j}\). Entonces, si el número observado de muertes \(y_{j}\) no está en las colas de la distribución predictiva, diríamos que la observación es consistente con el modelo observado.

La figura muestra los histogramas, obtenidos con simulación, de la distribución predictiva posterior de \(10\) hospitales (se tomó una muestra de los \(94\)). Para cada hospital, se escribió el número de expuestos, \(e\), y se grafica una línea vertical en el número de muertes observado. Notemos que en muchos de los histogramas el número de muertes observado se encuentra en las colas de las distribuciones, lo que sugiere que nuestras observaciones son inconsistentes con el modelo ajustado.

Ahora examinaremos la consistencia de las muertes observadas, \(y_{j}\), para todos los hospitales. Para cada distribución predictiva posterior calculamos la probabilidad de que una observación futura \(y_{j}^*\) sea al menos tan extrema como \(y_{j}\), estas probabilidades son comúnmente llamadas valores \(p\) predictivos:

\[ \begin{eqnarray} \nonumber P(extremos) = min\{P(y_{j}^* \leq y_{j}),P(y_{j}^*\geq y_{j})\} \end{eqnarray} \]

En la figura anterior graficamos las probabilidades de extremos (calculadas con simulación) en el eje vertical, contra el número de exposición en escala logarítmica de cada hospital en el eje horizontal. Notemos que muchas de estas probabilidades son pequeñas, el \(28\%\) son menores a \(0.15\), lo que indica que para el \(28\%\) de los hospitales el número de muertes observado \(y_{j}\) está en la cola de la distribución y por consiguiente el modelo es inadecuado.

Modelo de unidades independientes

Consideremos ahora que los hospitales son independientes, estimaremos la tasa de mortalidad para cada hospital con el modelo

\[ \begin{eqnarray} \nonumber y_{j}\sim Poisson(e_{j}\lambda_{j}), \end{eqnarray} \]

donde \(y_{j}\) es el número de muertes observadas por transplante corazón en el hospital \(j\), \(e_{j}\) es el número de expuestos, y \(\lambda_{j}\) es la tasa de mortalidad, medida en número de muertes por unidad de exposición, con \(j=1,...,94\). Utilizamos el subíndice \(j\) para enfatizar que son parámetros diferentes, cada uno estimado con sus propios datos.

Por facilidad asignamos a \(\lambda_{j}\) una distribución inicial \(Gamma(\alpha_{0}, \beta_{0})\) que es conjugada de la distribución Poisson,

\[ \begin{eqnarray} \nonumber g(\lambda_{j}|\alpha_{0},\beta_{0}) = \frac{\beta_{0}^{\alpha_{0}} \lambda_{j}^{\alpha_{0}-1}exp(-\lambda_{j} \beta_{0})}{\Gamma(\alpha_{0})}, \lambda_{j}>0. \end{eqnarray} \]

Calculamos la densidad posterior de \(\lambda_{j}\) usando el Teorema de Bayes,

\[ \begin{eqnarray} \nonumber g(\lambda_{j}|y_{j},\alpha_{0},\beta_{0}) &\propto& g(\lambda_{j}|\alpha_{0},\beta_{0})f(y_{j}|\lambda_{j})\\ \nonumber &=& \frac{\beta_{0}^{\alpha_{0}} \lambda_{j}^{\alpha_{0}-1}exp(-\lambda_{j} \beta_{0})}{\Gamma(\alpha_{0})}\frac{exp(-e_{j}\lambda_{j})(e_{j}\lambda_{j})^{y_{j}}}{y_{j}!}\\ \nonumber &\propto& \lambda_{j}^{\alpha_{0}+y_{j}-1}exp(-\lambda_{j}(\beta_{0}+e_{j})), \end{eqnarray} \]

y obtenemos que \(\lambda_{j}|y_{j} \sim Gamma(\alpha_{0}+y_{j}, \beta_{0}+e_{j})\), con media

\[ \begin{eqnarray} \nonumber E(\lambda_{j}|y_{j},\alpha_{0},\beta_{0}) &=& \frac{\alpha_{0}+y_{j}}{\beta_{0}+e_{j}}\\ \label{eqn:pond.indep} &=& (1-A_{j})\frac{y_{j}}{e_{j}}+A_{j}\frac{\alpha_{0}}{\beta_{0}} \end{eqnarray} \]

donde

\[ \begin{eqnarray} A_{j}=\frac{\beta_{0}}{\beta_{0}+e_{j}}, \end{eqnarray} \] y varianza

\[ \begin{eqnarray} \nonumber Var(\lambda_{j} |y_{j},\alpha_{0},\beta_{0}) = \frac{\alpha_{0}+y_{j}}{(\beta_{0}+e_{j})^2}. \end{eqnarray} \]

Notemos que podemos escribir la media posterior como un promedio ponderado de la tasa observada \(y_{j}/e_{j}\) y la media inicial \(\alpha_{0} / \beta_{0}\). El factor \(A_{j}\) es el encogimiento hacia la media inicial y depende del número de exposición de cada hospital y del parámetro de escala \(\beta_{0}\).

Efecto de \({\beta_{0}}\) y \({e_{j}}\) en la distribución posterior

Consideremos la media y varianza posteriores de \(\lambda_{j}\) para un hospital particular. Tomando \(\alpha_{0}\) fija, valores mayores de \(\beta_{0}\) corresponden a una menor varianza en la distribución inicial pues \(Var(\lambda_{j} | \alpha_{0},\beta_{0}) = \alpha_{0}/\beta_{0}^2\). La varianza de la distribución inicial se refleja en la media de la distribución posterior de \(\lambda_{j}\) a través de \(\beta_{0}\), menor varianza (mayor \(\beta_{0}\)) corresponde a mayor encogimineto hacia la media inicial. Este efecto de \(\beta_{0}\) concuerda con la incertidumbre de nuestro conocimiento, pues menor varianza en la distribución inicial indica menor incertidumbre y la aportación de la media inicial a la media posterior es más importante. La elección de \(\beta_{0}\) también afecta la varianza, pues un menor valor de \(\beta_{0}\) implica una mayor varianza tanto en la distribución inicial como en la final.

Por su parte, tomando \(\alpha_{0}\) y \(\beta_{0}\) fijas, el efecto del número de expuestos \(e_{j}\) sobre la media posterior de \(\lambda_{j}\) actúa de manera contraria a \(\beta_{0}\). Un mayor número de expuestos tiene como consecuencia un menor encogimiento hacia la media inicial, dando mayor importancia a la tasa observada \(y_{j}/e_{j}\). Esto refleja que más expuestos implican más información proveniente de la muestra, restando importancia a nuestro conocimiento inicial. En cuanto a la varianza posterior, es inversamente proporcional al número de expuestos indicando nuevamente que más expuestos implica más conocimiento y por tanto menor incertidumbre.

En la siguiente fugura mostraremos el encogimiento hacia la media inicial bajo distintos escenarios de \(\beta_{0}\), cada punto representa un hospital y el color indica a que valor de \(\beta_{0}\) corresponde. En esta gráfica podemos apreciar el mayor encogimiento para valores mayores de \(\beta_{0}\) y el decaimiento en el encogimiento conforme aumenta el número de expuestos.

Distribución inicial

Establecimos que por facilidad se utilizará una inicial \(Gamma(\alpha_{0},\beta_{0})\), sin embargo hace falta asignar valores a los parámetros iniciales. Analizaremos \(3\) distintas combinaciones de parámetros iniciales y su efecto en las estimaciones posteriores.

Para cada combinación de parámetros la siguiente tabla muestra los deciles de \(2000\) simulaciones de una distribución \(Gamma(\alpha_0,\beta_{0})\) y los deciles de las muertes observadas por cada \(1000\) expuestos considerando todos los hospitales; multiplicamos las simulaciones por \(1000\) para que indiquen tasas de mortalidad medidas en número de muertes por \(1000\) expuestos. El propósito de la tabla es describir la forma de la distribución \(Gamma\) al cambiar los parámetros, por ejemplo, una \(Gamma(0.01,0.01)\) es muy plana y podríamos describirla como poco informativa, mientras que una \(Gamma(1,1000)\) tiene una forma más cercana a las tasas de mortalidad observadas.

decil \((0.01,0.01)\) \((0.5,0.01)\) \((1,1000)\) Observados
min 0 0 0 0
10 0 890.3 0.1 0
20 0 3134.7 0.2 0.3
30 0 7769.0 0.3 0.5
40 0 15036.0 0.5 0.7
50 0 24293.9 0.7 0.8
60 0 37306.0 0.9 1.1
70 0 57640.2 1.2 1.4
80 0 84255.5 1.6 1.7
90 4.2 140721.9 2.3 2.0
max 287307.6 592073.2 7.0 3.9

Ahora realizamos una gráfica para cada pareja de parámetros
\((\alpha_0,\beta_{0})\) donde mostramos en negro las tasas observadas, y en color las estimaciones posteriores de las tasas de mortalidad con intervalos del \(95\%\) de probabilidad, ambas medidas en número de muertes por cada \(1000\) expuestos. Cada punto representa un hospital y el color corresponde al número de muertes observadas, \(\{y_{j}\}\). Regresaremos a esta gráfica más adelante pero por lo pronto observemos que tanto las estimaciones como los intervalos de probabilidad son muy diferentes al cambiar los parámetros de la distribución inicial.

Ahora justificaremos la elección de las parejas de parámetros iniciales. La única información con la que contamos para definir una distribución inicial es que la tasa de mortalidad por transplante de corazón debe ser positiva y no muy grande. Debido a que no tenemos más información nuestro primer modelo utiliza una inicial vaga.

No informativas

  • Asignamos a \(\lambda_{j}\) una distribución inicial \(Gamma(0.01, 0.01)\). La tabla de deciles indica que es una distribución muy plana y si observamos la gráfica de encogimiento notamos que para una inicial con este valor en el parámetro de escala \(\beta_{0}\), el encogimiento de la media posterior hacia la media inicial es muy cercano a cero, dando poca importancia a la media inicial \(\alpha_{0}/\beta_{0}=1\). Como consecuencia, las estimaciones posteriores de \(\lambda_{j}\) son casi iguales a las tasas observadas \(\{y_{j}/e_{j}\}\), y se presentan los problemas de tamaño de muestra notados al usar las tasas observadas como estimaciones de las tasas de mortalidad. Además, los intervalos de confianza para los hospitales que no experimentaron muertes son poco creíbles pues son mucho más chicos que el resto.

  • Consideramos una inicial \(Gamma(0.5,0.01)\). En este caso, la elección de los parámetros se basó en la distribución inicial no informativa de Jeffreys, consiste en una \(Gamma(0.5,0)\). Es una distribución impropia por lo que cambiamos el parámetro de escala por \(0.01\) para obtener una inicial propia con varianza grande. Los resultados no son muy razonables, pues la media inicial es \(50\), mucho mayor a las tasas observada para todos los hospitales, provocando que las estimaciones del modelo sean mayores a las tasas observadas en todos los hospitales. Este efecto es contrario al que buscábamos al usar una inicial vaga pues la distribución inicial tiene un impacto muy fuerte en las estimaciones posteriores.

Informativa
  • Consideramos una distribución inicial \(Gamma(1,1000)\), ésta inicial es informativa. Para obtener sus parámetros igualamos los media y varianza teórica de la distribución \(Gamma\) con la media y varianza observadas en el conjunto de tasas de mortalidad tomando en cuenta todos los hospitales. Las estimaciones de las tasas de mortalidad que obtenemos parecen razonables, sin embargo, especificar la distribución inicial con la muestra tiene problemas lógicos y prácticos:
    • Los datos se están usando \(2\) veces, primero la información de todos los hospitales se usa para especificar la distribución inicial, y después la información de cada hospital se usa para estimar su tasa de mortalidad \(\lambda_{j}\), lo que ocasiona que sobreestimemos nuestra precisión.
    • De acuerdo a la lógica bayesiana no tiene sentido estimar \(\alpha_{0}\) y \(\beta_{0}\), pues son parte de la distribución inicial y no deben depender de los datos.

A pesar de los problemas señalados parece ser conveniente intentar mejorar las estimaciones individuales de \(\lambda_{j}\) usando la información de todos los hospitales. La manera correcta de hacerlo es establecer un modelo de probabilidad para todo el conjunto de parámetros
\(\{\alpha,\beta,\lambda_{1},...,\lambda_{94}\}\) y después realizar un análisis de la distribución conjunta. Se llevará a cabo un análisis completo en la siguiente sección en donde se usará un modelo jerárquico.

Podemos concluir que el modelo de unidades independientes no es robusto para nuestros datos pues las estimaciones posteriores de las tasas de mortalidad son muy sensibles a la elección de los parámetros de la distribución inicial.

Modelo jerárquico

En este ejemplo se destaca el modelo jerárquico como una estrategia intermedia entre el modelo de unidades iguales y el modelo de unidades independientes. Nos permite reflejar un escenario en donde la información de cada estudio aporta información acerca del parámetro de interés \(\lambda_j\) de los demás estudios sin considerarlos idénticos, de manera que se estima un parámetro \(\lambda_j\) diferente para cada hospital.

Enumeramos algunas de las ventajas potenciales de usar un modelo jerárquico.

  1. Modelo unificado. El problema de elegir entre un modelo de unidades iguales o uno de unidades independientes se resuelve al modelar explícitamente la variabilidad entre las unidades.

  2. Unir fuerzas. Debido al supuesto de intercambiabilidad, al estimar el parámetro de cada unidad se usa información de las demás unidades, conllevando a un encogimiento de la estimación individual hacia la media poblacional, y resulta en una mejor precisión de las estimaciones. La magnitud del encogimiento depende de la varianza entre las unidades, y su efecto resulta particularmente benéfico cuando hay pocas observaciones dentro de una unidad, en cuyo caso hay una gran reducción de la incertidumbre ya que las estimaciones posteriores incorporan la información de otras unidades con menor variabilidad.

  3. Incertidumbre en los parámetros. Asignar una distribución a los hiperparámetros, \(\theta\), nos permite incorporar nuestra incertidumbre sobre la distribución inicial de \(\lambda\).

  4. Cómputo. La estructura jerárquica facilita los cálculos posteriores, debido a que la distribución posterior se factoriza en distribuciones condicionales más sencillas que facilitan la implementación de un muestreador de Gibbs.

Retomando el problema de estimación de tasas de mortalidad por transplante de corazón, modelaremos las \(\lambda_{j}\) con un modelo jerárquico, suponemos intercambiabilidad para reflejar que no contamos con información que nos permita distinguir entre los hospitales.

Primero definimos la distribución de los datos,

\[ \begin{eqnarray} \nonumber y_{j} \sim Poisson(e_{j}\lambda_{j}), \end{eqnarray} \]

donde \(y_{j}\) es el número de muertes observadas, \(e_{j}\) es el número de expuestos y \(\lambda_{j}\) es la tasa de mortalidad para el hospital \(j\), con \(j=1,...,94\).

Después asignamos una distribución al vector de parámetros \(\lambda=(\lambda_{1},...,\lambda_{94})\), para ello suponemos que las tasas de mortalidad \(\{\lambda_{1},...,\lambda_{94}\}\) son una muestra aleatoria de una distribución \(Gamma(\alpha,\alpha/\mu)\) con la forma

\[ \begin{eqnarray} \nonumber g(\lambda_j|\alpha,\mu)=\frac{(\alpha/\mu)^\alpha\lambda_j^{\alpha-1}exp(-\alpha\lambda_j/\mu)}{\Gamma(\alpha)}, \lambda_j>0. \end{eqnarray} \]

La media y varianza iniciales de \(\lambda_{j}\) están dadas por \(\mu\) y \(\mu^2/\alpha\), para toda \(j\). Las llamaremos la media y varianza poblacionales ya que son comunes para todos los hospitales. En la segunda etapa, los hiperparámetros \(\mu\) y \(\alpha\) se suponen independientes y les asignaremos iniciales vagas. Al parámetro de media,

\[ \begin{eqnarray} \nonumber h(\mu)\propto\frac{1}{\mu}, \mu>0. \end{eqnarray} \] Al parámetro de precisión \(\alpha\) le asignamos una densidad inicial propia pero plana, de la forma, \[ \begin{eqnarray} \nonumber h(\alpha)=\frac{z_{0}}{(\alpha+z_0)^2}, \alpha>0. \end{eqnarray} \]

El valor \(z_0\) es la mediana de \(\alpha\), no contamos con información inicial de forma que por ahora usaremos \(z_0=0.5\).

Debido a la estructura de independencia condicional del modelo jerárquico y a que se eligió una inicial conjugada, el análisis posterior es relativamente sencillo. Utilizamos el Teorema de Bayes para calcular la distribución posterior de \(\lambda_{j}\) condicional a los valores de los hiperparámetros \(\mu\) y \(\alpha\),

\[ \begin{eqnarray} \nonumber g(\lambda_{j}|y_{j},\alpha,\mu) &\propto& g(\lambda_{j}|\alpha,\mu)f(y_{j}|\lambda_{j})\\ \nonumber &=& \frac{(\alpha/\mu)^{\alpha} \lambda_{j}^{\alpha-1}exp(-\lambda_{j} \alpha/\mu)}{\Gamma(\alpha)}\frac{exp(-e_{j}\lambda_{j})(e_{j}\lambda_{j})^{y_{j}}}{y_{j}!}\\ \nonumber &\propto& \lambda_{j}^{y_{j}+\alpha-1}exp(-\lambda_{j}(e_{j}+\alpha/\mu)) \end{eqnarray} \]

obtenemos así que las tasas \(\{\lambda_{1},..., \lambda_{94}\}\) tienen distribuciones posteriores independientes \(Gamma(y_{j}+\alpha, e_{j}+\alpha/\mu)\), con media:

\[ \begin{eqnarray} E(\lambda_{j}|y,\alpha,\mu) &=& \frac{y_{j}+\alpha}{e_{j}+\alpha/\mu}\\ \label{eqn:pond} &=& (1-A_{j})\frac{y_{j}}{e_{j}}+A_{j}\mu, \end{eqnarray} \] donde \[ \begin{eqnarray} \label{eqn:factor} A_{j}=\frac{\alpha}{\alpha+e_{j}\mu}, \end{eqnarray} \]

Denominamos al factor \(A_{j}\) como el encogimiento hacia la media poblacional \(\mu\), más adelante derivamos la distribución posterior de \(\mu\), pero por ahora la enunciamos con el propósito de mostrar que su distribución incorpora información de todos los hospitales,

\[ \begin{eqnarray} \nonumber f(\mu|y)=K \int \prod_{i=1}^{94} \left[ \frac {(\alpha/\mu)^\alpha\Gamma(y_{i}+\alpha)} {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}} \right ]\frac{z_{0}}{(\alpha+z_0)^2} \frac{1}{\mu} d\alpha \end{eqnarray} \] donde \(K\) es una constante.

Al escribir la media posterior de \(\lambda_{j}\) como un promedio ponderado, podemos ver el efecto de unir fuerzas mencionado en las observaciones del modelo jerárquico: hay un encogimiento hacia \(\mu\) que depende del número de expuestos, \(e_{j}\), para los hospitales con menor número de expuestos el encogimiento hacia \(\mu\) es mayor, mientras que para aquellos con mayor número de expuestos, es más importante la tasa observada, \(y_{j}/e_{j}\). De esta manera, mayor encogimiento corresponde a las observaciones con mayor incertidumbre.

Notemos también que la factorización de la media posterior es similar a la que obteníamos en el modelo de medias independientes. La diferencia radica en que ahora es un sólo modelo (opuesto a \(94\)), y los parámetros de la distribución inicial de \(\lambda_{j}\) forman parte del modelo de probabilidad pues les asignamos una distribución inicial.

Distribuciones posteriores

Sea \(\lambda=(\lambda_{1},...,\lambda_{94})\) y \(y=(y_{1},...,y_{94})\), calculamos la densidad posterior conjunta de los parámetros,

\[ \begin{align} \nonumber f(\lambda,\alpha,\mu|y) &\propto f(y|\lambda,\alpha,\mu)f(\lambda,\alpha,\mu)\\ \nonumber &= f(y|\lambda) f(\lambda|\alpha,\mu) p(\alpha,\mu)\\ \nonumber &= \prod_{i=1}^{94}f(y_{i}|\lambda_{i}) \prod_{i=1}^{94} f(\lambda_{i}|\alpha,\mu) f(\alpha)f(\mu)\\ \nonumber &\propto \prod_{i=1}^{94} \frac {exp(-e_{i}\lambda_{i})(e_{i}\lambda_{i})^{y_{i}}} {y_{i}!} \prod_{i=1}^{94} \frac {(\alpha/\mu)^{\alpha}\lambda_{i}^{\alpha-1}exp(-\lambda_{i}(\alpha/\mu))} {\Gamma(\alpha)} \frac{z_{0}}{(\alpha+z_0)^2} \frac{1}{\mu}\\ \nonumber &\propto \prod_{i=1}^{94}\frac {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}\lambda_{i}^{y_{i}+\alpha-1}exp(-\lambda_{i}(e_{i}+\alpha/\mu))} {\Gamma(y_{i}+\alpha)} \prod_{i=1}^{94}\frac {(\alpha/\mu)^\alpha\Gamma(y_{i}+\alpha)} {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}}\frac{z_{0}}{(\alpha+z_0)^2} \frac{1}{\mu}, \end{align} \]

de aquí podemos integrar las tasas de mortalidad, \(\lambda_{j}\), para obtener la distribución posterior de los hiperparámetros \(f(\alpha,\mu|y)\),

\[ {\normalsize \begin{align} \nonumber f(\alpha,\mu|y) \propto \int_{\lambda_{1}} ...\lambda_{y_{94}} \prod_{i=1}^{94}(\frac {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}\lambda_{i}^{y_{i}+\alpha-1}exp(-\lambda_{i}(e_{i}+\alpha/\mu))} {\Gamma(y_{i}+\alpha)}) \prod_{i=1}^{94}k_i d\lambda_{1}...d\lambda_{94},\\ \nonumber \end{align} } \]

donde,

\[ \normalsize{ \begin{align} \nonumber k_i = \frac {(\alpha/\mu)^\alpha\Gamma(y_{i}+\alpha)} {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}}\frac{z_{0}}{(\alpha+z_0)^2} \frac{1}{\mu} \end{align}} \]

observemos que las \(\{k_j\}\) no dependen de \(y\) por lo que son constantes en la integral, además para cada \(i\) de la primera multiplicación tenemos una distribución \(Gamma(y_{i}+\alpha, e_{i}+\alpha/\mu)\) por lo que integran \(1\).

Resultando,

\[ {\normalsize \begin{align} \nonumber f(\alpha,\mu|y)=K\prod_{i=1}^{94} \frac {(\alpha/\mu)^\alpha\Gamma(y_{i}+\alpha)} {(e_{i}+\alpha/\mu)^{y_{i}+\alpha}}\frac{z_{0}}{(\alpha+z_0)^2} \frac{1}{\mu} \end{align}} \]

donde \(K\) es la constante de proporcionalidad.

Simulemos de las distribuciones posteriores, para ello procederemos como sigue:

  1. Simulamos \((\mu, \alpha)\) de la distribución marginal posterior.

  2. Simulamos \(\lambda_1,...,\lambda_{94}\) de la distribución posterior condicional a los valores simulados \((\mu, \alpha)\).

Para el primer paso, notamos que ambos parámetros son positivos por lo que es conveniente transformarlos: \(\theta_1 = log(\alpha\), \(\theta_2 = log(\mu)\).

Definimos ahora la distribución posterior en términos de los parámetros transformados

Podemos usar las simulaciones de \(\alpha\) y \(\mu\) para ver el encogimiento de las estimaciones de cada hosiptal hacia la media poblacional. Notamos un mayor encogimiento para los hospitales con menor número de expuestos.

Ahora simualmos observaciones de la distribución posterior de \(\lambda\):

Ya que tenemos las distribuciones posteriores podemos hacer inferencia acerca de la tasa de mortalidad \(\lambda\). A continuación graficamos los intervalos posteriores del \(95\%\) de probabilidad para las tasas \(\lambda_{j}\), el color representa el número de muertes observadas \(y_{j}\), en gris se graficaron las tasas observadas, la gráfica se dividió en \(3\) páneles de acuerdo al número de muertes observadas en los hospitales.

Analizamos la distribución predictiva posterior para la misma muestra de \(10\) hospitales que se utilizó en el modelo de unidades iguales.

Observemos que únicamente en uno de los histogramas el número de muertes observadas se encuentra cerca de la cola de la distribución, lo que indica concordancia de las observaciones con el modelo ajustado.

Finalmente, revisamos la consistencia de los valores observados \(y_{j}\) con la distribución predictiva posterior para todos los hospitales, para ello calculamos la probabilidad de que una observación futura \(y_{j}^*\) sea al menos tan extrema como \(y_{j}\) para todas las observaciones:

\[ \begin{eqnarray} \nonumber P(extremos) = min\{P(y_{j}^*\leq y_{j}),P(y_{j}^*\geq y_{j})\} \end{eqnarray} \]

A continuación graficamos las probabilidades de extremos (calculadas con simulación). Con el modelo jerárquico solamente el \(6\%\) de las probabilidades son menores al \(0.15\), una disminución considerable al \(28\%\) obtenido con el modelo de unidades iguales.

Comparemos ahora las probabilidades de al menos tan extremo usando el modelo jerárquico contra las probabilidades de al menos tan extremo usando el modelo de unidades iguales, los puntos se colorearon de acuerdo al número de expuestos de cada hospital. Observemos que las probabilidades bajo el modelo jerárquico son mayores en la mayoría de los casos.

Stan

Ahora veremos como hacer la estimación usando Stan. Primero hacemos un par de cambios en la definición del modelo, esto es porque la distribución inicial de \(\mu\) no es propia (i.e. no integra uno) y preferimos un ejemplo con iniciales propia, usaremos una \(Gamma\) para \(\alpha\) y \(\mu\) eligiendo parámetros de manera que sean iniciales vagas.

\[\mu, \alpha \sim Gamma(0.01, 0.01).\]

En el modelo definimos una distribución de probabilidad para cada hospital, las tasas de mortalidad se modelan como observaciones de una distribución \(Gamma(\alpha, \alpha/\mu)\), y especificamos la distribución de los hiperparámetros.

Repetimos la gráfica de intervalos posteriores para \(\lambda_j\) (expresada por 1000 expuestos) y tasas observadas.

Referencias

Albert, J. 2009. Bayesian Computation with R. Use R! Springer New York. https://books.google.com.mx/books?id=AALhk\_mt7SYC.