From 64cbc8e5927b6082807f51d1fc9c3adc384a0fe7 Mon Sep 17 00:00:00 2001 From: Felipe Gonzalez Date: Tue, 16 Jan 2024 20:04:45 -0600 Subject: [PATCH] Correcciones menores --- notas/02-flujo-basico-2.qmd | 31 ++++++++++++------------------- notas/02-flujo-basico.qmd | 16 +++++++++++++--- notas/src/sclara.stan | 6 +++--- 3 files changed, 28 insertions(+), 25 deletions(-) diff --git a/notas/02-flujo-basico-2.qmd b/notas/02-flujo-basico-2.qmd index e192be9..af68b9e 100644 --- a/notas/02-flujo-basico-2.qmd +++ b/notas/02-flujo-basico-2.qmd @@ -13,7 +13,7 @@ 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 +El primer aspecto importante es persistir **siguiendo 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 @@ -48,7 +48,7 @@ digraph { node [shape=plaintext] N Npos [label = +>] - Nobs + Nobs [label = obs>] #Nneg [label = ->] #sens #esp @@ -126,7 +126,8 @@ 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 +de 0's y 1's, tenemos que la verosimilitud de la muestra dada cada +conjetura $\theta$, y con $sens$ y $esp$ fijas, es: $$p(D|\theta, sens, esp) = \theta_{obs}^{N_{+}}(1-\theta_{obs})^{N_{-}}$$ Suponiendo que la distribución apriori de $\theta$ es uniforme, tenemos entonces @@ -215,6 +216,8 @@ Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación: ```{r} +#| label: fig-apriori-falla +#| fig-cap: "Verificación a priori" ggplot(simulacion_rep_error, aes(x = theta)) + geom_histogram(bins = 50) + labs(x = "theta", y = "Prob posterior") + @@ -245,6 +248,8 @@ Ahora usamos histogramas por ejemplo para mostrar cómo luce la posterior, y comparamos con los valores de la simulación: ```{r} +#| label: fig-apriori-falla +#| fig-cap: "Verificación a priori fallida (modelo incorrecto)" ggplot(simulacion_rep, aes(x = theta)) + geom_histogram(bins = 50) + labs(x = "theta", y = "Prob posterior") + @@ -378,8 +383,8 @@ 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")) +sims <- ajuste$draws(c("theta", "sens", "esp"), format = "df") +resumen <- ajuste$summary(c("theta")) ``` ```{r} @@ -389,7 +394,7 @@ 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)) + +ggplot(sims, aes(x = theta)) + geom_histogram() ``` @@ -413,18 +418,6 @@ para explicar la prevalencia observada, pero si no pensamos con cuidado podríam 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 +el rango de $\theta$ 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 e747591..ced24af 100644 --- a/notas/02-flujo-basico.qmd +++ b/notas/02-flujo-basico.qmd @@ -701,7 +701,7 @@ ggplot(datos_graf, aes(x=theta, y = densidad, group = n)) + --- -En este punto, podríamos pasar al siguiente paso, que es +En este punto, podríamos ir 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 @@ -784,13 +784,23 @@ 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, +**Nota 1**: 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. +**Nota 2**: En nuestro caso, las integrales de interés usualmente son de +la forma +$$I = \int f(\theta)p(\theta|D) d\theta,$$ +donde $D$ es la información de la muestra, $\theta$ es un vector de parámetros +del modelo, y $f(\theta)$ es una función de $\theta$ que nos interesa. Por ejemplo, +para la media posterior de $\theta_i$, usaríamos $f(\theta) = \theta_i$. Podemos +aproximar cualquier integral si tenemos simulaciones de la posterior: + +$$\theta_i \sim p(\theta|D) \implies \frac{1}{M} \sum_{i=1}^{M} f(\theta_i) \to I.$$ + +--- Finalmente, checamos todo nuestra construcción de estimación como hicimos arriba, la diferencia es que ahora usamos simulaciones para entender el comportamiento diff --git a/notas/src/sclara.stan b/notas/src/sclara.stan index 51ddb27..812224f 100644 --- a/notas/src/sclara.stan +++ b/notas/src/sclara.stan @@ -8,7 +8,7 @@ data { } parameters { - real p; //seroprevalencia + real theta; //seroprevalencia real sens; //sensibilidad real esp; //especificidad } @@ -16,7 +16,7 @@ parameters { transformed parameters { real prob_pos; - prob_pos = p * sens + (1 - p) * (1 - esp); + prob_pos = theta * sens + (1 - theta) * (1 - esp); } model { @@ -26,7 +26,7 @@ model { 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); + theta ~ beta(1.0, 10.0); sens ~ beta(2.0, 1.0); esp ~ beta(2.0, 1.0); }