diff --git a/notas/02-flujo-basico.qmd b/notas/02-flujo-basico.qmd index 9bcdae2..7e0b221 100644 --- a/notas/02-flujo-basico.qmd +++ b/notas/02-flujo-basico.qmd @@ -1,4 +1,4 @@ -# Flujo de trabajo básico +# Motivación: Flujo de trabajo básico ```{r} #| include: false @@ -14,18 +14,30 @@ trabajo que seguiremos. El objetivo en este ejemplo es estimar la proporción de personas que es seropositiva de una enfermedad en una población dada, usando una muestra de la población de interés a la que se le aplicó una prueba de seropositivdad. -## Modelo generativo +Recordamos que el flujo básico es: + +1. Definir un modelo generativo para la muestra de datos. +2. Definir la cantidad que queremos estimar en relación al fenómeno de interés. +3. Definir un proceso estadístico para hacer una estimación. +4. Probar el proceso 3 usando 1 y 2. +5. (Usar datos) Analizar los datos, resumir resultados. +6. Checar cómputos y desempeño del modelo. + + + +## Paso 1: Modelo generativo Consideremos primero qué variables de interés tenemos: $p$, la proporción de seropositivos en la población, $N$ que es el número de personas a las que les hicimos la prueba, y $N_{+}$ y $N_{-}$ que cuentan el número de positivos y seronegativos -en la muestra. Supondremos que la prueba da resultados exactos. - +en la muestra. Supondremos que la prueba da resultados exactos. Denotaremos +por $\theta$ a la proporción de seropositivos en la muestra. Comenzamos construyendo el diagrama que indica cómo influye cada variable en otra (nota: no son asociaciones, sino que indican qué variables "escuchan" a otras para determinar su valor). En este caso, $N$ y $\theta$ son variable que no depende de ninguna -otra, mientras que $N_{+}$ y $N_{-}$ dependen de $N$ y $\theta$. +otra, mientras que $N_{+}$ y $N_{-}$ dependen de $N$ y $\theta$. Como $\theta$ es una +cantidad que no observamos directamente, mostramos su nodo como un círculo. ```{r} #| out-width: 100% @@ -56,10 +68,49 @@ digraph { } -")#, width = 200, height = 50) +", width = 300, height = 100) ``` -Y ahora construimos el modelo generativo: +Que también podríamos simplificar (suponiendo la $N$ fija y conocida, pues +$N_+$ y $M$ dan $N_$) como: + + +```{r} +#| out-width: 100% +#| code-fold: true +grViz(" +digraph { + graph [ranksep = 0.3, rankdir = LR] + node [shape=circle] + theta [label = <θ>] + node [shape=plaintext] + N + Npos [label = +>] + #sens + #esp + edge [minlen = 3] + theta -> Npos + N -> Npos + #esp -> Pos + #sens -> Pos + #esp -> Neg + #sens -> Neg +{ rank = same; theta; N } +{ rank = same; Npos} +#{ rank = max; sens; esp} + + +} +", width = 300, height = 100) +``` + +Y ahora construimos el modelo generativo. Supondremos que la muestra de $N$ personas +se toma de manera aleatoria de la población (una población grande, así que podemos +ignorar el efecto de muestreo). Supondremos provisionalmente, además, que la prueba +es perfecta, es decir, no hay falsos positivos o negativos. + +La siguiente función simula una muestra de $N$ personas, y regresa el número de +Positivos y Negativos en la muestra. ```{r} sim_pos_neg <- function(theta = 0.01, N = 20, sens = 1, esp = 1) { @@ -74,22 +125,60 @@ sim_pos_neg <- function(theta = 0.01, N = 20, sens = 1, esp = 1) { } ``` -Podemos en primer lugar -hacer algunas pruebas del modelo generativo en casos extremos: +::: 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: ```{r} set.seed(8212) -sim_pos_neg(theta = 0.995, N = 10) -sim_pos_neg(theta = 0.005, N = 10) +sim_pos_neg(theta = 1.0, N = 10) +sim_pos_neg(theta = 0.0, N = 10) sim_pos_neg(theta = 0.1, N = 1e7) |> pull(Pos) |> mean() |> round(4) ``` -## Análisis de datos bayesiano: motivación +## Paso 2: Definir estimando + +Ahora podemos definir en términos de nuestro modelo el valor que queremos estimar. +En este caso, coincide con un párametro del modelo $\theta$, pero no necesariamente +es así siempre: como veremos más adelante, puede ser una cantidad que se deriva +de otras variables y parámetros del modelo. + + +## Introdcción a 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 +bayesiana esta incertidumbre la expresamos mediante una distribución de probabilidades +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 +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$. + +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, + +$$P(\theta <= 0.01) = \int_0^{0.01} p(\theta|D) d\theta$$ +Nota: la integral la interpretamos como suma en el caso discreto. -En estadística bayesiana podemos consideramos las posibles explicaciones -(estados de $\theta$) de los datos, y consideramos como más creíbles aquellas explicaciones -que pueden ocurrir de manera más probable. +### 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 @@ -99,27 +188,24 @@ Para algunos valores de $\theta$ (posibles conjeturas acerca del valor de $\thet escribir una tabla como sigue (Nota: discretizamos por el momento a un número finito de valores de $p$ para hacer el argumento más simple): + + + ```{r} theta <- seq(0, 1, length.out = 11) tibble(conjetura_theta = theta, verosimiltud = theta^2 * (1 - theta)^3) |> - kable(col.names = c("Conjetura $\\theta$", "$p(D|\\theta)$"), escape = FALSE) |> kable_paper() + kbl(col.names = c("Conjetura θ", "p(D|θ)"), + escape = FALSE) ``` En la tabla vemos que hay algunas conjeturas, o posibles valores de $\theta$, que tienen probabilidad considerablemente más alta que otra. La notación $$p(D|\theta)$$ -significa: la probabilidad de los datos $D$ dado el valor de $\theta$. - -Ahora observemos que si queremos hacer **afirmaciones probabilísticas** acerca de $\theta$, -esta cantidad no es la que debemos considerar, pues no es una distribución de probabilidad. -Afirmaciones probabilísticas son afirmaciones de la forma: +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). -- ¿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) - -Para afirmaciones probabilísticas, realmente quisiéramos tener una distrubución de $\theta$ dado que observamos los datos $D$: -$$p(\theta|D)$$ Usando reglas de probabilidad (en particular la regla de Bayes), observamos que $$p(\theta | D) = \frac{p(D|\theta)p(\theta)} { p(D)}$$ @@ -132,20 +218,19 @@ 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 puede estar $\theta$, antes de observar ningún dato. -Podríamos quizá evitar este problema poniendo $p(\theta)$ constante, de manera -que sólo tendríamos que normalizar como sigue: +Por el momento, podríamos quizá evitar este problema poniendo $p(\theta)$ constante, +de manera que es parte de la constante de normalización, y sólo tendríamos que normalizar como sigue: ```{r} p <- seq(0, 1, length.out = 11) prob_post <- tibble(conjetura = p, probablidad = p^2 * (1 - p)^3) |> mutate(prob_posterior = probablidad / sum(probablidad)) prob_post |> - kable(col.names = c("Conjetura $\\theta$", "$p(D|\\theta)$", - "$p(\\theta|D)$"), escape = FALSE) |> kable_paper() + kable(col.names = c("Conjetura θ", "p(D|θ)","p(θ|D)")) |> + kable_paper() ``` - 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 @@ -163,7 +248,7 @@ p <- seq(0, 1, length.out = 11) prob_priori <- tibble(conjetura = p, prob_priori = (1-p) * (1 -p)) |> mutate(prob_priori = prob_priori / sum(prob_priori)) prob_priori |> - kable(col.names = c("Conjetura $\\theta$", "$p(\\theta)$"), escape = FALSE) |> kable_paper() + kable(col.names = c("Conjetura θ", "p(θ)")) |> kable_paper() ``` @@ -177,8 +262,8 @@ prob_post <- prob_priori |> mutate(prod = prob_priori * prob_post_ant) prob_post|> - kable(col.names = c("Conjetura $\\theta$", "$p(\\theta)$", "$p(D|\\theta)$", - "$p(D|\\theta)p(\\theta)$"), escape = FALSE) |> kable_paper() + kable(col.names = c("Conjetura θ", "p(θ)", "p(D|θ)", + "p(D|θ)p(θ)")) |> kable_paper() ``` @@ -190,9 +275,8 @@ prob_post <- prob_post |> mutate(prob_posterior = prod / sum(prod)) prob_post|> - kable(col.names = c("Conjetura $\\theta$", "$p(\\theta)$", "$p(D|\\theta)$", - "$p(D|\\theta)p(\\theta)$", "$p(\\theta|D)$"), escape = FALSE) |> kable_paper() - + kable(col.names = c("Conjetura θ", "p(θ)", "p(D|θ)", + "p(D|θ)p(θ)", "p(θ|D)"), escape = FALSE) |> kable_paper() ``` La última columna nos da el resultado final de la inferencia bayesiana. Podemos @@ -211,26 +295,224 @@ prob_post |> ``` -Podemos hacer esto también de manera continua. Para los datos que obtuvimos, -tenemos igualmente +## Paso 3: definir un proceso estadístico -$$p(\theta | D) \propto p(D|\theta)p(\theta) $$ -por el momento pondremos $p(\theta) = 1$ para $\theta\in [0,1]$ (densidad uniforme), entonces +Ahora podemos definir, para nuestro ejemplo discretizado, la función +que calcula la posterior dados los pasos 1 y 2: + +```{r} +calcular_posterior <- 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") |> + 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) +} +``` + +```{r} +muestra <- c(1,0,0,1,0) +``` + + +```{r} +calcular_posterior(muestra, prob_priori) +``` + +Procedemos ahora a hacer algunas pruebas simples de nuestra función: + +```{r} +calcular_posterior(rep(0, 50)) |> round(3) +calcular_posterior(rep(1, 50)) |> round(3) +calcular_posterior(c(rep(0, 100), rep(1, 100))) |> round(3) +``` + +## Paso 4: Probar el proceso de estimación + +Antes de utilizar datos, verificamos cómo se comporta nuestro proceso de estimación +de acuerdo a los supuestos de nuestro modelo generativo. + +::: callout-note +# Verificación a priori + +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: -$p(\theta|D) \propto \theta^2(1-\theta^2)$ +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) + +```{r} +theta <- 0.2 +N <- 30 +# simular +set.seed(9914) +datos_sim <- sim_pos_neg(theta = theta, N = N) +posterior <- calcular_posterior(datos_sim$Pos) +ggplot(posterior, aes(x = theta, y = prob_posterior)) + + geom_col() + + labs(x = "theta", y = "Prob posterior") + + geom_vline(xintercept = theta, color = "red", linetype = "dashed") +``` +En este caso, la estimación parece correcta. Podemo repetir el proceso +con distintos valores de $\theta$: + +```{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(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) +``` +Y vemos que en general nuestro método parece funcionar correctamente + +::: callout-tip +# Observaciones + +- Más adelante veremos cómo comparar valores a estimar con la posterior a +través de varias simulaciones de manera más rigurosa. Por el momento, recuerda +que **incluso pruebas simples o limitadas son mejores que ninguna prueba**. + +- Típicamente los valores iniciales se toman de la distribución a priori, como +hicimos arriba. Esta +prueba es en general más apropiada, pues no nos interesan configuración de parámetros +con probabilidad inicial extremadamente baja (imposibles según nuestros supuestos), pero +también es posible tomar algunos valores fijos de interés. + +::: + +Este paso también es importante para entender si, bajo nuestros propios supuestos, +es factible obtener información útil bajo el diseño que propongamos. Por ejemplo, +alguien podría proponer un diseño de muestra que sólo tome 5 personas. Podemos +probar cómo se comportan nuestras estimaciones: + +```{r} +simulacion_rep <- map_df(1:20, + function(rep){ + 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 = 3) + posterior <- calcular_posterior(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) +``` +Nuestra respuesta en este caso es que quizá con 3 personas la información obtenida +no será suficiente para tomar decisiones útiles: nótese que la posterior +está muy poco concentrada alrededor del verdadero valor de $\theta$. + + +## Paso 5: Analizar los datos y resumir resultados. + +Con este trabajo hecho (ojo: para modelos grandes es un trabajo considerable, +pero importante), podemos proceder a analizar los datos. + +Supongamos que se tomó una muestra de $N=20$ personas, con 17 negativos y +3 positivos. Calculamos la posterior: + +```{r} +# en nuestro modelo *no* importa el orden, verifica: +datos_tbl <- tibble(Pos = c(rep(1, 3), rep(0, 17))) +posterior <- calcular_posterior(muestra = datos_tbl$Pos) +ggplot(posterior, aes(x = theta, y = prob_posterior)) + + geom_col() + + labs(x = "theta", y = "Prob posterior") +``` +Y hay varias maneras de resumir esta posterior. Por ejemplo, podemos +calcular (ojo: veremos más detalles de esto más adelante): + +```{r} +# Media +posterior |> + mutate(theta = theta) |> + summarise(media = sum(theta * prob_posterior)) +# Intervalo de alta probabilidad 90% +posterior |> + mutate(theta = theta) |> + arrange(desc(prob_posterior)) |> + mutate(cumsum = cumsum(prob_posterior)) |> + filter(cumsum <= 0.9) |> + pull(theta) |> + range() +``` + +## Paso 6: Evaluar el modelo y cómputos + +En este ejemplo, el modelo es muy simple, y los cómputos son sencillos. Para modelos +más complejos es necesario checar que los cómputos sean correctos, y que el modelo +ajusta razonablemente bien a los datos en los aspectos que nos interesan, de modo +que *dejaremos esta discusión cuando veamos el flujo bayesiano más avanzado*. + + +## Versión continua + +En el ejemplo anterior utilizamos una variable aleatoria discreta para modelar +la seroprevalencia, pero esto generalmente no es conveniente. Ahora repetimos el +ejercicio considerando más naturalmente que $\theta$ puede tomar cualquier +valor en $[0,1]$. + +Para el paso 1 y 2 (definir modelo generativo y cantidad a estimar), utilizamos el mismo diagrama de arriba y la misma función que simula datos. Igual que antes, +para cualquier muestra $D$ compuesta de 0 y 1's (negativos y positivos), la probabilidad +de observar la muestra $D$ dada una conjetura $\theta$ es: + +$$ p(D|\theta) = \theta^{N_+}(1-\theta)^{N_-}$$ +Y recordamos que desde el punto de vista bayesiano, queremos resumir nuestra +información obtenida +con la distribución posterior $p(\theta|D)$, e igual que antes tenemos que: + +$$p(\theta | D) \propto p(D|\theta)p(\theta).$$ +Por el momento pondremos la densidad continua uniforme $p(\theta) = 1$ para $\theta\in [0,1]$ (densidad uniforme), entonces + +$$p(\theta|D) \propto \theta^{N_+}(1-\theta)^{N_-}$$ En este caso, para normalizar tenemos que hacer la integral de la expresión de la derecha, y dividir por el resultado. En general, escribiremos $$B(a,b) = \int_{0}^1 \theta^{a-1}(1-\theta)^{b-1} d\theta$$ -de modo que en nuestro caso +así que en nuestro caso, la posterior es: -$$p(\theta|D) = \frac{1}{B(3,4)} \theta^2(1-\theta)^3$$ +$$p(\theta|D) = \frac{1}{B(N_+ + 1,N_{-}+1)} \theta^{N_+}(1-\theta)^{N_-}$$ Es posible demostrar con cálculo que $B(a,b) = \frac{(a-1)!(b-1)!}{(a+b-1)!}$, -pero eso no es importante. Este tipo de densidades pertenecen a la familia -beta con parámetros $(a,b)$, donde $a>0, b>0$. +pero eso no es importante. Este tipo de densidades pertenecen a la **familia +beta** con parámetros $(a,b)$, donde $a>0, b>0$. -De modo que nuestra posterior es una beta con parámetros $(3,4)$, y se ve así: +Por ejemplo, si observamos 2 positivos y tres negativos, +nuestra posterior es una beta con parámetros $(3,4)$, y se ve así: ```{r} library(tidyverse) @@ -240,10 +522,23 @@ tibble(theta = theta, densidad = dbeta(theta, 3, 4)) |> geom_line() + labs(x = "theta", y = "Densidad posterior") ``` -### Ejercicio +Notamos adicionalmente que es posible seleccionar otra distribución inicial +que no sea la uniforme. En este caso particular es *conveniente* (aunque no +siempre tiene sentido) usar una distribución beta, de manera que es fácil +ver que si ponemos por ejemplo + +$$p(\theta) \propto \theta^{a}(1-\theta)^{b}$$ + +entonces la posterior, por la fórmula de Bayes, es: + +$$p(\theta|D) \propto \theta^{N_+ +a }(1-\theta)^{N_{-}+b}$$ +que también es de la familia beta, pero con parámetros +$(N_+ +a+1, N_- +b+1)$. + +### Ejercicio: actualizaciones de posterior Podemos examinar la posterior para dados distintos datos. Supondremos que -la a priori es uniforme. +la distribución a priori es uniforme. ```{r} set.seed(92192) @@ -264,23 +559,52 @@ ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) + facet_wrap(~ muestra) + geom_abline(slope = 0, intercept = 1, color = "gray") ``` -- En esta gráfica vemos que - - - - +Ahora repetimos con una inicial beta $(0,2)$ (que equivale a observar dos negativos +y ningún positivo en una muestra de 3 personas), de modo que +$p(\theta) = 2(1-\theta)$: +```{r} +set.seed(92192) +theta_seq <- seq(0,1, 0.001) +datos_sim <- sim_pos_neg(theta = 0.25, N = 12) |> + mutate(obs = ifelse(Pos==1, "P", "N")) |> + mutate(n = 1:12) +# graficar posteriores para cada n +datos_graf <- datos_sim |> + mutate(n_pos = cumsum(Pos), n_neg = cumsum(Neg)) |> + mutate(muestra = accumulate(obs, ~ paste0(.x, .y))) |> + group_by(n) |> + mutate(dens_graf = + list(tibble(theta = theta_seq, + densidad = dbeta(theta_seq, 1 + n_pos + 0, 1 + n_neg + 2)))) |> + unnest(dens_graf) +ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) + + geom_line() + + facet_wrap(~ muestra) + + geom_abline(slope = -2, intercept = 2, color = "gray") +``` +--- +En este punto, podríamos pasar al siguiente paso, que es +escribir una función para calcular la posterior. En realidad ya +sabemos su función de densidad, pero cualquier resumen que hagamos de esta +distribución requerirá de integrales (¿por qué? piensa en cómo calcular la probabilidad +de ser menor que un valor, o cómo se calcula la media). +Aunque en este ejemplo simple la posterior +tiene una forma conocida y hay manera de calcular (analíticamente o con rutinas numéricas +ya implementadas) esos resúmenes de interés +(media, cuantiles, etc.), en general calcular integrales no es una estrategia +que podamos llevar muy lejos. -### Simulación de la posterior +### Métodos Monte Carlo Una vez que tenemos la densidad posterior podemos mostrarla o resumirla de varias maneras. Si tenemos una expresión analítica, esos resúmen típicamente @@ -293,51 +617,140 @@ varias maneras. Si tenemos una expresión analítica, esos resúmen típicamente Este proceso puede ser no trivial incluso para densidades posteriores conocidas. La alternativa a integrar es **simular** de la posterior y calcular las cantidades de interés a partir de las simulaciones. En general, esto es más fácil que integrar. -En nuestro ejemplo: +En nuestro ejemplo, en lugar de usar una función de *calcular_posterior*, +construimos una que es *simular_posterior*. + +Esta función será simple porque simular de una beta es un problema estándar, +y existen muchas implementaciones. Podríamos escribir, por ejemplo: + +```{r} +simular_posterior <- function(muestra, n){ + tibble(theta = + rbeta(n, sum(muestra) + 1, length(muestra) - sum(muestra) + 1)) +} +``` + + +```{r} +muestra +sims_post <- simular_posterior(muestra, 10000) +``` +```{r} +sims_post |> + ggplot(aes(x = theta)) + + geom_histogram(bins = 50) +``` +Si queremos calcular la media, por ejemplo, hacemos + +```{r} + sims_post |> pull(theta) |> mean() +``` + +Si queremos la probabilidad de que la prevalencia esté por debajo de 20% hacemos: + +```{r} +sims_post |> + summarise(prob = mean(theta < 0.2)) +``` + +::: callout-tip +# Métodos Monte Carlo +Los métodos Monte Carlo están basados en simulación de variables aleatorias. +Las cantidades que nos interesan son integrales bajo una densidad de probabilidad. +Si queremos calcular en general +$$I = \int f(x)p(x)dx,$$ +simulamos una gran cantidad de observaciones $x_1,\ldots, x_M$ bajo $p(x)$, +y entonces (Ley de los grandes números): +$$\frac{1}{M} \sum_{i=1}^{M} x_i \to I$$ +cuando $M\to \infty$. De este modo, podemos aproximar con la precisión que +requiramos la integral $I$. +::: +**Nota importante**: +Notamos, sin embargo, que sin más información del proceso de simulación, +no es posible *demostrar* que una aproximación es "suficientemente" buena, no +importa que tan grande sea $M$. Más +adelante veremos una batería de diagnósticos para al menos excluir los casos comunes +en los que la aproximación es mala. +Finalmente, checamos todo nuestra construcción de estimación como hicimos +arriba, la diferencia es que ahora usamos simulaciones para entender el comportamiento +de la posterior. En este caso, el proceso es como sigue: +1. Generamos un valor de la apriori $\theta_{sim} \sim \text{Beta}(1,3)$ +2. Simulamos datos de la muestra ($N=25$) con el valor simulado de $\theta$ +3. Simulamos 10000 valores de la posterior +4. Repetimos 1-3 -Esta función no es una densidad de probabilidad sobre $p$, pues no integra a uno. -Sin embargo, podemos normalizarla para que integre a uno. Usando cálculo, podemos -mostrar que la **probabilidad posterior** para $p$ está dada por la densidad -$$ \frac{1}{B(E,F)} p^E(1-p)^F$$ -donde $B(E,F) = \frac{(E+F+1)!}{E!F!}$. Esta densidad se llama **beta** con -parámetros $E+1$ y $F+1$. ```{r} -simular_posterior <- function(muestra, n){ - rbeta(n, sum(muestra$Pos) + 1, sum(muestra$Neg) + 1) -} +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 = 25) + # 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} -una_muestra <- sim_pos_neg(theta = 0.1, N = 100) -simular_posterior(una_muestra, 1000) |> - tibble(theta = _) |> - ggplot(aes(x = theta)) + - geom_histogram() +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) ``` +## Observaciones +El proceso de arriba lo refinaremos considerablemente en el resto del curso. +- En primer lugar, los modelos generativos serán más complicados, y estarán +basados en teoría más compleja (que expresamos con diagramas causales) +- Usaremos más herramientas y componentes +para construir modelos estadísticos apropiados, ya +sea que construyamos un modelo completo para todo el proceso de generación de datos, +o que usemos modelos estándar como regresión para aproximar respuestas, cuando +es apropiado +- Refinaremos el proceso de checar que el cómputo (checar Monte Carlo) y +la inferencia (verificación apriori) es correcta bajo nuestros supuestos. +- Finalmente, veremos qué hacer después de hacer la estimación y que los puntos de +arriba están resueltos, para tener confianza en nuestras conclusiones. -## Modelo generativo: errores de medición 1 +## Modelo generativo: errores de medición -Ahora supongamos que la prueba no es perfecta, y que tiene una sensibilidad y -una especificidad conocida. -el diagrama es ahora: +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% @@ -347,40 +760,47 @@ digraph { graph [ranksep = 0.3, rankdir = LR] node [shape=circle] p + Npos node [shape=plaintext] N Npos [label = +>] - Nneg [label = ->] + Nobs + #Nneg [label = ->] #sens #esp edge [minlen = 3] p -> Npos - p -> Nneg + #p -> Nneg N -> Npos - N -> Nneg - esp -> Npos - sens -> Npos - esp -> Nneg - sens -> Nneg + Npos -> Nobs + #N -> Nneg + esp -> Nobs + sens -> Nobs + #esp -> Nneg + #sens -> Nneg { rank = same; p; N } -{ rank = same; Npos; Nneg} +{ rank = same; Npos} { rank = max; sens; esp} } ")#, width = 200, height = 50) ``` -#Y para constuir el modelo generativo notamos que -#la probabilidad de que un individuo sea positivo en la muestra sale -#positivo es: +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. -#$$P(E) = P(E|Pos)P(Pos) + P(E|Neg)P(Neg) = sens*p + (1-esp)*(1-p)$$ +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: -Y el modelo generativo es: ```{r} -sim_pos_neg <- function(p = 0.01, N = 20, sens = 1, esp = 1) { +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, p) + 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 @@ -396,30 +816,192 @@ Hacemos unas pruebas: ```{r} set.seed(8212) -sim_pos_neg(p = 0.3, N = 1e7, sens = 0.9, esp = 0.99) |> pull(Pos) |> mean() |> +sim_pos_neg(theta = 0.3, N = 1e7, sens = 0.7, esp = 1) |> pull(Pos) |> mean() |> round(4) -sim_pos_neg(p = 0.3, N = 1e7, sens = 1, esp = 1) |> pull(Pos) |> mean() |> +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 <- function(muestra, n, sens = 1, esp = 1){ - p_obs <- rbeta(n, sum(muestra$Pos) + 1, sum(muestra$Neg) + 1) - p <- (p_obs + esp -1 ) / (sens + esp - 1) - p +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(88) -una_muestra <- sim_pos_neg(p = 0.3, N = 600, sens = 0.5, esp = 0.99) +set.seed(328) +una_muestra <- sim_pos_neg(theta = 0.2, N = 600, sens = 0.6, esp = 0.999) mean(una_muestra$Pos) -simular_posterior(una_muestra, 1000, sens = 0.5, esp = 0.99) |> - tibble(p = _) |> - ggplot(aes(x = p)) + +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 @@ -429,9 +1011,6 @@ 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} - -```{r} ```{r} #| out-width: 100% @@ -440,24 +1019,26 @@ grViz(" digraph { graph [ranksep = 0.3, rankdir = LR] node [shape=circle] - p + theta esp sens node [shape=plaintext] N Npos [label = +>] - Nneg [label = ->] + # Nneg [label = ->] edge [minlen = 3] - p -> Npos - p -> Nneg + theta -> Npos + #p -> Nneg N -> Npos - N -> Nneg + #N -> Nneg esp -> Npos sens -> Npos - esp -> Nneg - sens -> Nneg -{ rank = same; p; N } -{ rank = same; Npos; Nneg} + #esp -> Nneg + #sens -> Nneg + esp -> Cal + sens -> Cal +{ rank = same; theta; N } +#{ rank = same; Npos; Nneg} { rank = max; sens; esp} } ")#, width = 200, height = 50) @@ -484,29 +1065,79 @@ sim_pos_neg <- function(p = 0.01, N = 20, pos_gold = c(103,122), neg_gold = c(39 } ``` +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} -simular_posterior <- function(muestra, n, pos_gold = c(103,122), neg_gold = c(399,401)){ - # simulamos sens y especificidad - sens <- rbeta(n, pos_gold[1] + 1, pos_gold[2] - pos_gold[1] + 1) - esp <- rbeta(n, neg_gold[1] + 1, neg_gold[2] - neg_gold[1] + 1) - # simulamos muestra - p_obs <- rbeta(n, sum(muestra$Pos) + 1, sum(muestra$Neg) + 1) - p <- (p_obs + esp -1 ) / (sens + esp - 1) - ifelse(p > 0, p, 0) -} +library(cmdstanr) +mod_sc <- cmdstan_model("./src/sclara.stan") +print(mod_sc) ``` + ```{r} -set.seed(828) -una_muestra <- sim_pos_neg(p = 0.01, N = 3300) -mean(una_muestra$Pos) -simular_posterior(una_muestra, 1000) |> - tibble(p = _) |> - ggplot(aes(x = p)) + +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/src/golf-logistico.stan b/notas/src/golf-logistico.stan new file mode 100644 index 0000000..1146315 --- /dev/null +++ b/notas/src/golf-logistico.stan @@ -0,0 +1,15 @@ +data { + int N; + array[N] int n; + vector[N] d; + array[N] int y; +} +parameters { + real alpha; + real beta; +} +model { + y ~ binomial_logit(n, alpha + beta * d); + alpha ~ normal(5, 2); + beta ~ normal(0, 0.02); +} diff --git a/notas/src/sclara.stan b/notas/src/sclara.stan new file mode 100644 index 0000000..51ddb27 --- /dev/null +++ b/notas/src/sclara.stan @@ -0,0 +1,32 @@ +data { + int N; + int n; + int kit_pos; + int n_kit_pos; + int kit_neg; + int n_kit_neg; +} + +parameters { + real p; //seroprevalencia + real sens; //sensibilidad + real esp; //especificidad +} + +transformed parameters { + real prob_pos; + + prob_pos = p * sens + (1 - p) * (1 - esp); + +} +model { + // modelo de número de positivos + n ~ binomial(N, prob_pos); + // modelos para resultados del kit + kit_pos ~ binomial(n_kit_pos, sens); + kit_neg ~ binomial(n_kit_neg, esp); + // iniciales para cantidades no medidas + p ~ beta(1.0, 10.0); + sens ~ beta(2.0, 1.0); + esp ~ beta(2.0, 1.0); +}