1 Modelo general

Si Su estado de información acerca de las secuencia de variables de conteo \(y_1,\ldots,y_n\) es intercambiable, entonces el modelamiento \(y_1,\ldots,y_n\) admite representación jerárquica de la forma: \[\begin{align} y_i\mid\theta &\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\,,\quad i = 1,\ldots,n \\ \theta &\sim p(\theta) \end{align}\] Este modelo es potencialmente restrictivo por la relación media-varianza: si \(\textsf{E}(y_i\mid\theta) = \textsf{Var}(y_i\mid\theta) = \theta\). Por tal motivo, se recomienda chequear la calidad del modelo en términos de bondad de ajuste (distribución predictiva posterior) y predicción (validación cruzada).

Alternativas:

Así, la distribución muestral (distribución condicional conjunta) de \(\boldsymbol{y} = (y_1,\ldots,y_n)\) dado \(\theta\) está dada por \[ p(\boldsymbol{y}\mid\theta) = \frac{\theta^{s}e^{-n\theta}}{\prod_{i=1}^n y_i!}\,, \] donde \(s = \sum_{i=1}^n y_i\), lo cual sugiere que \(s\) es un estadístico suficiente para \(\theta\).

Por lo tanto la distribución posterior es \[ p(\theta\mid\boldsymbol{y}) \propto \theta^{s}e^{-n\theta}p(\theta)\,. \]

Dado que las \(y_i\)’s son condicionalmente i.i.d. dado \(\theta\) y \(s\) es un estadístico suficiente para \(\theta\), entonces el modelo es equivalente a \[\begin{align*} s\mid\theta &\sim \textsf{Poisson}(n\theta) \\ \theta &\sim p(\theta) \end{align*}\]

2 Modelo Gamma-Poisson

La familia de distribuciones Gamma es conjugada para la distribución muestral Poisson.

La variable aleatoria \(X\) tiene distribución Gamma con parámetros \(\alpha,\beta > 0\), i.e., \(X\mid\alpha,\beta\sim\textsf{Gamma}(\alpha,\beta)\), si su función de densidad de probabilidad es \[ p(x\mid\alpha,\beta) = \frac{\beta^\alpha}{\Gamma(\alpha)}\,x^{\alpha-1}\,e^{-\beta x}\,,\quad x>0\,. \]

Así, el modelo Gamma-Poisson es \[\begin{align*} y_i\mid\theta&\stackrel{\text{iid}}{\sim}\textsf{Poisson}(y_i\mid\theta)\,,\quad i = 1,\ldots,n \\ \theta &\sim \textsf{Gamma}(a,b) \end{align*}\] donde \(a\) y \(b\) con los hiperparámetros del modelo.

Bajo el modelo Gamma-Poison se tiene que la distribución posterior es \[ \theta \mid \boldsymbol{y} \sim \textsf{Gamma}(\theta\mid a + s, b+n)\,, \] donde \(s=\sum_{i=1}^n y_i\), y por lo tanto la media posterior es \[ \textsf{E}(\theta\mid \boldsymbol{y}) = \frac{a+s}{b+n} = \frac{b}{b+n}\cdot \frac{a}{b}+\frac{n}{b+n}\cdot \frac{s}{n}\,, \] la cual es un promedio ponderado del valor esperado previo y la media muestral.

También, la distribución predictiva posterior es Binomial Negativa con parámeros \(a+s\) y \(b+n\), i.e., \(y^*\mid \boldsymbol{y}\sim \textsf{BN}(a+s,b+n)\): \[ p(y^*\mid \boldsymbol{y}) = \frac{\Gamma(y^* +a+s)}{\Gamma(a+s)\Gamma(y^*+1)}\left[\frac{b+n}{b+n+1}\right]^{a+s} \left[\frac{1}{b+n+1}\right]^{y^*}\,,\quad y^*=0,1,\,\ldots. \]

La variable aleatoria \(X\) tiene distribución Binomial Negativa con parámetros \(\alpha,\beta > 0\), i.e., \(X\mid\alpha,\beta\sim\textsf{BN}(\alpha,\beta)\), si su función de masa de probabilidad es \[ p(x\mid\alpha,\beta) = \frac{\Gamma(x+\alpha)}{\Gamma(\alpha)\,\Gamma(x+1)}\,\left[\frac{\beta}{\beta+1}\right]^{\alpha}\,\left[\frac{1}{\beta+1}\right]^x\,,\quad x=0,1,\,\ldots. \]

Por medio de la distribución predictiva posterior se caracterizan diversos aspectos acerca de una observación futura. Por ejemplo, la varianza predictiva \(\textsf{Var}(y^*\mid \boldsymbol{y})\) se puede interpretar como una medida de la incertidumbre posterior acerca de una observación futura \(y^*\). Esto motiva un contraste interesante entre inferencia y predicción: \(\theta\) (el objetivo inferencial) y \(y^*\) (el objetivo predictivo) tienen la misma media posterior, pero la varianza posterior de \(y^*\) es mayor: \[ \textsf{E}(\theta\mid\boldsymbol{y}) = \textsf{E}(y^*\mid\boldsymbol{y}) = \frac{a+s}{b+n}\,, \] mientras que \[ \textsf{Var}(\theta\mid\boldsymbol{y}) = \frac{a+s}{b+n}\left(0 + \frac{1}{b+n}\right) \qquad\text{y}\qquad \textsf{Var}(y^*\mid\boldsymbol{y}) = \frac{a+s}{b+n}\left(1 + \frac{1}{b+n}\right)\,. \]

3 Pruebas de hipótesis

Bajo el paradigma Bayesiano, las hipótesis \(H_0\) y \(H_1\) se consideran como cantidades aleatorias de forma que \(\textsf{Pr}(H_0) + \textsf{Pr}(H_1) = 1\) (típicamente \(\textsf{Pr}(H_0) = \textsf{Pr}(H_1) = 0.5\)). Para probar el sistema simplemente se calculan las probabilidades posteriores de cada una de las hipótesis por medio del teorema de Bayes: \[ \textsf{Pr}(H_k\mid\mathbf{D})=\frac{\textsf{Pr}(\mathbf{D}\mid H_k)\textsf{Pr}(H_k)}{\textsf{Pr}(\mathbf{D}\mid H_0)\textsf{Pr}(H_0) + \textsf{Pr}(\mathbf{D}\mid H_1)\textsf{Pr}(H_1)}\,,\qquad k=0,1\,, \] donde \(\mathbf{D}\) denota los datos disponibles. De esta forma, se tiene que \[ \frac{\textsf{Pr}(H_0\mid\mathbf{D})}{\textsf{Pr}(H_1\mid\mathbf{D})} = \frac{\textsf{Pr}(\mathbf{D}\mid H_0)}{\textsf{Pr}(\mathbf{D}\mid H_1)}\times\frac{\textsf{Pr}(H_0)}{\textsf{Pr}(H_1)}\,, \] es decir, \[ \text{Posibilidades relativas a posteriori} = \text{Factor de Bayes ($B_{01}$)}\times\text{Posibilidades relativas a priori}\,. \] En https://saludpublica.mx/index.php/spm/article/view/5678/6216 se presenta una discusión acerca de la traducción correcta de odds al español.

La cantidad \(\textsf{Pr}(\mathbf{D}\mid H_k)\) se denomina verosimilitud marginal o distribución predictiva previa bajo el modelo especificado por \(H_k\) y se calcula integrando sobre el espacio de parámetros, \[ \textsf{Pr}(\mathbf{D}\mid H_k) = \int p(\mathbf{D}\mid\boldsymbol{\theta}_k, H_k)\,p(\boldsymbol{\theta}_k\mid H_k)\,\text{d}\boldsymbol{\theta}_k\,. \] donde \(\boldsymbol{\theta}_k\) representa los parámetros del modelo bajo \(H_k\).

Se recomienda interpretar el factor de Bayes \(B_{10}\) porque sopesar evidencia en contra de \(H_0\) es más común, pero también se puede hacer la interpretación de términos de evidencia a favor. Siguiendo a Kass (1995), \(B_{10}\) se puede interpretar de la siguiente manera:

\(B_{10}\) Evidencia en contra de \(H_0\)
1 a 3 No vale más que una mención
3 a 20 Moderada
20 a 150 Fuerte
> 150 Decisiva

4 Ejemplo

La Encuesta Social General de 1998 recopiló datos sobre el nivel educativo y el número de hijos de una muestra de mujeres que tenían 40 años de edad. En este ejemplo, se comparan las mujeres con títulos universitarios con aquellas que no los tienen, en función del número de hijos.

# cargar data
Y <- read.csv("GSS1998.txt")

# y1 : numero de hijos, mujeres de 40 años, sin pregrado o menos
# y2 : numero de hijos, mujeres de 40 años, con pregrado o mas
y1 <- Y[(Y$FEMALE == 1) & (Y$YEAR >= 1990) & (Y$AGE == 40) & (Y$DEG <3 ),]$CHILDS
y2 <- Y[(Y$FEMALE == 1) & (Y$YEAR >= 1990) & (Y$AGE == 40) & (Y$DEG >=3),]$CHILDS
y1 <- y1[!is.na(y1)]
y2 <- y2[!is.na(y2)]

# tamaños de muestra
n1 <- length(y1)
print(n1)
## [1] 111
n2 <- length(y2)
print(n2)
## [1] 44
# estadisticos suficientes
s1 <- sum(y1)
print(s1)
## [1] 217
s2 <- sum(y2)
print(s2)
## [1] 66
# relacion media-varianza
r1 <- mean(y1)/var(y1)
print(r1)
## [1] 1.030034
r2 <- mean(y2)/var(y2)
print(r2)
## [1] 1.09322
# distribucion de frecuencias
par(mfrow = c(1,2), mar=c(3,3,1.4,1.4), mgp = c(1.75,0.75,0))
plot(100*table(y1)/n1, type = "h", xlab = "y", ylab = "F. Relativa", col = gray(.5), 
     lwd = 3, ylim = c(0, 50), main = "Menos que pregrado")
plot(100*table(y2)/n2, type = "h", xlab = "y", ylab = "F. Relativa", col = gray(0), 
     lwd = 3, ylim = c(0, 50), main = "Pregrado o más")

4.1 Distribuciones posterior y predictiva posterior

# previa Gamma(2,1)
a <- 2
b <- 1

# parametros de la posterior
ap1 <- a + s1
bp1 <- b + n1
ap2 <- a + s2
bp2 <- b + n2
# grafico
par(mfrow = c(1,2), mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
theta <- seq(0, 5, length = 1000)
plot(NA, NA, xlim = c(0,5), ylim = c(0,3), xlab=expression(theta),
     ylab = expression(paste("p","(",theta," | ",y,")",sep="")),
     main = "Posterior")
lines(theta, dgamma(theta, shape = ap1, rate = bp1), col = 2, lwd = 2)
lines(theta, dgamma(theta, shape = ap2, rate = bp2), col = 4, lwd = 2)
lines(theta, dgamma(theta, shape = a  , rate = b  ), col = 1, lwd = 1)
abline(h = 0, col = 1)
legend("topright", legend = c("Menos que preg.", "Preg. o más", "Previa"), 
       bty = "n", lwd = 2, col = c(2, 4, 1))
y <- 0:12
plot(y - .07, dnbinom(y, size = ap1, mu = ap1/bp1), col = 2, type = "h",
     ylab = "p(y* | y )", xlab = "y*", ylim = c(0, .35), lwd = 3,
     main = "Predictiva posterior")
points(y + .07, dnbinom(y, size = ap2, mu = ap2/bp2), col = 4, type = "h", lwd = 3)

# media posterior e intervalo de credibilidad
tab <- cbind(c(ap1/bp1, qgamma(p = c(.025,.975), shape = ap1, rate = bp1)),
             c(ap2/bp2, qgamma(p = c(.025,.975), shape = ap2, rate = bp2)))
colnames(tab) <- c("Menos que pregrado", "Pregrado o más")
rownames(tab) <- c("Media", "Q2.5%", "Q97.5%")
print(round(t(tab), 3))
##                    Media Q2.5% Q97.5%
## Menos que pregrado 1.955 1.705  2.223
## Pregrado o más     1.511 1.173  1.891
# varianza posterior de theta y de y^*
tab <- cbind(c((a+s1)/(b+n1)*(0+1/(b+1)), ((a+s1)/(b+n1))*(1+1/(b+1))),
             c((a+s2)/(b+n2)*(0+1/(b+1)), ((a+s2)/(b+n2))*(1+1/(b+1))))
colnames(tab) <- c("Menos que pregrado", "Pregrado o más")
rownames(tab) <- c("Var. Parámetro", "Var. Predictiva")
print(round(t(tab), 3))
##                    Var. Parámetro Var. Predictiva
## Menos que pregrado          0.978           2.933
## Pregrado o más              0.756           2.267

4.2 Cálculo de probabilidades

# probabilidad posterior de que theta_j > 2
# ¿como se lleva a cabo el calculo de manera analitica?
set.seed(1234)
th1_mc <- rgamma(n = 10000, shape = ap1, rate = bp1)
th2_mc <- rgamma(n = 10000, shape = ap2, rate = bp2)
print(mean(th1_mc > 2))
## [1] 0.363
print(mean(th2_mc > 2))
## [1] 0.0064
# probabilidad posterior de que y_j^* > 2
set.seed(1234)
y1_mc <- rpois(n = 10000, lambda = th1_mc)
y2_mc <- rpois(n = 10000, lambda = th2_mc)
print(mean(y1_mc > 2))
## [1] 0.3084
print(mean(y2_mc > 2))
## [1] 0.1955

Una evidencia fuerte de una diferencia entre dos poblaciones no implica que la diferencia en términos prácticos también sea grande.

4.3 Prueba de hipótesis

Se quiere probar el sistema de hipótesis \[ H_0: \theta_1 = \theta_2\qquad\text{frente a}\qquad H_1: \theta_1\neq\theta_2\,. \] En este caso se tiene que \[ p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid H_0) = \frac{1}{\prod_{i=1}^{n_1} y_{1,i}!}\,\frac{1}{\prod_{i=1}^{n_2} y_{2,i}!}\,\frac{b^a}{\Gamma(a)}\,\frac{\Gamma(a+s_1+s_2)}{(b+n_1+n_2)^{a+s_1+s_2}} \] mientras que \[ p(\boldsymbol{y}_1,\boldsymbol{y}_2\mid H_1) = \frac{1}{\prod_{i=1}^{n_1} y_{1,i}!}\,\frac{1}{\prod_{i=1}^{n_2} y_{2,i}!}\,\frac{b^a}{\Gamma(a)}\,\frac{b^a}{\Gamma(a)}\,\frac{\Gamma(a+s_1)}{(b+n_1)^{a+s_1}}\,\frac{\Gamma(a+s_2)}{(b+n_2)^{a+s_2}} \] donde \(\boldsymbol{y}_j=(y_{j,1},\ldots,y_{j,n_j})\) para \(j=1,2\), y \(a\) y \(b\) son los hiperparámetros del modelo, y por lo tanto el factor de Bayes correspondiente es \[ B_{01} = \frac{\Gamma(a)}{b^a}\,\frac{\Gamma(a+s_1+s_2)}{\Gamma(a+s_1)\Gamma(a+s_2)}\,\frac{(b+n_1)^{a+s_1}(b+n_2)^{a+s_2}}{(b+n_1+n_2)^{a+s_1+s_2}}\,. \]

# factor de Bayes 10
# calcular en escala log y exponenciar
B01 <- exp(lgamma(a)-a*log(b)+lgamma(a+s1+s2)-lgamma(a+s1)-lgamma(a+s2)+
                   (a+s1)*log(b+n1)+(a+s2)*log(b+n2)-(a+s1+s2)*log(b+n1+n2))
B10 <- 1/B01
print(B10)
## [1] 1.128268
# estimacion puntual de theta_1 - theta_2 
mean(th1_mc - th2_mc)
## [1] 0.4472504
# intervalo de credibilidad para theta_1 - theta_2
quantile(x = th1_mc - th2_mc, probs = c(.025, .975))
##        2.5%       97.5% 
## -0.01438616  0.87496075