From 7e183d18c11c3615ed9afe06482cb3682f6989b6 Mon Sep 17 00:00:00 2001 From: Felipe Gonzalez Date: Tue, 16 Jan 2024 13:00:41 -0600 Subject: [PATCH] Separar secciones, y algunos agregados --- notas/02-flujo-basico-2.qmd | 430 +++++++++++++++++++++++++ notas/02-flujo-basico.qmd | 608 ++++++++++-------------------------- notas/_quarto.yml | 3 +- 3 files changed, 590 insertions(+), 451 deletions(-) create mode 100644 notas/02-flujo-basico-2.qmd diff --git a/notas/02-flujo-basico-2.qmd b/notas/02-flujo-basico-2.qmd new file mode 100644 index 0000000..e192be9 --- /dev/null +++ b/notas/02-flujo-basico-2.qmd @@ -0,0 +1,430 @@ +# Flujo de trabajo básico: refinando el modelo + +```{r} +#| include: false +library(tidyverse) +library(kableExtra) +library(DiagrammeR) +library(rethinking) +ggplot2::theme_set(ggplot2::theme_light()) +``` + +En esta parte veremos cómo continuaríamos refinando el modelo tomando +en cuenta otros aspectos importantes del problema de estimar la seropositividad +de una población. + +El primer aspecto importante es persistir **seguiendo el flujo de trabajo**. En +segundo lugar, este ejemplo ilustra que la construcción de modelos preferiblemente +se hace de manera iterativa: una vez que probamos nuestro código, entendemos cómo +funciona nuesto modelo, podemos continuar agregando componentes para hacerlo +más útil. + +## Prevalencia con error conocido + + +Nuestro ejemplo de la sección es poco realista. Usualmente, las pruebas que +son utilizadas para medir la prevalencia no son perfectas. Bajo condiciones +muy controladas, el perfil de desempeño de las pruebas se miden, y los resultados +son del siguiente tipo: + +- En pruebas de *gold standard*, el kit identificó correctamente como positivos a 103 de 122 personas infectadas, e identificó correctamente como negativos a 399 de 401 personas no infectadas. +- Sin considerar la incertidumbre, esto implica que la prueba tiene una sensibilidad de 84% y una especificidad de 99.5%. + +### Paso 1: modelo generativo + +Primero supondremos que estos porcentajes de error son fijos. Nuestro modelo +que incluye el error de medición se como sigue: + + +```{r} +#| out-width: 100% +#| code-fold: true +grViz(" +digraph { + graph [ranksep = 0.3, rankdir = LR] + node [shape=circle] + p + Npos + node [shape=plaintext] + N + Npos [label = +>] + Nobs + #Nneg [label = ->] + #sens + #esp + edge [minlen = 3] + p -> Npos + #p -> Nneg + N -> Npos + Npos -> Nobs + #N -> Nneg + esp -> Nobs + sens -> Nobs + #esp -> Nneg + #sens -> Nneg +{ rank = same; p; N } +{ rank = same; Npos} +{ rank = max; sens; esp} +} +")#, width = 200, height = 50) +``` + +Donde vemos ahora que el estado real de cada persona de la prueba +es desconocido, aunque el resultado de la prueba depende de ese estado, y +la cantidad de positivos que observamos es ahora $N_{obs}$, que depende también +de la sensibilidad y especificidad de la prueba. + +Y para constuir el modelo generativo notamos que +la probabilidad de que un individuo infectado salga positivo +es $\text{sens}$, y la probabilidad de que un individuo no infectado +salga positivo es $(1-\text{esp})$. De este modo, el modelo generativo +es: + + +```{r} +sim_pos_neg <- function(theta = 0.01, N = 20, sens = 0.84, esp = 0.995) { + # verdaderos positivos que capturamos en la muestra + Pos_verdadero <- rbinom(N, 1, theta) + Neg_verdadero <- 1 - Pos_verdadero + # positivos observados en la muestra: si es positivo, calculamos + # la probabilidad de que realmente sea positivo + sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> + mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> + mutate(Neg = 1 - Pos) + # Observaciones + sim_tbl |> select(Pos, Neg) +} +``` + +Hacemos unas pruebas: + +```{r} +set.seed(8212) +sim_pos_neg(theta = 0.3, N = 1e7, sens = 0.7, esp = 1) |> pull(Pos) |> mean() |> + round(4) +sim_pos_neg(theta = 0.3, N = 1e7, sens = 1, esp = 1) |> pull(Pos) |> mean() |> + round(4) +``` + +### Paso 2: cantidad a estimar + +En este punto hay que tener cuidado, porque *no* queremos estimar la proporción +de positivos potenciales en la población (pues la prueba es imperfedta), sino +la proporción de verdaderos +positivos en la población. Esta cantidad sigue siendo representada +por $\theta$ en nuestro modelo generativo. + + +### Paso 3: modelo estadístico + +El modelo estadístico es ahora diferente. Vamos a plantear primero +$p(D|\theta, sens, esp)$, que es la probabilidad de observar los datos $D$ dado +que $\theta$ es el parámetro de interés, y $sens$ y $esp$ (que en este caso +suponemos conocidos). Es fácil ver que la probabilidad de obtener +un positivo ahora es: + +$\theta_{obs} = P(Positivo | \theta, sens, esp) = \theta \cdot sens + (1-\theta) \cdot (1-esp)$ + +Si llamamos a esta cantidad $\theta_{obs}$, de forma que dada una muestra +de 0's y 1's, tenemos + +$$p(D|\theta, sens, esp) = \theta_{obs}^{N_{+}}(1-\theta_{obs})^{N_{-}}$$ +Suponiendo que la distribución apriori de $\theta$ es uniforme, tenemos entonces +que la distribución posterior cumple: + +$$p(\theta|D, sens, esp) \propto \theta_{obs}^{N_{+}}(1-\theta_{obs})^{N_{-}}$$ +donde $\theta_obs$ está dada por la fórmula de arriba. Sustituyendo: + +$$p(\theta|D, sens, esp) \propto (\theta \cdot sens + (1-\theta) \cdot (1-esp))^{N_{+}}(\theta(1-sens) + (1-\theta)esp)^{N_{-}}$$ + +Esta posterior tiene la estructura de una distribución beta, pero es +un poco más complicada. En este punto, utilizaremos una técnica que funciona +para problemas chicos (de unos cuantos parámetros), y que consiste en +hacer una aproximación discreta de la distribución posterior: + + +::: callout-note +# Método de aproximación de rejilla + +1. Dividimos el intervalo $[0,1]$ en $m$ partes iguales, y calculamos +el valor de la expresión proporcional a la posterior en cada uno de estos +2. Normalizamos estos valores para que sumen 1, y obtenemos una distribución discreta +que aproxima la posterior +3. Muestreamos de esta distribución discreta para obtener una muestra de la posterior + +Este método sólo es factible en modelos simples +cuando hay solamente unos cuantos parámetros por estimar, +pues su complejidad crece exponencialmente con el número de parámetros. Rara vez +se usa en la práctica por esta razón. +::: + + +Aquí implementamos esta técnica de aproximación por rejilla. Incluimos +también una Beta(1,3) como a priori: + +```{r} +simular_posterior_error <- function(muestra, n, sens = 1, esp = 1){ + theta <- seq(1e-12, 1-1e-12, by = 0.0001) + p_obs <- theta * sens + (1 - theta) * (1 - esp) + # verosimilitud (en logaritmo) + log_dens_sin_norm <- log(p_obs) * sum(muestra) + + log(1-p_obs) * (length(muestra) - sum(muestra)) + # a priori + log_dens_sin_norm <- log_dens_sin_norm + dbeta(theta, 1, 3, log = TRUE) + # normalizar + log_dens_norm <- log_dens_sin_norm - log_sum_exp(log_dens_sin_norm) + densidad_post <- exp(log_dens_norm) + tibble(theta = sample(theta, size = n, replace = TRUE, prob = densidad_post)) +} +``` + +Y ahora podemos ver cómo se ve la posterior: + +```{r} +set.seed(328) +una_muestra <- sim_pos_neg(theta = 0.2, N = 600, sens = 0.6, esp = 0.999) +mean(una_muestra$Pos) +sims_post_error <- + simular_posterior_error(una_muestra$Pos, 5000, sens = 0.6, esp = 0.999) +sims_post_error |> + ggplot(aes(x = theta)) + + geom_histogram() +``` + +Ahora seguimos el flujo. Agregaremos la verificación a priori para entender +si nuestro modelo recupera los parámetros. + +```{r} +set.seed(8112) +simulacion_rep_error <- map_df(1:20, + function(rep){ + # simular de la apriori + theta_sim <- rbeta(1, 1, 3) + # simular datos según modelo + datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999) + # simulaciones montecarlo para la posterior + posterior <- simular_posterior_error(datos_sim$Pos, 10000, sens = 0.6, esp = 0.999) + # junta todo + posterior |> mutate(n_sim = n()) |> + mutate(rep = rep) |> + mutate(theta_sim = theta_sim) + }) +``` + +Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, +y comparamos con los valores de la simulación: + +```{r} +ggplot(simulacion_rep_error, aes(x = theta)) + + geom_histogram(bins = 50) + + labs(x = "theta", y = "Prob posterior") + + geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + + facet_wrap(~rep) +``` + +Contrasta con lo que pasaría si usaramos el modelo sin considerar fuentes de error: + +```{r} +set.seed(812) +simulacion_rep <- map_df(1:20, + function(rep){ + # simular de la apriori + theta_sim <- rbeta(1, 1, 3) + # simular datos según modelo + datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999) + # simulaciones montecarlo para la posterior + posterior <- simular_posterior_error(datos_sim$Pos, 10000, 1, 1) + # junta todo + posterior |> mutate(n_sim = n()) |> + mutate(rep = rep) |> + mutate(theta_sim = theta_sim) + }) +``` + +Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, +y comparamos con los valores de la simulación: + +```{r} +ggplot(simulacion_rep, aes(x = theta)) + + geom_histogram(bins = 50) + + labs(x = "theta", y = "Prob posterior") + + geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + + facet_wrap(~rep) +``` +Este resultado está lejos de ser aceptable. + +Comparamos esta densidad con lo que obtendríamos sin considerar +el error de medición, con los mismos datos: + +```{r} +set.seed(8) +sims_post <- + simular_posterior_error(una_muestra$Pos, 5000, 1, 1) +ambas_sims_tbl <- + sims_post_error |> + mutate(tipo = "Con error") |> + bind_rows(sims_post |> + mutate(tipo = "Sin error")) +ambas_sims_tbl |> ggplot(aes(x = theta, fill = tipo)) + + geom_histogram(position = "identity", alpha = 0.5, bins = 50) + + scale_fill_manual(values = c("red", "blue")) + + geom_vline(xintercept = 0.2, linetype = "dashed", color = "black") + +``` +Y vemos que la diferencia entre las distribuciones es considerable. En primer +lugar, la distribución con error de medición es más ancha (hay más incertidumbre). +En segundo lugar, como estimador de el parámetro de interés, nuestro +modelo que no considera el error parece dar estimaciones sesgadas hacia +abajo. Esto es porque la prevalencia no es tan baja, y la sensibilidad de la +prueba no es muy buena, de manera que con el modelo con error inferimos +correctamente que hay más prevalencia que lo que indicaría la proporción +de positivos en las pruebas. + + + +Aunque este ejemplo es claro, prevalencia, sensibilidad y especificidad interactúan de maneras +a veces poco intuitivas. + +## Prevalencia con datos de referencia + +Ahora haremos un paso adicional: los valores de sensibilidad y especificidad +generalmente no son conocidos con certeza, sino que son estimados a partir de +una muestra de "estándar de oro". En esta prueba particular, el kit identificó correctamente como positivos a 103 de 122 personas infectadas, +e identificó correctamente como negativos a 399 de 401 personas no infectadas. +Consideraremos 122 y 401 como tamaños de muestra fijos y conocidos (las personas +fueron extraídas de otra población). + +Denotamos como $Ref$ a los datos de referencia de "estándar de oro". + + +```{r} +#| out-width: 100% +#| code-fold: true +grViz(" +digraph { + graph [ranksep = 0.3, rankdir = LR] + node [shape=circle] + theta + esp + sens + Npos [label = +>] + node [shape=plaintext] + Nobs [label = obs>] + # Nneg [label = ->] + edge [minlen = 3] + theta -> Npos + #p -> Nneg + N -> Npos + Npos -> Nobs + #N -> Nneg + esp -> Nobs + sens -> Nobs + #esp -> Nneg + #sens -> Nneg + esp -> Ref + sens -> Ref +{ rank = same; theta; N } +#{ rank = same; Npos; Nneg} +{ rank = max; sens; esp} +} +")#, width = 200, height = 50) +``` + +Por los argumentos de arriba, las distribuciones de esp y sens son beta +y podemos incorporarlas en la simulación de la posterior. Nuestra nueva +función para simular el proceso generativo es: + +```{r} +sim_pos_neg <- function(p = 0.01, N = 20, pos_gold = c(103,122), neg_gold = c(399,401)) { + # Simular especificidad y sensibilidad + sens <- rbeta(1, pos_gold[1] + 1, pos_gold[2] - pos_gold[1] + 1) + esp <- rbeta(1, neg_gold[1] + 1, neg_gold[2] - neg_gold[1] + 1) + # verdaderos positivos que capturamos en la muestra + Pos_verdadero <- rbinom(N, 1, p) + Neg_verdadero <- 1 - Pos_verdadero + # positivos observados en la muestra: si es positivo, calculamos + # la probabilidad de que realmente sea positivo + sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> + mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> + mutate(Neg = 1 - Pos) + # Observaciones + sim_tbl |> select(Pos, Neg) +} +``` + +Considerando que tenemos tres parámetros, en este punto decidimos no +hacer la aproximación de rejilla. Es posible hacer otro tipo de aproximaciones +(por ejemplo cuadráticas), pero en lugar de esto veremos cómo lo haríamos +con [Stan](https://mc-stan.org/). Más adelante discutiremos los algoritmos que Stan utiliza para +simular de la posterior de modelos muy generales. Por el momento, notamos +que está basado en un algoritmo de simulación MCMC (Markov Chain Montecarlo), que es el estándar +para modelos que no son muy simples. **Este ejemplo es para ilustrar cómo +resolveríamos el problema más general**, no es necesario que en este punto +entiendas cómo funciona o los detalles de la implementación. + + +```{r} +library(cmdstanr) +mod_sc <- cmdstan_model("./src/sclara.stan") +print(mod_sc) +``` + + + +```{r} +n <- 50 +N <- 3300 +datos_lista <- list(N = 3300, n = 50, + kit_pos = 103, n_kit_pos = 122, + kit_neg = 399, n_kit_neg = 401) +ajuste <- mod_sc$sample(data = datos_lista, refresh = 1000) +sims <- ajuste$draws(c("p", "sens", "esp"), format = "df") +resumen <- ajuste$summary(c("p")) +``` + +```{r} +resumen |> select(variable, mean, q5, q95) +``` + +Y podemos graficar la posterior de la seroprevalencia: + +```{r, fig.width = 5, fig.height=3.5} +ggplot(sims, aes(x = p)) + + geom_histogram() +``` + +Y vemos que los datos son consistentes con el dato reportado por los autores (alrededor +de 1.2%), pero que no podemos excluir valores de prevalencia muy bajos (por abajo de 0.3% por ejemplo). Por otro lado, también son consistentes valores muy altos de seroprevalencia, de manera que este estudio resultó ser poco informativo de la IFR del COVID. + +Podemos hacer diagnósticos adicionales acerca de la razón de esta variabilidad alta, +si graficamos la relación entre especificidad de la prueba y estimación de prevalencia: + +```{r, fig.width = 5, fig.height=4} +ggplot(sims, aes(x = esp, y = p)) + geom_point() + + xlab("Especificidad del kit") + ylab("Prevalencia") + geom_smooth() +``` + +La asociación entre estas dos cantidades es interesante porque conceptualmente +(y desde punto de vista del modelo), no hay relación entre estas dos variables: su asociación +aparece porque son causas que compiten para explicar una observación. + +Nótese que dada la prevalencia baja, la especificidad del kit es un factor importante +para explicar la prevalencia observada, pero si no pensamos con cuidado podríamos +concluir que los falsos positivos no deberían ser problema por que la especificidad +para ser muy buena. + + + + + + + + + + + + + +Y notamos que aún con una muestra relativamente grande, +el rango de $p$ es considerable: va desde valores cercanos a 0 hasta valores +alrededor de 0.025-0.03. diff --git a/notas/02-flujo-basico.qmd b/notas/02-flujo-basico.qmd index 7e0b221..e747591 100644 --- a/notas/02-flujo-basico.qmd +++ b/notas/02-flujo-basico.qmd @@ -1,4 +1,4 @@ -# Motivación: Flujo de trabajo básico +# Flujo de trabajo básico: motivación ```{r} #| include: false @@ -72,7 +72,7 @@ digraph { ``` Que también podríamos simplificar (suponiendo la $N$ fija y conocida, pues -$N_+$ y $M$ dan $N_$) como: +$N_+$ y $M$ dan $N_{-}$) como: ```{r} @@ -125,19 +125,9 @@ sim_pos_neg <- function(theta = 0.01, N = 20, sens = 1, esp = 1) { } ``` -::: callout-tip -# Pruebas unitarias -La práctica estándar de pruebas unitarias debe utilizarse también, en -la medida de los posible, en estadística. Para evitar pruebas *flaky* podemos -usar semillas del -generador de números aleatorios, probar casos extremos, o -muestras muy grandes. -::: - -Por el momento no seremos muy rigurosos en este aspecto, pero al menos -podemos hacer algunas pruebas del modelo generativo en casos extremos: +Podemos hacer algunas pruebas del modelo generativo en casos extremos: ```{r} set.seed(8212) @@ -147,6 +137,23 @@ sim_pos_neg(theta = 0.1, N = 1e7) |> pull(Pos) |> mean() |> round(4) ``` +En la práctica podemos definir pruebas más exhaustivas si es necesario. En +este caso, se trata principalmente de *pruebas unitarias* que se utilizan +comunmente en desarrollo de software. + +::: callout-tip +# Pruebas unitarias + +La práctica estándar de pruebas unitarias consiste en probar +unidades relativamente pequeñas de código (por ejemplo funciones) +para verificar que funcionan correctamente. + +Esta estrategia debe utilizarse también, en +la medida de los posible, en estadística. + +::: + + ## Paso 2: Definir estimando Ahora podemos definir en términos de nuestro modelo el valor que queremos estimar. @@ -155,7 +162,7 @@ es así siempre: como veremos más adelante, puede ser una cantidad que se deriv de otras variables y parámetros del modelo. -## Introdcción a Paso 3: definir un proceso estadístico +## Paso 3: definir un proceso estadístico Dada la información limitada que tenemos acerca de la población, esperamos tener cierta incertidumbre en nuestra estimación del valor de $\theta$. En estadística @@ -164,7 +171,8 @@ sobre posibles valores del $\theta$. Si denotamos por $D$ a los datos observados nuestro objetivo es calcular o aproximar $$p(\theta|D)$$ -que es una distribución sobre los posibles valores de $\theta$, y que pone más +que es una distribución sobre los posibles valores de $\theta$, una vez +que tenemos información de la muestra, y que pone más masa de probabilidad sobre las conjeturas de $\theta$ que son más probables o creíbles. A esta distribución le llamamos la **distribución posterior** de $\theta$. @@ -173,15 +181,15 @@ Con esta posterior podemos hacer afirmaciones probabilísticas de la forma: - ¿Cuál es la probabilidad de que $\theta$ sea menor a 1%? (Muy pocos seropositivos) - ¿Cuál es la probabildad de que $\theta$ sea mayor a 80%? (Población cerca de saturación) -que se calculan integrando $p(\theta|D)$ sobre los valores de $\theta$ que nos interesan, por ejemplo, +Estas cantidades se calculan, al menos teóricamente, +integrando $p(\theta|D)$ sobre los valores de $\theta$ que nos interesan, por ejemplo, $$P(\theta <= 0.01) = \int_0^{0.01} p(\theta|D) d\theta$$ Nota: la integral la interpretamos como suma en el caso discreto. -### Análisis de datos bayesiano: motivación - Supongamos entonces una $\theta$ dada, y que observamos la muestra -$1,0,0,1,0$. La probabilidad de observar esta muestra es entonces +$1,0,0,1,0$. La probabilidad de observar esta muestra es (suponiendo observaciones +independientes): $$\theta(1-\theta)(1-\theta)\theta(1-\theta) = \theta^2(1-\theta)^3$$ Para algunos valores de $\theta$ (posibles conjeturas acerca del valor de $\theta$) podemos @@ -205,20 +213,24 @@ $$p(D|\theta)$$ significa: la probabilidad de los datos $D$ dado el valor de $\theta$. Nótese que esta distribución **no** es la posterior que describimos arriba, y no es una distribución de probabilidad sobre $\theta$ (las probabilidades no suman uno). +Esta función se llama usualmente **verosimilitud** de los datos, e incorpora +supuestos concretos del proceso generador de los datos. Usando reglas de probabilidad (en particular la regla de Bayes), observamos que -$$p(\theta | D) = \frac{p(D|\theta)p(\theta)} { p(D)}$$ +$$p(\theta | D) = \frac{p(D|\theta)p(\theta)} { p(D)}.$$ Como $p(\theta|D)$ debe dar una distribución de probabilidad (suma o integra a 1), entonces $p(D)$ debe ser una constante de normalización para el numerador de la derecha, es decir, basta escribir $$p(\theta | D) \propto p(D|\theta)p(\theta) $$ Ahora es donde encontramos que tenemos que tener $p(\theta)$ para poder calcular -la cantidad que nos interesa. $p(\theta)$ es simplemente una afirmación de dónde +la cantidad que nos interesa, que es la distribución posterior $p(\theta|D)$. $p(\theta)$, +la **distribución a priori** o **distribución inicial** +es simplemente una afirmación de dónde puede estar $\theta$, antes de observar ningún dato. -Por el momento, podríamos quizá evitar este problema poniendo $p(\theta)$ constante, +Por el momento, podríamos poner $p(\theta)$ constante, de manera que es parte de la constante de normalización, y sólo tendríamos que normalizar como sigue: ```{r} @@ -232,34 +244,39 @@ prob_post |> Con esto, expresamos nuestro conocimiento acerca de $\theta$, después de observar los -datos, con una *distribución de probabilidad sobre las posibles conjecturas*. A esta -distribución de probabilidad le llamamos **probabilidad posterior**. Este es el +datos, con una *distribución posterior de probabilidad sobre las posibles conjecturas*. Este es el resultado principal de inferencia bayesiana, y es la base para tomar decisiones -relativas a $p$. +relativas a $\theta$. -### Usando información adicional +### Usando información adicional {-} Supongamos que tenemos información adicional acerca de $\theta$, por ejemplo, que en un experimento similar anterior alguien tomó una muestra de dos personas, -y encontró dos negativos. Tenemos entonces como creencias inciales: +y encontraron dos negativos. Tenemos entonces como creencias inciales: ```{r} -p <- seq(0, 1, length.out = 11) -prob_priori <- tibble(conjetura = p, prob_priori = (1-p) * (1 -p)) |> +theta <- seq(0, 1, length.out = 11) +prob_priori <- tibble(conjetura = theta) |> + mutate(prob_priori = (1 - theta) * (1 - theta)) |> mutate(prob_priori = prob_priori / sum(prob_priori)) prob_priori |> kable(col.names = c("Conjetura θ", "p(θ)")) |> kable_paper() ``` +Por ejemplo, al probabilidad inicial de que $\theta$ sea muy grande es cercana a cero, +pues observamos dos negativos y ningún positivo. Ahora regresamos a considerar nuestra +fórmula -Ahora reconsideramos la posterior de nuestra muestra de 5 personas, y calculamos +$$p(\theta | D) \propto p(D|\theta)p(\theta), $$ + +En este caso, la apriori o inicial tiene un efecto sobre la posterior. +Reconsideramos entonces la posterior de nuestra muestra de 5 personas, y calculamos el producto de $P(D|\theta)$ por $p(\theta)$: ```{r} prob_post <- prob_priori |> - mutate(prob_post_ant = p^2 * (1-p)^3) |> - mutate(prob_post_ant = prob_post_ant / sum(prob_post_ant)) |> - mutate(prod = prob_priori * prob_post_ant) + mutate(verosimilitud = theta^2 * (1 - theta)^3) |> + mutate(prod = verosimilitud * prob_priori) prob_post|> kable(col.names = c("Conjetura θ", "p(θ)", "p(D|θ)", @@ -283,7 +300,7 @@ La última columna nos da el resultado final de la inferencia bayesiana. Podemos resumir algunas de sus características, por ejemplo: - Es muy poco probable que la seropositividad sea mayor o igual a 0.7 -- Un intervalo de 98% de probabilidad para la seropositividad es $[0.1, 0.5]$ +- Un intervalo de 90% de probabilidad para la seropositividad es $[0.1, 0.5]$ La gráfica de la posterior es: @@ -295,8 +312,6 @@ prob_post |> ``` -## Paso 3: definir un proceso estadístico - Ahora podemos definir, para nuestro ejemplo discretizado, la función que calcula la posterior dados los pasos 1 y 2: @@ -344,21 +359,21 @@ de acuerdo a los supuestos de nuestro modelo generativo. Lo mínimo que esperamos de nuestro método es que, bajo nuestros propios supuestos acerca del proceso generador de datos y nuestro procedimiento de estimación -definido, nuestra función de estimación: - -- No tenga problemas numéricos o de programación -- Las estimaciones que arroja son apropiadas para la cantidad que nos interesa - -::: - -El segundo punto es más difícil de definir, y más adelante lo consideraremos con -detalle. El procedimiento general es: +definido, nuestra función de estimación no tenga problemas numéricos o de programación, +y que las estimaciones que arroja son apropiadas para la cantidad que nos interesa +estimar. El procedimiento a grandes rasgos es: 1. Establecer valores de los parámetros a estimar 2. Simular datos observados (con una $N$ apropiada, dependiendo del tamaño de muestra que esperamos, aunque se puede explorar hacer más grande o más chico este valor). 3. Calcular posterior de las cantidades de interés 4. Compara los valores de 1) con la posterior de 3) +::: + +Definir que las posteriores son apropiadas para la cantidad que nos interesa +estimar es delicado, y más adelante veremos algunos criterios para evaluar este +aspecto. Por lo pronto, haremos algunas pruebas simples que pueden diagnosticar +errores graves: ```{r} theta <- 0.2 @@ -376,6 +391,7 @@ En este caso, la estimación parece correcta. Podemo repetir el proceso con distintos valores de $\theta$: ```{r} +set.seed(21) simulacion_rep <- map_df(1:20, function(rep){ # simular valor inicial @@ -393,7 +409,7 @@ ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) + geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + facet_wrap(~rep) ``` -Y vemos que en general nuestro método parece funcionar correctamente +Y vemos que en general nuestro método parece funcionar correctamente. ::: callout-tip # Observaciones @@ -408,6 +424,11 @@ prueba es en general más apropiada, pues no nos interesan configuración de par con probabilidad inicial extremadamente baja (imposibles según nuestros supuestos), pero también es posible tomar algunos valores fijos de interés. +- Veremos más adelante también *pruebas predictivas a priori*, que son una +manera de checar supuestos, y también +nos sirven para evaluar si nuestras distribuciones apriori son consistentes +con conocimiento previo o teoría. + ::: Este paso también es importante para entender si, bajo nuestros propios supuestos, @@ -436,6 +457,97 @@ Nuestra respuesta en este caso es que quizá con 3 personas la información obte no será suficiente para tomar decisiones útiles: nótese que la posterior está muy poco concentrada alrededor del verdadero valor de $\theta$. +### Introduciendo un bug + +Supongamos que tenemos un error en el cálculo de la posterior: + +```{r} +calcular_posterior_bug <- function(muestra, prob_priori){ + # distribución inicial o a prior + theta <- seq(0, 1, length.out = 11) + priori <- tibble(theta = theta, prob_priori = (1 - theta) * (1 - theta)) |> + mutate(prob_priori = prob_priori / sum(prob_priori)) + # calcular la probabilidad posterior + N <- length(muestra) + Npos <- sum(muestra) + prob_post <- tibble(theta = theta) |> + left_join(priori, by = "theta") |> + # la siguiente línea tiene un error! + mutate(prob_posterior = theta ^ Npos * (1 - theta)^((N - Npos * prob_priori))) |> + mutate(prob_posterior = prob_posterior / sum(prob_posterior)) + prob_post |> select(theta, prob_posterior) +} +``` + +Nuestro chequeo apriori se ve entonces: + +```{r} +simulacion_rep <- map_df(1:20, + function(rep){ + # simular valor inicial + theta_sim <- sample(seq(0, 1, length.out = 11), + prob = prob_priori$prob_priori, size = 1) + datos_sim <- sim_pos_neg(theta = theta_sim, N = 30) + posterior <- calcular_posterior_bug(datos_sim$Pos) + posterior |> mutate(theta = theta) |> + mutate(rep = rep) |> + mutate(theta_sim = theta_sim) + }) +ggplot(simulacion_rep, aes(x = theta, y = prob_posterior)) + + geom_col() + + labs(x = "theta", y = "Prob posterior") + + geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + + facet_wrap(~rep) +``` +Donde vemos en varios casos que la "posterior" está lejos de ser consistente +con los valores simulados de prueba para $\theta$. + +::: callout-tip +# Aspectos numéricos + +Es importante notar que los cálculos que hicimos arriba ingoran un principio +importante al hacer cálculos de productos de probabilidades: generalmente es +mejor utilizar la **escala logarítmica** para hacer los cálculos, y sólo al +final convertir a probabilidades. Esto es porque es fácil tener +subflujos numéricos al multiplicar muchas probabilidades pequeñas. + +::: + + +Aunque en este caso no es crítico, la siguiente función sigue esta práctica +que en general es necesario seguir: + +```{r} +# Evitar desbordes al sumar exponenciales +log_sum_exp <- function(x){ + max_x <- max(x) + max_x + log(sum(exp(x - max_x))) +} +calcular_posterior <- function(muestra, prob_priori){ + # evitar 0 o 1 exactos + theta <- seq(1e-12, 1 - 1e-12, length.out = 11) + # no es necesario normalizar esta distribución apriori + log_priori <- tibble(theta = theta, log_prob_priori = 2 * log(1 - theta)) + # calcular la probabilidad posterior + N <- length(muestra) + Npos <- sum(muestra) + prob_post_tbl <- tibble(theta = theta) |> + left_join(log_priori, by = "theta") |> + # log verosimilitud + mutate(log_prob_posterior = + Npos * log(theta) + log(1 - theta) * (N - Npos)) |> + # sumar log apriori + mutate(log_prob_posterior = log_prob_posterior + log_prob_priori) |> + mutate(log_prob_posterior_norm = + log_prob_posterior - log_sum_exp(log_prob_posterior)) |> + mutate(prob_posterior = exp(log_prob_posterior_norm)) + prob_post_tbl |> select(theta, prob_posterior) +} +``` + +Ejercicio: corre las pruebas para esta versión de la función como hicimos arriba. +Este es un cambio correcto, y desde el punto de vista de desarrollo, +si nuestra batería de pruebas es apropiado podemos hacerlo con más confianza. ## Paso 5: Analizar los datos y resumir resultados. @@ -737,407 +849,3 @@ arriba están resueltos, para tener confianza en nuestras conclusiones. -## Modelo generativo: errores de medición - - -Nuestro ejemplo de arriba es poco realista. Usualmente, las pruebas que -son utilizadas para medir la prevalencia no son perfectas. Bajo condiciones -muy controladas, el perfil de desempeño de las pruebas se miden, y los resultados -son del siguiente tipo: - -- En pruebas de *gold standard*, el kit identificó correctamente como positivos a 103 de 122 personas infectadas, e identificó correctamente como negativos a 399 de 401 personas no infectadas. -- Sin considerar la incertidumbre, esto implica que la prueba tiene una sensibilidad de 84% y una especificidad de 99.5%. - -Primero supondremos que estos porcentajes de error son fijos. Nuestro modelo -que incluye el error de medición se como sigue: - - -```{r} -#| out-width: 100% -#| code-fold: true -grViz(" -digraph { - graph [ranksep = 0.3, rankdir = LR] - node [shape=circle] - p - Npos - node [shape=plaintext] - N - Npos [label = +>] - Nobs - #Nneg [label = ->] - #sens - #esp - edge [minlen = 3] - p -> Npos - #p -> Nneg - N -> Npos - Npos -> Nobs - #N -> Nneg - esp -> Nobs - sens -> Nobs - #esp -> Nneg - #sens -> Nneg -{ rank = same; p; N } -{ rank = same; Npos} -{ rank = max; sens; esp} -} -")#, width = 200, height = 50) -``` - -Donde vemos ahora que el estado real de cada persona de la prueba -es desconocido, aunque el resultado de la prueba depende de ese estado, y -la cantidad que observamos es ahora $N_{obs}$, que depende también -de la sensibilidad y especificidad de la prueba. - -Y para constuir el modelo generativo notamos que -la probabilidad de que un individuo infectado salga positivo -es $\text{sens}$, y la probabilidad de que un individuo no infectado -salga positivo es $(1-\text{esp})$. De este modo, el modelo generativo -es: - - -```{r} -sim_pos_neg <- function(theta = 0.01, N = 20, sens = 0.84, esp = 0.995) { - # verdaderos positivos que capturamos en la muestra - Pos_verdadero <- rbinom(N, 1, theta) - Neg_verdadero <- 1 - Pos_verdadero - # positivos observados en la muestra: si es positivo, calculamos - # la probabilidad de que realmente sea positivo - sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> - mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> - mutate(Neg = 1 - Pos) - # Observaciones - sim_tbl |> select(Pos, Neg) -} -``` - -Hacemos unas pruebas: - -```{r} -set.seed(8212) -sim_pos_neg(theta = 0.3, N = 1e7, sens = 0.7, esp = 1) |> pull(Pos) |> mean() |> - round(4) -sim_pos_neg(theta = 0.3, N = 1e7, sens = 1, esp = 1) |> pull(Pos) |> mean() |> - round(4) -``` - -### Paso 2: cantidad a estimar - -En este punto hay que tener cuidado, porque no queremos estimar la proporción -de positivos en la muestra (pues la muestra tiene error de medición), sino -la proporción de positivos en la población. Esta cantidad sigue siendo representada -por $\theta$ en nuestro modelo generativo. - - -### Paso 3: modelo estadístico - -El modelo estadístico es ahora diferente. Vamos a plantear primero -$p(D|\theta, sens, esp)$, que es la probabilidad de observar los datos $D$ dado -que $\theta$ es el parámetro de interés, y $sens$ y $esp$ (que en este caso -suponemos conocidos). Es fácil ver que la probabilidad de obtener -un positivo ahora es: - -$p_{obs} = P(N_{obs} = 1 | \theta, sens, esp) = \theta \cdot sens + (1-\theta) \cdot (1-esp)$ - -Si llamamos a esta cantidad $\theta_{obs}$, de forma que dada una muestra -de 0's y 1's, tenemos - -$$p(D|\theta, sens, esp) = p_{obs}^{N_{+}}(1-p_{obs})^{N_{-}}$$ -Suponiendo que la distribución apriori de $\theta$ es uniforme, tenemos entonces -que la distribución posterior cumple: - -$$p(\theta|D, sens, esp) \propto p_{obs}^{N_{+}}(1-p_{obs})^{N_{-}}$$ -para aquellos valores de $\theta$ que esté en $[0,1]$. - -$$p(\theta|D, sens, esp) \propto (\theta \cdot sens + (1-\theta) \cdot (1-esp))^{N_{+}}(\theta(1-sens) + (1-\theta)esp)^{N_{-}}$$ -Esta posterior tiene la estructura de una distribución beta, pero es -un poco más complicada. En este punto, utilizaremos una técnica que funciona -para problemas chicos (de unos cuantos parámetros), y que consiste en -hacer una aproximación discreta de la distribución posterior: - - -::: callout-note -# Método de aproximación de rejilla - -1. Dividimos el intervalo $[0,1]$ en $m$ partes iguales, y calculamos -el valor de la expresión proporcional a la posterior en cada uno de estos -2. Normalizamos estos valores para que sumen 1, y obtenemos una distribución discreta -que aproxima la posterior -3. Muestreamos de esta distribución discreta para obtener una muestra de la posterior - -Este método sólo es factible en modelos simples -cuando hay solamente unos cuantos parámetros por estimar, -pues su complejidad crece exponencialmente con el número de parámetros. Rara vez -usa en la práctica por esta razón. -::: - - -Aquí implementamos esta técnica de aproximación por rejilla. Incluimos -también una Beta(1,3) como a priori: - -```{r} -simular_posterior_error <- function(muestra, n, sens = 1, esp = 1){ - theta <- seq(10e-6, 1-10e-6, by = 0.0001) - p_obs <- theta * sens + (1 - theta) * (1 - esp) - # verosimilitud (en logaritmo) - log_dens_sin_norm <- log(p_obs) * sum(muestra) + - log(1-p_obs) * (length(muestra) - sum(muestra)) - # a priori - log_dens_sin_norm <- log_dens_sin_norm + dbeta(theta, 1, 3, log = TRUE) - # normalizar - dens_sin_norm <- exp(log_dens_sin_norm) - densidad_post <- dens_sin_norm / sum(dens_sin_norm) - tibble(theta = sample(theta, size = n, replace = TRUE, prob = dens_sin_norm)) -} -``` - -Y ahora podemos ver cómo se ve la posterior: - -```{r} -set.seed(328) -una_muestra <- sim_pos_neg(theta = 0.2, N = 600, sens = 0.6, esp = 0.999) -mean(una_muestra$Pos) -sims_post_error <- - simular_posterior_error(una_muestra$Pos, 5000, sens = 0.6, esp = 0.999) -sims_post_error |> - ggplot(aes(x = theta)) + - geom_histogram() -``` - -Ahora seguimos el flujo. Agregaremos la verificación a priori para entender -si nuestro modelo recupera los parámetros. - -```{r} -set.seed(8112) -simulacion_rep_error <- map_df(1:20, - function(rep){ - # simular de la apriori - theta_sim <- rbeta(1, 1, 3) - # simular datos según modelo - datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999) - # simulaciones montecarlo para la posterior - posterior <- simular_posterior_error(datos_sim$Pos, 10000, sens = 0.6, esp = 0.999) - # junta todo - posterior |> mutate(n_sim = n()) |> - mutate(rep = rep) |> - mutate(theta_sim = theta_sim) - }) -``` - -Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, -y comparamos con los valores de la simulación: - -```{r} -ggplot(simulacion_rep_error, aes(x = theta)) + - geom_histogram(bins = 50) + - labs(x = "theta", y = "Prob posterior") + - geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + - facet_wrap(~rep) -``` - -Contrasta con lo que pasaría si usaramos el modelo sin considerar fuentes de error: - -```{r} -set.seed(812) -simulacion_rep <- map_df(1:20, - function(rep){ - # simular de la apriori - theta_sim <- rbeta(1, 1, 3) - # simular datos según modelo - datos_sim <- sim_pos_neg(theta = theta_sim, N = 150, sens = 0.6, esp = 0.999) - # simulaciones montecarlo para la posterior - posterior <- simular_posterior(datos_sim$Pos, 10000) - # junta todo - posterior |> mutate(n_sim = n()) |> - mutate(rep = rep) |> - mutate(theta_sim = theta_sim) - }) -``` - -Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, -y comparamos con los valores de la simulación: - -```{r} -ggplot(simulacion_rep, aes(x = theta)) + - geom_histogram(bins = 50) + - labs(x = "theta", y = "Prob posterior") + - geom_vline(aes(xintercept = theta_sim), color = "red", linetype = "dashed") + - facet_wrap(~rep) -``` -Este resultado está lejos de ser aceptable. - -Comparamos esta densidad con lo que obtendríamos sin considerar -el error de medición, con los mismos datos: - -```{r} -set.seed(8) -sims_post <- - simular_posterior(una_muestra$Pos, 5000) -ambas_sims_tbl <- - sims_post_error |> - mutate(tipo = "Con error") |> - bind_rows(sims_post |> - mutate(tipo = "Sin error")) -ambas_sims_tbl |> ggplot(aes(x = theta, fill = tipo)) + - geom_histogram(position = "identity", alpha = 0.5, bins = 50) + - scale_fill_manual(values = c("red", "blue")) + - geom_vline(xintercept = 0.2, linetype = "dashed", color = "black") - -``` -Y vemos que la diferencia entre las distribuciones es considerable. En primer -lugar, la distribución con error de medición es más ancha (hay más incertidumbre). -En segundo lugar, como estimador de el parámetro de interés, nuestro -modelo que no considera el error parece dar estimaciones sesgadas hacia -abajo. Esto es porque la prevalencia no es tan baja, y la sensibilidad de la -prueba no es muy buena, de manera que con el modelo con error inferimos -correctamente que hay más prevalencia que lo que indicaría la proporción -de positivos en las pruebas. - - - -Aunque este ejemplo es claro, prevalencia, sensibilidad y especificidad interactúan de maneras -a veces poco intuitivas. - - - - -## Modelo generativo: errores de medición 2 - -Ahora haremos un paso adicional: los valores de sensibilidad y especificidad -generalmente no son conocidos con certeza, sino que son estimados a partir de -una muestra de "estándar de oro". En esta prueba particular, el kit identificó correctamente como positivos a 103 de 122 personas infectadas, -e identificó correctamente como negativos a 399 de 401 personas no infectadas. -Consideraremos 122 y 401 como tamaños de muestra fijos y conocidos (las personas -fueron extraídas de otra población). - - -```{r} -#| out-width: 100% -#| code-fold: true -grViz(" -digraph { - graph [ranksep = 0.3, rankdir = LR] - node [shape=circle] - theta - esp - sens - node [shape=plaintext] - N - Npos [label = +>] - # Nneg [label = ->] - edge [minlen = 3] - theta -> Npos - #p -> Nneg - N -> Npos - #N -> Nneg - esp -> Npos - sens -> Npos - #esp -> Nneg - #sens -> Nneg - esp -> Cal - sens -> Cal -{ rank = same; theta; N } -#{ rank = same; Npos; Nneg} -{ rank = max; sens; esp} -} -")#, width = 200, height = 50) -``` - -Por los argumentos de arriba, las distribuciones de esp y sens son beta, -y podemos incorporarlas en la simulación de la posterior. Nuestra nueva -función para simular de la posterior es: - -```{r} -sim_pos_neg <- function(p = 0.01, N = 20, pos_gold = c(103,122), neg_gold = c(399,401)) { - sens <- rbeta(1, pos_gold[1] + 1, pos_gold[2] - pos_gold[1] + 1) - esp <- rbeta(1, neg_gold[1] + 1, neg_gold[2] - neg_gold[1] + 1) - # verdaderos positivos que capturamos en la muestra - Pos_verdadero <- rbinom(N, 1, p) - Neg_verdadero <- 1 - Pos_verdadero - # positivos observados en la muestra: si es positivo, calculamos - # la probabilidad de que realmente sea positivo - sim_tbl <- tibble(Pos_verdadero, Neg_verdadero) |> - mutate(Pos = rbinom(N, 1, Pos_verdadero * sens + Neg_verdadero * (1-esp))) |> - mutate(Neg = 1 - Pos) - # Observaciones - sim_tbl |> select(Pos, Neg) -} -``` - -Considerando que tenemos tres parámetros, en este punto decidimos no -hacer la aproximación de rejilla. Es posible hacer otro tipo de aproximaciones -(por ejemplo cuadráticas), pero en este punto veremos cómo lo haríamos -con Stan. Más adelante discutiremos los algoritmos que Stan utiliza para -simular de la posterior de modelos muy generales. Por el momento, notamos -que está basado en un algoritmo de simulación MCMC (Markov Chain Montecarlo), que es el estándar -para modelos que no son muy simples. **Este ejemplo es para ilustrar cómo -resolveríamos el problema más general**, no es necesario que en este punto -entiendas cómo funciona o los detalles de la implementación. - - -```{r} -library(cmdstanr) -mod_sc <- cmdstan_model("./src/sclara.stan") -print(mod_sc) -``` - - - -```{r} -n <- 50 -N <- 3300 -datos_lista <- list(N = 3300, n = 50, - kit_pos = 103, n_kit_pos = 122, - kit_neg = 399, n_kit_neg = 401) -ajuste <- mod_sc$sample(data = datos_lista, refresh = 1000) -sims <- ajuste$draws(c("p", "sens", "esp"), format = "df") -resumen <- ajuste$summary(c("p")) -``` - -```{r} -resumen |> select(variable, mean, q5, q95) -``` - -Y podemos graficar la posterior de la seroprevalencia: - -```{r, fig.width = 5, fig.height=3.5} -ggplot(sims, aes(x = p)) + - geom_histogram() -``` - -Y vemos que los datos son consistentes con el dato reportado por los autores (alrededor -de 1.2%), pero que no podemos excluir valores de prevalencia muy bajos (por abajo de 0.3% por ejemplo). Por otro lado, también son consistentes valores muy altos de seroprevalencia, de manera que este estudio resultó ser poco informativo de la IFR del COVID. - -Podemos hacer diagnósticos adicionales acerca de la razón de esta variabilidad alta, -si graficamos la relación entre especificidad de la prueba y estimación de prevalencia: - -```{r, fig.width = 5, fig.height=4} -ggplot(sims, aes(x = esp, y = p)) + geom_point() + - xlab("Especificidad del kit") + ylab("Prevalencia") + geom_smooth() -``` - -La asociación entre estas dos cantidades es interesante porque conceptualmente -(y desde punto de vista del modelo), no hay relación entre estas dos variables: su asociación -aparece porque son causas que compiten para explicar una observación. - -Nótese que dada la prevalencia baja, la especificidad del kit es un factor importante -para explicar la prevalencia observada, pero si no pensamos con cuidado podríamos -concluir que los falsos positivos no deberían ser problema por que la especificidad -para ser muy buena. - - - - - - - - - - - - - -Y notamos que aún con una muestra relativamente grande, -el rango de $p$ es considerable: va desde valores cercanos a 0 hasta valores -alrededor de 0.025-0.03. diff --git a/notas/_quarto.yml b/notas/_quarto.yml index dfabe3a..5826360 100644 --- a/notas/_quarto.yml +++ b/notas/_quarto.yml @@ -7,7 +7,8 @@ book: chapters: - index.qmd - 01-introduccion.qmd - + - 02-flujo-basico.qmd + - 02-flujo-basico-2.qmd bibliography: referencias/book.bib format: