Uno de los procedimientos más estándar para estimar cantidades desconocidas es el método de máxima verosimilitud. Los estimadores de máxima verosimilitud tienen propiedades convenientes, y dan resultados buenos si los modelos probabilísticos
Máxima verosimilitud es un proceso intuitivo, y consiste en aprender o estimar valores de parámetros desconocidos suponiendo para los datos su explicación más probable. Para esto, usando supuestos y modelos requeriremos calcular la probabilidad de un conjunto de observaciones.
Supongamos que una máquina produce dos tipos de bolsas de 25 galletas: la mitad de las veces produce una bolsa con 5 galletas de avena y 20 de chispas de chocolate, y la otra mitad produce galletas con 23 galletas de avena y 2 de chispas de chocolate.
Tomamos una bolsa, y no sabemos qué tipo de bolsa es (parámetro desconocido). Extraemos al azar una de las galletas, y es de chispas de chocolate (observación).
Por máxima verosimilitud, inferimos que la bolsa que estamos considerando tiene 5 galletas de avena. Esto es porque es más probable observar una galleta de chispas en las bolsas que contienen 5 galletas de avena que en las bolsas que contienen 23 galletas de avena.
Cómo se aprecia en el ejemplo anterior, el esquema general es:
\[x_1, x_2, \ldots, x_n\]
\[P(x_1,x_2,\ldots, x_n|\theta)\]
En principio, los modelos que consideramos pueden ser complicados y tener varias partes o parámetros. Veamos primero un ejemplo clásico con un solo parámetro, y cómo lo resolveríamos usando máxima verosimilitud.
Nota: Cuando decimos muestra en general nos referimos a observaciones independientes obtenidas del mismo proceso (ver la sección anterior para ver qué significa que sea independientes). Este esquema es un supuesto que simplifica mucho los cálculos, como discutimos antes. Muchas veces este supuesto sale del diseño de la muestra o del estudio, pero en todo caso es importante considerar si es razonable o no para nuestro problema particular.
Supongamos que queremos saber qué proporción de registros de una base de datos tiene algún error menor de captura. No podemos revisar todos los registros, así que tomamos una muestra de 8 registros, escogiendo uno por uno al azar de manera independiente. Revisamos los 8 registros, y obtenemos los siguientes datos:
\[x_1 = 0, x_2 = 1, x_3 = 0, x_4 = 0, x_5 =1, x_6 =0, x_7 =0, x_8 =0\]
donde 1 indica un error menor. Encontramos dos errores menores. ¿Cómo estimamos el número de registros con errores leves en la base de datos?
Ya sabemos una respuesta razonable para nuestro estimador puntual, que sería \(\hat{p}=2/8=0.25\). Veamos cómo se obtendría por máxima verosimilitud.
Según el proceso con el que se construyó la muestra, debemos dar una probabilidad de observar los 2 errores en 8 registros. Supongamos que en realidad existe una proporción \(p\) de que un registro tenga un error. Entonces calculamos la probabilidad de observar la muestra dado p:
\[P(X_1 = 0, X_2 = 1, X_3 = 0, X_4 = 0, X_5 =1, X_6 =0, X_7 =0, X_8 =0)\]
que es igual a
\[ \begin{aligned} &P(X_1 = 0)P(X_2 = 1)P(X_3 = 0)P( X_4 = 0) \times \\ &P(X_5 =1)P(X_6 =0)P(X_7 =0)P(X_8 =0) \end{aligned}\]
pues la probabilidad de que cada observación sea 0 o 1 no depende de las observaciones restantes (la muestra se extrajo de manera independiente).
Nótese entonces que para cualquier \(X_i\),
\[P(X_i=1) = p, P(X_i=0) = 1-p,\]
así que la cantidad de arriba es igual a
\[(1-p)p(1-p)(1-p)p(1-p)(1-p)(1-p) \]
que se simplifica a
\[ L(p) = p^2(1-p)^6\]
Ahora la idea es encontrar la p que maximiza la probabilidad de lo que observamos. En este caso se puede hacer con cálculo, pero vamos a ver una gráfica de esta función y cómo resolverla de manera numérica.
verosimilitud <- function(p){
p^2 * (1-p)^6
}
dat_verosim <- tibble(x = seq(0,1, 0.001)) |> mutate(prob = map_dbl(x, verosimilitud))
ggplot(dat_verosim, aes(x = x, y = prob)) + geom_line() +
geom_vline(xintercept = 0.25, color = "red") +
xlab("p")
Nótese que esta gráfica:
Obsérvese que para hacer esto usamos dos partes del proceso generador de datos:
y resolvimos el problema de estimación convirtiéndolo en uno de optimización.
Para confirmar nuestro resultado, escribimos la función de verosimilitud
crear_verosim <- function(datos_error){
verosim <- function(p){
n <- length(datos_error)
err <- sum(datos_error)
# nota: es mejor usar el logaritmo de la verosimilitud
res <- (p^err) * ((1 - p)^(n - err))
res
}
verosim
}
datos_error <- c(0, 1, 0, 0, 1, 0, 0 ,0)
verosim <- crear_verosim(datos_error)
verosim(0.50)[1] 0.00390625verosim(0.25)[1] 0.01112366Y ahora maximizamos:
estimador_mv <- function(datos){
verosim <- crear_verosim(datos)
res <- optimize(verosim, interval = c(0.0001, 0.9999), maximum = TRUE)
res$maximum
}
est_mv <- estimador_mv(datos_error)
est_mv[1] 0.2499972Y obtenemos el resultado esperado. En este ejemplo particular, no es necesario optimizar, puede obtenerse analíticamente el valor del estimador.
Cuando calculamos la verosimilitud arriba, nótese que estamos multiplicando números que pueden ser muy chicos (por ejemplo \(p^6\), etc). Esto puede producir desbordes numéricos fácilmente. Por ejemplo para un tamaño de muestra de 1000, podríamos tener que calcular
p <- 0.1
proba <- (p ^ 800)*(1-p)^200
proba[1] 0En estos casos, es mejor hacer los cálculos en la escala logarítmica. El logaritmo convierte productos en sumas, y baja exponentes multiplicando. Si calculamos en escala logaritmica la cantidad de arriba, no tenemos problema:
log_proba <- 800 * log(p) + 200 * log(1-p)
log_proba[1] -1863.14Ahora notemos que
Maximizar la verosimilitud es lo mismo que maximizar la log-verosimilitud, pues el logaritmo es una función creciente. Si \(x_{max}\) es el máximo de \(f\), tenemos que \(f(x_{max})>f(x)\) para cualquier \(x\), entonces tomando logaritmo, \[log(f(x_{max}))>log(f(x))\] para cualquier \(x\), pues el logaritmo respeta la desigualdad por ser creciente.
Usualmente usamos la log verosimilitud para encontrar estimador de máxima verosimilitud
Hay razónes teóricas y de interpretación por las que también es conveniente hacer esto.
En 2020 se hizo un estudio de seroprevalencia de COVID en Santa Clara, California.
Se tomó una muestra de individuos (ver los detalles en el artículo original). Para propósitos de este análisis supondremos que la muestra puede considerarse como aleatoria simple.
Se obtuvieron 3,300 individuos, y 50 de ellos resultaron con prueba positiva (1.5%). Sin embargo, el kit de prueba que estaban utilizando tenía una especificidad de 99.5% (probabilidad de dar negativo para una persona que no está infectada) y una sensibilidad de 83%. Estas cantidades también fueron con estimadas con muestras gold standard no muy grandes, y reportan un intervalo de 99.2% a 99.7% para la primera cantidad y de 76% a 88.4% para la segunda. Por el momento usaremos las estimaciones puntuales.
Construiremos una estimación de máxima verosimilitud para el % de personas \(\theta\) que ha tenido COVID en Santa Clara al momento del estudio.
El parámetro de interés principal \(p\), la proporción de personas seroprevalentes.
La probabilidad de observar un resultado positivo es, por probabilidad total:
\[ p\theta_{pos} + (1-p)(1-\theta_{neg}),\]
es decir, podemos observar un verdadero positivo o un falso positivo. Escribimos la función de (log) verosimilitud:
crear_log_verosim <- function(n, positivos, theta_pos, theta_neg){
log_verosim <- function(p){
prob_pos <- theta_pos * p + (1 - theta_neg) * (1-p)
res <- positivos * log(prob_pos) + (n - positivos) * log(1 - prob_pos)
res
}
log_verosim
}
log_verosim <- crear_log_verosim(3300, 50, 0.83, 0.995)
log_verosim(0.10)[1] -419.3992Y ahora maximizamos y reportamos el estimador de máxima verosimilitud en porcentaje:
estimador_mv <- function(n, positivos, theta_pos, theta_neg){
log_verosim <- crear_log_verosim(n, positivos, theta_pos, theta_neg)
res <- optim(par = 0.2,
log_verosim,
method = "L-BFGS-B",
lower = c(0.00001),
upper = c(0.99999),
control = list(fnscale = -1))
res$par
}
est_mv <- estimador_mv(3300, 50, 0.83, 0.995)
100 * est_mv[1] 1.231146Como vemos, el resultado tiene un ajuste con respecto a la muestra debido a las características del kit de prueba (especialmente por los falsos positivos).
Las estimaciones de arriba son puntuales, y generalmente requerimos medidas de incertidumbre. Podemos usar una variación del bootstrap que supone que simula de nuestro modelo de probabilidad para entender cuánto puede variar nuestra estimación
Consideramos el ejemplo del estudio de Santa Clara. Ya tenemos el código para encontrar el estimador de máxima verosimilitud. Pero necesitaremos simular datos de nuestro modelo de probabilidad cuando conocemos el parámetro \(\theta\) (la tasa de seroprevalencia). En este caso, la función de simulación es simple:
simular_muestra <- function(n, theta, theta_pos, theta_neg){
verdaderos_pos <- rbernoulli(n, theta)
positivos <- map_int(verdaderos_pos, function(valor){
ifelse(valor, rbernoulli(1, theta_pos), rbernoulli(1, 1 - theta_neg))
})
sum(positivos)
}
simular_muestra(3300, 0.01, 0.83, 0.995)[1] 44Y ahora podemos escribir una iteración del bootstrap paramétrico:
simular_boot <- function(est_mv, theta_pos, theta_neg){
# simular muestra bootstrap
muestra_boot <- simular_muestra(3300, est_mv, theta_pos, theta_neg)
# calcular nuevo estimador max verosim
est_boot_mv <- estimador_mv(3300, muestra_boot, theta_pos, theta_neg)
est_boot_mv
}
simular_boot(est_mv, 0.83, 0.995)[1] 0.01231146Repetimos un número grande de veces (de mil a 10 mil por ejemplo):
sims_boot <- map_df(1:1000, function(rep){
tibble(rep = rep, p = simular_boot(est_mv, 0.83, 0.99))
})
ggplot(sims_boot, aes(x = p)) + geom_histogram()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
quantile(sims_boot$p, c(0.05, 0.95)) |> round(3) 5% 95%
0.008 0.017 Y este sería nuestro intervalo de confianza para la prevalencia bajo los supuestos que hicimos. Sin embargo, también tenemos que considerar la incertidumbre en las cantidades de sensibilidad y especificidad como veremos a continuación.
Podemos probar qué pasaria si la especificidad estuviera de lado bajo según los intervalos reportados:
est_bajo_mv <- estimador_mv(3300, 50, 0.83, 0.99)
100 * est_bajo_mv[1] 0.6300095que es considerablemente más bajo que el número anterior. Nótese que parece ser que este análisis es muy dependiente de la verdadera especificidad. Por ejemplo, si hacemos la especificidad 98%, obtenemos que el estimador de máxima verosimilitud es casi cero, es decir: lo más probable es que estemos observando solamente falsos positivos.
100 * estimador_mv(3300, 50, 0.83, 0.98)[1] 0.001Veremos cómo tratar esto de otra manera en la sección de inferencia bayesiana. Sin embargo, consideremos que pasa con nuestra estimación si la especificidad está del lado bajo:
sims_boot <- map_df(1:2000, function(rep){
tibble(rep = rep, p = simular_boot(est_bajo_mv, 0.83, 0.99))
})
ggplot(sims_boot, aes(x = p)) + geom_histogram(binwidth = 0.0005)
quantile(sims_boot$p, c(0.05, 0.95)) |> round(3) 5% 95%
0.002 0.011 Y el resultado que obtuvimos es consistente con una seroprevalencia de hasta 0.2%. Claramente, las conclusiones acerca del subreporte de resultados y de el riesgo de COVID cambia mucho del análisis anterior a este que considera el extremo inferior de especificidad del kit de prueba.