3.5 Modelo Exponencial

Suponga que \(\mathbf{Y}=\{Y_1,\ldots,Y_n\}\) corresponde a una muestra de variables aleatorias con distribución Exponencial. Luego, la función de distribución conjunta o verosimilitud está dada por

\[\begin{align} p(\mathbf{Y} \mid \theta)&=\prod_{i=1}^n\theta e^{(-\theta y)}I_{(0,\infty)}(y_i) \notag \\ &=\theta^n e^{(-\theta \sum_{i=1}^ny_i)}I_{(0,\infty)^n}(y_1,\ldots,y_n) \end{align}\]

Donde \(\{0,1\ldots\}^n\) denota el producto cartesiano \(n\) veces sobre el intervalo \((0,\infty)\). Por otro lado, como el parámetro \(\theta\) está restringido al espacio \(\Theta=(0,\infty)\), entonces es posible formular varias opciones para la distribución previa del parámetro, al igual que en la distribución Poisson. Así mismo, suponga que la distribución previa para el parámetro de interés es la distribución Gamma tal como aparece en la expresión (3.7). Bajo este marco de referencia se tienen los siguientes resultados

Resultado 3.11 La distribución posterior del parámetro \(\theta\) sigue una distribución \[\begin{equation*} \theta \mid \mathbf{Y} \sim Gamma\left(\alpha+n,\beta+\sum_{i=1}^ny_i\right) \end{equation*}\]

Prueba. \[\begin{align*} p(\theta \mid \mathbf{Y})&\propto p(\mathbf{Y} \mid \theta)p(\theta \mid \alpha,\beta)\\ &=\theta^n e^{(-\theta \sum_{i=1}^ny_i)}I_{(0,\infty)^n}(y_1,\ldots,y_n)\frac{\beta^\alpha \theta^{\alpha-1} e^{-\beta\theta}}{\Gamma(\alpha)}I_{(0,\infty)}(\theta)\\ &\propto \theta^{\alpha+n-1}e^{-(\beta+\sum_{i=1}^ny_i)}I_{(0,\infty)}(\theta) \end{align*}\]

Por lo tanto, factorizando convenientemente, se encuentra una expresión idéntica a la función de distribución de una variable aleatoria con distribución \(Gamma(\alpha+n,\beta+\sum_{i=1}^ny_i)\).

Resultado 3.12 La distribución predictiva previa para una observación \(\mathbf{y}=\{y_1,\ldots,y_n\}\) de la muestra aleatoria está dada por

\[\begin{equation} p(\mathbf{Y})=\frac{\Gamma(\alpha+n)}{\Gamma(\alpha)}\frac{\beta^\alpha}{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}} I_{(0,\infty)^n}(y_1,\ldots,y_n) \end{equation}\]

y define una auténtica función de densidad de probabilidad continua.

Prueba. De la definición de función de distribución predictiva se tiene que

\[\begin{align*} p(\mathbf{Y})&=\int p(\mathbf{Y} \mid \theta)p(\theta \mid \alpha,\beta)\ d\theta\\ &=\int_0^{\infty}\theta^n e^{(-\theta \sum_{i=1}^ny_i)}I_{(0,\infty)^n}(y_1,\ldots,y_n)\frac{\beta^\alpha \theta^{\alpha-1} e^{-\beta\theta}}{\Gamma(\alpha)} \ d\theta\\ &=\frac{\Gamma(n+\alpha)}{\Gamma(\alpha)}\frac{\beta^\alpha}{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}I_{(0,\infty)^n}(y_1,\ldots,y_n)\\ &\hspace{2cm}\times \int_0^{\infty} \frac{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{\Gamma(n+\alpha)} \theta^{\alpha+n-1}e^{-(\beta+\sum_{i=1}^ny_i)\theta} \ d\theta\\ &=\frac{\Gamma(\alpha+n)}{\Gamma(\alpha)}\frac{\beta^\alpha}{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}I_{(0,\infty)^n}(y_1,\ldots,y_n) \end{align*}\]

Por ejemplo, en el caso en que la muestra aleatoria estuviera constituida por una sola variable aleatoria, entonces no es difícil ver, utilizando la definición de la función matemática Gamma, que la función de distribución predictiva estaría dada por

\[\begin{align*} p(Y)&=\frac{\Gamma(\alpha+1)}{\Gamma(\alpha)}\frac{\beta^\alpha}{(\beta+y)^{\alpha+1}} I_{(0,\infty)}(y)\notag \\ &=\frac{\alpha \beta^\alpha}{(\beta+y)^{\alpha+1}} I_{(0,\infty)}(y) \end{align*}\]

Para chequear la convergencia de la anterior distribución es necesario recurrir a los resultados del cálculo integral. Dado que el espacio de muestreo de la variable aleatoria \(Y\) es el intervalo \((0,\infty)\), entonces la integral definida es igual a uno, lo que conlleva a que en este caso particular \(P(Y)\) sea una auténtica función de densidad de probabilidad.

\[\begin{align*} \int_0^{\infty}p(Y)\ dy=\int_0^{\infty}\frac{\alpha \beta^\alpha}{(\beta+y)^{\alpha+1}} \ dy = \beta^\alpha \left[\frac{(\beta+y)^{-\alpha}}{-\alpha} \right]_0^{\infty} =1 \end{align*}\]

Volviendo al caso general en donde se tiene una muestra aleatoria, se tiene el siguiente resultado.

Resultado 3.13 Después de la recolección de los datos, la distribución predictiva posterior para una conjunto de nuevas variables aleatorias \(\tilde{\mathbf{y}}=\{\tilde{y}_1,\ldots,\tilde{y}_{n^*}\}\), de tamaño \(n^*\), está dada por

\[\begin{align} p(\tilde{\mathbf{y}} \mid \mathbf{Y})&= \frac{\Gamma(n+\alpha+n^{*})}{\Gamma(n+\alpha)} \frac{(\beta+\sum_{i=1}^ny_i)^{n+\alpha}}{(\sum_{i=1}^{n^*}\tilde{y}_i+\beta+\sum_{i=1}^ny_i)^{n^*+\alpha+n}}\notag\\ &\hspace{4cm}\times I_{(0,\infty)^{n^*}}(\tilde{y}_1,\ldots,\tilde{y}_n) \end{align}\]

El anterior resultado permite calcular la distribución predictiva conjunta de variables aleatorias por observar. En algunas situaciones lo que se quiere pronosticar es el comportamiento probabilístico de promedio muestral de este conjunto de variables aleatorias; es decir, \(\bar{Y}^*=\sum_{i=1}^{n^*}\tilde{Y}_i\). En el siguiente resultado se presenta la distribución predictiva de esta variable aleatoria.

Resultado 3.14 Después de la recolección de los datos, la distribución predictiva posterior para el promedio muestral de un nuevo conjunto de variables aleatorias \(\bar{Y}^*=\sum_{i=1}^{n^*}\tilde{Y}_i\) está dada por

\[\begin{equation*} p(\bar{Y}^*)=\frac{n^*\Gamma(n^*+\alpha+n)}{\Gamma(n^*)\Gamma(\alpha+n)}\frac{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{(n^*\bar{Y}^*+\beta+\sum y_i)^{n^*+\alpha+n}}(n^*\bar{Y}^*)^{n^*-1}I_{(0,\infty)}(\bar{Y}^*) \end{equation*}\]

Prueba. En primer lugar se halla la distribución predictiva posterior de la variable \(\tilde{S}=\sum_{i=1}^{n^*}\tilde{Y}_i\), teniendo en cuenta que \(\tilde{S}|\theta\sim Gamma(n^*,\theta)\), de esta forma

\[\begin{align*} p(\tilde{S}|\mathbf{Y})&=\int p(\tilde{S}|\theta)p(\theta|\mathbf{Y})\ d\theta\\ &=\int_0^{\infty} \frac{\theta^{n^*}}{\Gamma(n^*)}\tilde{S}^{n^*-1}e^{-\theta\tilde{S}}I_{(0,\infty)}(\tilde{S})\frac{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{\Gamma(\alpha+n)}\theta^{{\alpha+n-1}}e^{-(\beta+\sum y_i)\theta}d\theta\\ &=\frac{\tilde{S}^{n^*-1}(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{\Gamma(n^*)\Gamma(\alpha+n)}I_{(0,\infty)}(\tilde{S})\int_0^{\infty} \theta^{n^*+\alpha+n-1}e^{-(\tilde{S}+\beta+\sum y_i)\theta}\ d\theta\\ &=\frac{\tilde{S}^{n^*-1}(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{\Gamma(n^*)\Gamma(\alpha+n)}\frac{\Gamma(n^*+\alpha+n)}{(\tilde{S}+\beta+\sum y_i)^{n^*+\alpha+n}}I_{(0,\infty)}(\tilde{S}) \end{align*}\]

Al aplicar el teorema de transformación a la distribución predictiva, se puede hallar la distribución de \(\bar{Y}^*\), dada por

\[\begin{align*} p(\bar{Y}^*|\mathbf{Y})=\frac{n^*\Gamma(n^*+\alpha+n)}{\Gamma(n^*)\Gamma(\alpha+n)}\frac{(\beta+\sum_{i=1}^ny_i)^{\alpha+n}}{(n^*\bar{Y}^*+\beta+\sum y_i)^{n^*+\alpha+n}}(n^*\bar{Y}^*)^{n^*-1}I_{(0,\infty)}(\bar{Y}^*) \end{align*}\]

En la práctica puede ocurrir que algunos de los valores de \(n\), \(n^*\), \(\sum_{i=1}^ny_i\) y \(n^*\bar{Y}^*\) sean muy grandes; por consiguiente, evaluar directamente la expresión anterior puede ocasionar problemas numéricos. Realizando algunas operaciones algebraicas, se encuentra la siguiente expresión equivalente para la distribución predictiva posterior de \(\bar{Y}^*\) que evita problemas numéricas:

\[\begin{equation} \tag{3.13} p(\bar{Y}^*|\mathbf{Y})=\frac{1}{\bar{Y}^*Beta(n,n^*)}\left(\frac{\beta+\sum_{i=1}^ny_i}{\beta+\sum_{i=1}^ny_i+n^*\bar{Y}^*}\right)^{\alpha+n}\left(\frac{n^*\bar{Y}^*}{\beta+\sum_{i=1}^ny_i+n^*\bar{Y}^*}\right)^{n^*}I_{(0,\infty)}(\bar{Y}^*) \end{equation}\]

Por otro lado, se puede ver que al utilizar la distribución previa no informativa de Jeffreys, la distribución predictiva posterior de \(\bar{Y}^*\) está dada por

\[\begin{equation} \tag{3.14} p(\bar{Y}^*|\mathbf{Y})=\frac{n^*\Gamma(n^*+n)}{\Gamma(n^*)\Gamma(n)}\frac{(\sum_{i=1}^ny_i)^n}{(n^*\bar{Y}^*+\sum y_i)^{n^*+n}}(n^*\bar{Y}^*)^{n^*-1}I_{(0,\infty)}(\bar{Y}^*) \end{equation}\]

La cual es equivalente a la siguiente expresión que en ocasiones puede ser útil para evitar problemas numéricos

\[\begin{equation} \tag{3.15} p(\bar{Y}^*|\mathbf{Y})=\frac{1}{\bar{Y}^*Beta(n,n^*)}\left(\frac{\sum_{i=1}^ny_i}{\sum_{i=1}^ny_i+n^*\bar{Y}^*}\right)^n\left(\frac{n^*\bar{Y}^*}{\sum_{i=1}^ny_i+n^*\bar{Y}^*}\right)^{n^*}I_{(0,\infty)}(\bar{Y}^*) \end{equation}\]

Ejemplo 3.8 Crowley y Hu (1977) reportan un conjunto de datos que da cuenta de los tiempos de sobrevivencia de \(n=69\) miembros del programa de transplante de corazón de Stanford (los tiempos se reportan en días después del transplante). Los datos pueden ser encontrados en el paquete survival Therneau y Lumley (2011) de R, mediante la implementación del siguiente código computacional.

library(survival)
library(dplyr)
data(heart)

sobrevida <- heart %>%
  filter(transplant == 1) %>%
  mutate(tiempo = stop - start)

A continuación, se muestran los primeros y últimos datos de este estudio. Se recuerda que el total de pacientes atendidos en este estudio fue de \(n=69\) y la suma de los tiempos de sobrevida es de \(\sum_{i=1}^ny_i=25998.5\).

sobrevida %>% {
  rbind(head(., 5), tail(., 5))
}

	start	stop	event	age	year	surgery	transplant	id	tiempo
1	1	16	1	6.297057	0.2655715	0	1	3	15
2	36	39	1	-7.737166	0.4900753	0	1	4	3
3	51	675	1	2.869268	0.7802875	0	1	7	624
4	12	58	1	-5.497604	0.8624230	0	1	10	46
5	26	153	1	-0.019165	0.8733744	0	1	11	127
65	2	16	1	-7.718001	5.9767283	0	1	95	14
66	13	180	0	-21.349760	6.0095825	0	1	96	167
67	21	131	0	-24.383299	6.1437372	0	1	97	110
68	96	109	0	-19.370294	6.2039699	0	1	98	13
69	38	39	0	-12.939083	6.3956194	1	1	100	1

Estos tiempos de sobrevivencia pueden ser modelados mediante una distribución exponencial. Además de inferir acerca del parámetro de esta distribución, también es posible inferir acerca del tiempo promedio de sobrevivencia de un individuo sometido a este tipo de transplantes. Luego, dadas las implicaciones del estudio, se debe ser muy cuidadoso en la asignación de los parámetros de la distribución previa. Una forma de hacerlo es asignar valores muy pequeños a estos parámetros. Otra forma de hacerlo es utilizando la distribución previa de Jeffreys, que corresponde a una distribución impropia y que conduce a resultados muy cercanos a los del enfoque anterior.

Utilizando parámetros previos muy cercanos a cero, la distribución posterior del parámetro de interés es \(Gamma(69, 25998.5)\). Como es bien sabido, una estimación bayesiana para el parámetro \(\theta\) está dada por la media de esta distribución posterior, la cual equivale a \(69/25998.5=0.0026\). Ahora, como la esperanza de la distribución exponencial es \(1/\theta\), entonces el tiempo promedio de sobrevivencia es de \(1/0.0026=376.78\) días. Sin embargo, en este tipo de estudios es común que se presenten muchos datos atípicos; por ende el promedio no es una medida de escala válida en este tipo de análisis, puesto que no es una medida robusta y se prefiere la utilización de la mediana. El siguiente código computacional en STAN puede ser usado para realizar inferencias sobre el parámetro \(\theta\), sobre el tiempo promedio y el tiempo mediano. De la misma forma, es posible obtener intervalos de credibilidad para estos parámetros.

Exponencial <- '
data {
  int<lower=0> n;
  vector[n] y;
}
parameters {
  real<lower=0> theta;
}
transformed parameters {
  real<lower=0> invtheta = 1/theta;
}
model {
  y ~ exponential(theta);
  theta ~ gamma(0.1, 0.1);
}
'
sample_data <- list(y = sobrevida$tiempo, 
                    n = nrow(sobrevida))
Expofit <- stan(model_code = Exponencial,
                data = sample_data, verbose = FALSE)

Después de las iteraciones, los resultados de este código muestran una estimación para \(\theta\) de 0.0026 con un intervalo de credibilidad de (0.00208, 0.00332). Para la media \(1/\theta\), se tiene una estimación puntual de 381.4 con un intervalo de credibilidad de (302.2, 485.8). La mediana se estimó en 376.8 días de sobrevivencia.

print(Expofit, digits = 4, 
      pars = c("theta", "invtheta"),
      probs = c(0.025, 0.975))

## Inference for Stan model: 1e6dccb88d5f711b516d61c512b5ced4.
## 4 chains, each with iter=2000; warmup=1000; thin=1; 
## post-warmup draws per chain=1000, total post-warmup draws=4000.
## 
##              mean se_mean      sd     2.5%    97.5% n_eff   Rhat
## theta      0.0027  0.0000  0.0003   0.0021   0.0033  1634 0.9999
## invtheta 381.9284  1.1482 46.7530 301.9569 483.5668  1658 0.9999
## 
## Samples were drawn using NUTS(diag_e) at Sun Jun  6 23:45:23 2021.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Las figuras 3.15 y 3.16 muestra la distribución posterior para este ejemplo, junto con la estimación puntual, correspondiente a la media.

bayesplot::mcmc_areas(Expofit, pars = "theta", 
                      prob = 0.95)

Figura 3.15: Distribución posterior.

bayesplot::mcmc_areas(Expofit, pars = "invtheta", 
                      prob = 0.95)

Figura 3.16: Distribución posterior.

Ejemplo 3.9 Suponga ahora que se va a realizar el trasplante de corazón a 5 pacientes, y se quiere conocer el comportamiento probabilístico del tiempo promedio de sobrevida en estos 5 pacientes. Aplicando la distribución predictiva y definiendo una distribución previa no informativa de Jeffreys, se tiene que

\[\begin{align*} p(\bar{Y}^*|\mathbf{Y})&=\frac{5\Gamma(5+69)}{\Gamma(5)\Gamma(69)}\frac{25998.5^{69}}{(5\bar{Y}^*+25998.5)^{5+69}}(5\bar{Y}^*)^4\\ &=\frac{1}{\bar{Y}^*Beta(5,69)}\left(\frac{25998.5}{5\bar{Y}^*+25998.5}\right)^{69}\left(\frac{5\bar{Y}^*}{5\bar{Y}^*+25998.5}\right)^5 \end{align*}\]

El cálculo de esta función predictiva se puede llevar a cabo con el siguiente código en R, además de comprobar que la integral de la función es 1.

pred_exp <- function(x){
  ((s/(s+x*n.mono))^n) * 
    ((x*n.mono/(s+x*n.mono))^n.mono) / 
    (x*beta(n,n.mono))
}


alfa <- beta <- 0
s <- 25998.5
n <- 69
n.mono <- 5
integrate(pred_exp, 0.0001, 10000)

## 1 with absolute error < 3.2e-10

La distribución predictiva de esta función se puede visualizar en la Figura 3.17, donde se puede ver que la mayor masa de la función se acumula alrededor de los 300 días. Usando el comando integrate(pred_exp, 800, 10000), también se puede observar que la probabilidad de que en promedio los cinco pacientes sobrevivan más de 800 días es de 0.0264485.

Figura 3.17: Distribución predictiva posterior para el tiempo promedio de sobrevivencia de transplante de corazón.

Referencias

Crowley, J., y M. Hu. 1977. «Covariance analysis of heart transplant survival data». Journal of the American Statistical Association 72: 27-36.

Therneau, T., y T. Lumley. 2011. survival: Survival analysis, including penalised likelihood.