3.3 Modelo Binomial negativo
La distribución binomial negativa describe el número de ensayos necesarios para alcanzar un número determinado y fijo de éxitos \(k\) en una secuencia independiente de experimentos tipo Bernoulli. Esta distribución es particularmente útil cuando el parámetro \(\theta\) que se quiere estimar es muy pequeño, como la proporción de una población que padece de alguna enfermedad rara. La razón por la que no se utiliza la distribución binomial es que al fijar el número de ensayos \(n\), con un aprobabilidad \(\theta\) muy pequeña, es muy probable que en en la muestra de tamaño \(n\) no se encuetre ningún paciente con la enfermedad; mientras que al utilizar la distribución binomial negativa, de antemano se garantiza que se obtendrá \(k\) pacientes con la enfermedad en la muestra.
Suponga que \(Y\) es una variable aleatoria cuya distribución es Binomial negativa, y que representa el número de ensayos necesarios \(y\) para alcanzar un número determinado y fijo de éxitos \(k\) en un experimento. La forma funcional de esta distribución es la siguiente \[\begin{equation} p(Y \mid \theta)=\binom{y-1}{k-1}\theta^k(1-\theta)^{y-k}I_{\{k,k+1,\ldots,\}}(y), \end{equation}\]
Así como en la distribución Bernoulli y Binomial, el parámetro \(\theta\) está restringido al espacio \(\Theta=[0,1]\). Luego, es admisible proponer que \(\theta\) siga una distribución Beta. Por tanto, la distribución previa del parámetro \(\theta\) está dada por la expresión (3.1). Bajo este marco de referencia se tienen los siguientes resultados
En algunas situaciones se puede encontrar una muestra de variables con distribución binomial negativa. Por ejemplo, la entrevista de pacientes para encontrar cierta enfermedad puede llevarse a cabo en diferentes puntos de atención médica o en diferentes ciudades del país. Así en cada punto de atención, se tendrá el dato correspondiente a una variable con distribución binomial negativa. El procedimiento inferencial bayesiano para estas situaciones se describe a continuación:
Ejemplo 3.4 Una franquicia de investigación farmacéutica ha desarrollado un nuevo tratamiento farmacológico sobre pacientes diabéticos que padezcan, a su vez, de enfermedades cardíacas o cardiopatías (angina de pecho, infarto de miocardio, insuficiencia mitral, estenosis mitral, entre otras). Para evaluar el nuevo tratamiento, es necesario seleccionar una muestra, mediante el diseño de un experimento clínico, de pacientes que tienen estas características.
Por otro lado, se sabe que la proporción de personas que padecen de diabetes y que además tienen algún tipo de condición cardiaca es muy baja y es necesario obtener una estimación precisa de la proporción de personas con estas condiciones. Con base en lo anteriormente expuesto, se puede pensar en seleccionar una muestra grande de personas y utilizar un acercamiento binomial para estimar esta proporción. Sin embargo, dado que la prevalencia de esta condición es bastante baja, es posible que el número de personas en la muestra que presenten estas enfermedades sea nulo; por consiguiente, la estimación binomial no será, de ninguna forma, precisa.
Por lo tanto, el diseño clínico está supeditado al uso de la distribución Binomial Negativa, en donde se entrevistarán pacientes, de una base de datos de un hospital de la ciudad asociado con la franquicia, hasta conseguir una muestra de cinco pacientes que padezcan de estas condiciones. Después de varios meses de entrevistas, se encontró el quinto paciente en la entrevista número 1106.
Mediante el análisis bayesiano, suponiendo una distribución previa \(Beta(0.5, 0.5)\), se llega a que la distribución posterior del parámetros \(\theta\) es \(Beta(0.5+5, 0.5+1106-5)=Beta(5.5, 1101.5)\). Por lo tanto, la estimación puntual del parámetro de interés, que corresponde a la media de la distribución posterior, es 0.0049, que equivale una proporción de 0.49% de personas con estas enfermedades.
El siguiente código computacional muestra cómo se puede llegar a las mismas conclusiones conSTAN, haciendo la salvedad de que STAN define esta distribución en términos del número de fracasos \(m = y - k\) necesarios para obtener \(k\) éxitos.
BinNegativa <- 'data {
int<lower=0> k;
int<lower=0> y;
}
transformed data {
int<lower=0> m;
m = y - k;
}
parameters {
real<lower=0> beta;
}
transformed parameters {
real<lower=0> theta;
theta = beta/(beta + 1);
}
model {
m ~ neg_binomial(k, beta);
theta ~ beta(0.5, 0.5);
}
'
y <- 1106
k <- 5
sample_data <- list(k = k, y = y)
BNfit <- stan(model_code = BinNegativa,
data = sample_data, verbose = FALSE)La siguiente salida de STAN permite conocer la estimación bayesiana posterior y los límites del intervalo de credibilidad al 95%.
print(BNfit, pars = "theta",
digits = 4, probs = c(0.025, 0.975))## Inference for Stan model: ca45da8b0e94b143378403f155105e09.
## 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.0049 1e-04 0.0021 0.0016 0.0099 1583 1.0009
##
## Samples were drawn using NUTS(diag_e) at Sun Jun 6 23:43:20 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).Después de las iteraciones necesarias, la salida del anterior código muestra la estimación puntual dada por 0.00498 y un intervalo de credibilidad al 95%, dado por \((0.00174, 0.01013)\).
Ejemplo 3.5 Continuando con la temática del ejemplo anterior, suponga que la franquicia llevó a cabo la misma investigación en las 31 ciudades con mayor densidad poblacional de país. En total, se tuvieron 29620 entrevistas para un total de éxitos de 152, tal como se muestra a continuación.
| Ciudad | y | k |
|---|---|---|
| BOGOTA | 1001 | 4 |
| MEDELLIN | 978 | 6 |
| CALI | 999 | 5 |
| BARRANQUILLA | 860 | 4 |
| CARTAGENA | 1155 | 4 |
| CUCUTA | 585 | 6 |
| BUCARAMANGA | 1030 | 3 |
| IBAGUE | 960 | 5 |
| SOLEDAD | 1002 | 6 |
| SANTA MARTA | 763 | 7 |
| SOACHA | 1036 | 5 |
| PASTO | 779 | 5 |
| MONTERIA | 1158 | 4 |
| VILLAVICENCIO | 1017 | 5 |
| BELLO | 888 | 6 |
| MANIZALES | 977 | 4 |
| VALLEDUPAR | 1256 | 6 |
| BUENAVENTURA | 1349 | 6 |
| NEIVA | 1047 | 5 |
| PALMIRA | 1088 | 5 |
| ARMENIA | 649 | 3 |
| POPAYAN | 765 | 4 |
| FLORIDABLANCA | 699 | 5 |
| SINCELEJO | 1042 | 4 |
| ITAGUI | 1212 | 5 |
| BARRANCABERMEJA | 660 | 5 |
| TULUA | 671 | 5 |
| ENVIGADO | 835 | 6 |
| DOSQUEBRADAS | 997 | 5 |
| RIOHACHA | 1146 | 4 |
| SINCELEJO | 1016 | 5 |
STAN
BinNegativa2 <- 'data {
int<lower=0> n;
int<lower=0> k[n];
int<lower=0> y[n];
}
transformed data {
int<lower=0> m[n];
for(i in 1:n){
m[i] = y[i] - k[i];
}
}
parameters {
real<lower=0> b;
}
transformed parameters {
real<lower=0> theta;
theta = b/(b + 1);
}
model {
for(i in 1:n){
m[i] ~ neg_binomial(k[i], b);
}
theta ~ beta(0.5, 0.5);
}
'
y <- c(1001, 978, 999, 860, 1155, 585, 1030,
960, 1002, 763, 1036, 779, 1158, 1017,
888, 977, 1256, 1349, 1047, 1088, 649,
765, 699, 1042, 1212, 660, 671, 835,
997, 1146, 1016)
k <- c(4, 6, 5, 4, 4, 6, 3, 5, 6, 7, 5, 5, 4,
5, 6, 4, 6, 6, 5, 5, 3, 4, 5, 4, 5, 5,
5, 6, 5, 4, 5)
sample_data <- list(k = k, y = y, n = length(y))
BNfit2 <- stan(model_code = BinNegativa2,
data = sample_data, verbose = FALSE)Después de cinco mil iteraciones, la salida del anterior código muestra la estimación puntual dada por 0.00515 y un intervalo de credibilidad al 95%, dado por \((0.00439, 0.00603)\), mucho más estrecho que el intervalo de credibilidad del anterior ejemplo
print(BNfit2, pars = "theta",
digits = 4, probs = c(0.025, 0.975))## Inference for Stan model: 756b8df22716f40c50aa8045790299a4.
## 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.0052 0 4e-04 0.0044 0.006 1302 1.0032
##
## Samples were drawn using NUTS(diag_e) at Sun Jun 6 23:43:53 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).La figura 3.9 muestra la distribución posterior para este ejemplo, junto con la estimación puntual, correspondiente a la media.
bayesplot::mcmc_areas(BNfit2, pars = "theta",
prob = 0.95)
Figura 3.9: Distribución posterior.
Una vez observados los datos actuales y encontrada la distribución posterior, se puede encontrar la distribución predictiva posterior de una nueva variable con distribución binomial negativa. Es decir, se puede definir el mecanismo probabilístico para el número de ensayos necesarios para encontrar \(\tilde{k}\) éxitos.
Resultado 3.6 Después de la recolección de datos, la distribución predictiva posterior para una nueva variable \(\tilde{Y}\) está dada por
\[\begin{equation*} p(\tilde{Y}|Y_1,\cdots,Y_n)=\binom{\tilde{y}-1}{\tilde{k}-1}\frac{Beta(\alpha+\tilde{k}+\sum k_i,\beta+\tilde{y}-\tilde{k}+\sum y_i-\sum k_i)}{Beta(\alpha+\sum k_i,\beta+\sum y_i-\sum k_i)}I_{\{\tilde{k},\tilde{k}+1,\cdots\}}(\tilde{y}) \end{equation*}\]Ejemplo 3.6 Siguiendo con los datos del ejemplo 3.5, suponga que se quiere recolectar información de tres pacientes con cardiopatía en cierta ciudad, y se quiere conocer acerca del número de entrevistas necesarias para . Utilizando la distribución previa \(Beta(0.5,0.5)\) y los datos de las 31 ciudades del ejemplo, se tiene que la distribución predictiva para el número de entrevistas necesarias para encontrar 3 pacientes está dada por
\[\begin{align*} &\ \ \ \ p(\tilde{Y}|Y_1,\cdots,Y_n)\\ &=\binom{\tilde{y}-1}{4}\frac{Beta(0.5+5+152,0.5+\tilde{y}-5+29620-152)}{Beta(0.5+152,0.5+29620-152)}I_{\{5,6,\cdots\}}(\tilde{y})\\ &=\binom{\tilde{y}-1}{4}\frac{Beta(157.5,\tilde{y}+29463.5)}{Beta(152.5,29468.5)}I_{\{5,6,\cdots\}}(\tilde{y}) \end{align*}\]
Con los siguientes códigos se puede calcular la anterior función predictiva.BNpred <- function(y, alfa, beta, s, n, k){
choose(y - 1, k - 1) *
exp(lbeta(alfa + k + s, beta + y - k + n - s) -
lbeta(alfa+s,beta+n-s))
}
alfa <- beta <- 0.5
s <- sum(k)
n <- sum(y)
k <- 5
fun <- rep(NA)
for(y in 5:5000){
fun[y - 4] <- BNpred(y, alfa, beta, s, n, k)
}
sum(fun)## [1] 0.9999994
Figura 3.10: Distribución predictiva posterior para el número de entrevistas necesarias para encontrar 5 pacientes.
Se puede ver que el número de entrevistas que tiene mayor probabilidad asociadas es el valor 768, usando el comando which(fun == max(fun)). También, se puede observar que la probabilidad de que en menos de 500 entrevistas se encuentren los 5 pacientes es de solo el 0.1200985 usando el comando sum(fun[1:(500 - 4)])
Nótese que es posible también asignar una previa informativa \(Beta(5.5, 1101.5)\), que da cuenta de la información del estudio del ejemplo anterior.↩︎